diff --git a/src/main/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2.java b/src/main/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2.java index 707eea94..28b9a9ee 100644 --- a/src/main/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2.java +++ b/src/main/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2.java @@ -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); } } diff --git a/src/main/java/org/jtransforms/fft/ReusableFloatInverseFFT.java b/src/main/java/org/jtransforms/fft/ReusableFloatInverseFFT.java new file mode 100644 index 00000000..abafce1e --- /dev/null +++ b/src/main/java/org/jtransforms/fft/ReusableFloatInverseFFT.java @@ -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 + * **************************************************************************** + */ +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 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; + } + } +} diff --git a/src/test/java/org/jtransforms/fft/ReusableFloatInverseFFTTest.java b/src/test/java/org/jtransforms/fft/ReusableFloatInverseFFTTest.java new file mode 100644 index 00000000..44b8a847 --- /dev/null +++ b/src/test/java/org/jtransforms/fft/ReusableFloatInverseFFTTest.java @@ -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 + * **************************************************************************** + */ +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"); + } + } +}