Extends JTransforms to reduce memory allocation on a per FFT iteration.

This commit is contained in:
denny
2026-11-02 04:35:18 -05:00
parent 0b3ce66cad
commit 415a7017f9
3 changed files with 294 additions and 4 deletions

View File

@@ -27,7 +27,7 @@ import java.util.Arrays;
import java.util.Queue;
import java.util.concurrent.ArrayBlockingQueue;
import org.apache.commons.math3.util.FastMath;
import org.jtransforms.fft.FloatFFT_1D;
import org.jtransforms.fft.ReusableFloatInverseFFT;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
@@ -85,7 +85,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
//Sized to process 40 times per second
private IFFTProcessorDispatcher mIFFTProcessorDispatcher = new IFFTProcessorDispatcher(25);
private FloatFFT_1D mFFT;
private ReusableFloatInverseFFT mFFT;
private ChannelDFT mChannelDFT;
private float[] mInlineSamples;
private float[] mInlineFilter;
@@ -453,7 +453,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
private void init(float[] coefficients)
{
int channelCount = getChannelCount();
mFFT = new FloatFFT_1D(channelCount);
mFFT = new ReusableFloatInverseFFT(channelCount);
mChannelDFT = ChannelDFTFactory.getChannelDFT(channelCount);
int bufferLength = getSubChannelCount() * mTapsPerChannel;
mSamplesPerBlock = channelCount; //Same as subChannelCount / 2
@@ -505,7 +505,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
for(int offset = 0; offset < samples.length; offset += blockLength)
{
//Rotate each of the channels to the correct phase using the IFFT
mFFT.complexInverse(samples, offset, true);
mFFT.complexInverse(samples, offset);
}
}

View File

@@ -0,0 +1,156 @@
/*
* *****************************************************************************
* Copyright (C) 2014-2026 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
* ****************************************************************************
*/
package org.jtransforms.fft;
import java.util.ArrayList;
import java.util.List;
/**
* Scaled complex inverse FFT that produces the same values as FloatFFT_1D.complexInverse(samples, offset, true) without
* allocating on each call.
*
* FloatFFT_1D allocates a scratch array of two floats per channel on every mixed-radix transform. This class runs the
* same mixed-radix passes using a single reused scratch array. It's in this package because those passes are package
* private. Sizes that FloatFFT_1D doesn't transform as mixed-radix (powers of two and sizes with other prime factors)
* are delegated to FloatFFT_1D.
*
* Note: this class isn't thread safe because the scratch array is shared across calls.
*/
public class ReusableFloatInverseFFT
{
private static final int[] SUPPORTED_RADIX = {3, 5};
private final FloatFFT_1D mFFT;
private final int mSize;
private final int[] mFactors;
private final float[] mScratch;
private final float mScale;
/**
* Constructs an instance
* @param size of the transform in complex samples
*/
public ReusableFloatInverseFFT(int size)
{
mSize = size;
mFFT = new FloatFFT_1D(size);
mFactors = getFactors(size);
mScratch = mFactors != null ? new float[2 * size] : null;
mScale = 1.0f / size;
}
/**
* Identifies the order of radix passes for the size, matching the order used by FloatFFT_1D.
* @param size of the transform
* @return radix values in pass order, or null if the size should be delegated to FloatFFT_1D.
*/
private static int[] getFactors(int size)
{
if(size < 2 || Integer.bitCount(size) == 1)
{
return null;
}
List<Integer> factors = new ArrayList<>();
int remaining = size;
while(remaining % 4 == 0)
{
factors.add(4);
remaining /= 4;
}
//FloatFFT_1D moves a radix-2 pass to the front
if(remaining % 2 == 0)
{
factors.addFirst(2);
remaining /= 2;
}
for(int radix : SUPPORTED_RADIX)
{
while(remaining % radix == 0)
{
factors.add(radix);
remaining /= radix;
}
}
if(remaining != 1)
{
return null;
}
return factors.stream().mapToInt(Integer::intValue).toArray();
}
/**
* Performs a scaled inverse transform in place.
* @param samples containing interleaved I/Q values
* @param offset to the start of the block in the samples array
*/
public void complexInverse(float[] samples, int offset)
{
if(mFactors == null)
{
mFFT.complexInverse(samples, offset, true);
return;
}
final float[] scratch = mScratch;
final int size = mSize;
boolean resultsInScratch = false;
int subTransforms = 1;
int twiddleOffset = 2 * size;
for(int radix : mFactors)
{
int next = radix * subTransforms;
int stride = 2 * (size / next);
float[] in = resultsInScratch ? scratch : samples;
int inOffset = resultsInScratch ? 0 : offset;
float[] out = resultsInScratch ? samples : scratch;
int outOffset = resultsInScratch ? offset : 0;
switch(radix)
{
case 2 -> mFFT.passf2(stride, subTransforms, in, inOffset, out, outOffset, twiddleOffset, 1);
case 3 -> mFFT.passf3(stride, subTransforms, in, inOffset, out, outOffset, twiddleOffset, 1);
case 4 -> mFFT.passf4(stride, subTransforms, in, inOffset, out, outOffset, twiddleOffset, 1);
default -> mFFT.passf5(stride, subTransforms, in, inOffset, out, outOffset, twiddleOffset, 1);
}
resultsInScratch = !resultsInScratch;
subTransforms = next;
twiddleOffset += (radix - 1) * stride;
}
if(resultsInScratch)
{
System.arraycopy(scratch, 0, samples, offset, 2 * size);
}
final float scale = mScale;
final int end = offset + 2 * size;
for(int x = offset; x < end; x++)
{
samples[x] *= scale;
}
}
}

View File

@@ -0,0 +1,134 @@
/*
* *****************************************************************************
* Copyright (C) 2014-2026 Dennis Sheirer
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>
* ****************************************************************************
*/
package org.jtransforms.fft;
import static org.junit.jupiter.api.Assertions.assertArrayEquals;
import static org.junit.jupiter.api.Assertions.assertEquals;
import java.lang.management.ManagementFactory;
import java.util.Random;
import org.junit.jupiter.api.Test;
/**
* Tests that the reusable inverse FFT is an exact replacement for FloatFFT_1D.
*/
class ReusableFloatInverseFFTTest
{
private static final int OFFSET = 6;
private static final int LARGEST_SIZE = 5000;
private static float[] noise(int length, Random random)
{
float[] samples = new float[length];
for(int x = 0; x < length; x++)
{
samples[x] = random.nextFloat() * 2.0f - 1.0f;
}
return samples;
}
/**
* Test: output values are bit-identical to FloatFFT_1D for every size, including powers of two and sizes with other
* prime factors, and samples outside of the block aren't changed. Each instance is used twice to show that
* scratch contents don't carry over between calls.
*/
@Test
void matchesFloatFFT1D()
{
Random random = new Random(0xC0FFEEL);
for(int size = 1; size <= LARGEST_SIZE; size += size < 512 ? 1 : 53)
{
ReusableFloatInverseFFT reusable = new ReusableFloatInverseFFT(size);
FloatFFT_1D reference = new FloatFFT_1D(size);
for(int pass = 0; pass < 2; pass++)
{
float[] expected = noise(2 * size + 2 * OFFSET, random);
float[] actual = expected.clone();
reference.complexInverse(expected, OFFSET, true);
reusable.complexInverse(actual, OFFSET);
assertArrayEquals(expected, actual, 0.0f, "Size " + size);
}
}
}
/**
* Test: every channel count from 2 to 4800 with only factors of 2, 3 and 5 matches FloatFFT_1D. These are the only
* sizes the channelizer can select.
*/
@Test
void matchesFloatFFT1DForChannelizerSizes()
{
Random random = new Random(0xBEEFL);
for(int size = 2; size <= 4800; size += 2)
{
int remaining = size;
for(int factor : new int[]{2, 3, 5})
{
while(remaining % factor == 0)
{
remaining /= factor;
}
}
if(remaining == 1)
{
float[] expected = noise(2 * size, random);
float[] actual = expected.clone();
new FloatFFT_1D(size).complexInverse(expected, 0, true);
new ReusableFloatInverseFFT(size).complexInverse(actual, 0);
assertArrayEquals(expected, actual, 0.0f, "Size " + size);
}
}
}
/**
* Test: transforms don't allocate, for a mixed-radix size with an odd pass count and an even pass count.
*/
@Test
void doesNotAllocate()
{
com.sun.management.ThreadMXBean threads = (com.sun.management.ThreadMXBean)ManagementFactory.getThreadMXBean();
long threadId = Thread.currentThread().threadId();
for(int size : new int[]{96, 400})
{
ReusableFloatInverseFFT fft = new ReusableFloatInverseFFT(size);
float[] samples = noise(2 * size, new Random(1L));
fft.complexInverse(samples, 0);
long before = threads.getThreadAllocatedBytes(threadId);
for(int x = 0; x < 1000; x++)
{
fft.complexInverse(samples, 0);
}
long allocated = threads.getThreadAllocatedBytes(threadId) - before;
//FloatFFT_1D would allocate at least 8 bytes per channel, per call. Allow for measurement overhead.
assertEquals(0, allocated / (1000L * size), "Size " + size + " allocated " + allocated + " bytes");
}
}
}