From 1c097491d96b29ebb3a83bae8c276be820f192c7 Mon Sep 17 00:00:00 2001 From: Denny Sheirer Date: Tue, 3 Nov 2026 05:49:40 -0500 Subject: [PATCH] #2567 Reduces memory overhead by extending the JTransforms FFT package to reuse its float array when calculaing the inverse FFT in the channelizer for a subset of the FFT size operations. (#2568) Co-authored-by: sheirerd --- .../dsp/fft/MixedRadixComplexFFT.java | 779 ++++++++++++++++++ .../ComplexPolyphaseChannelizerM2.java | 9 +- .../dsp/fft/MixedRadixComplexFFTTest.java | 286 +++++++ .../ComplexPolyphaseChannelizerM2Test.java | 8 +- 4 files changed, 1075 insertions(+), 7 deletions(-) create mode 100644 src/main/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFT.java create mode 100644 src/test/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFTTest.java diff --git a/src/main/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFT.java b/src/main/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFT.java new file mode 100644 index 00000000..38213e19 --- /dev/null +++ b/src/main/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFT.java @@ -0,0 +1,779 @@ +/* + * ***************************************************************************** + * 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 + * **************************************************************************** + * + * The transform kernels in this class are ported from the mixed radix complex FFT (cffti, cfftf and passf2/3/4/5) + * in JTransforms org.jtransforms.fft.FloatFFT_1D, which carries the following license: + * + * JTransforms + * Copyright (c) 2007 onward, Piotr Wendykier + * All rights reserved. + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, this + * list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR + * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ +package io.github.dsheirer.dsp.fft; + +/** + * Allocation-free mixed radix complex FFT for sizes that factor into 2, 3 and 5 only. + * + * This is a port of the JTransforms FloatFFT_1D mixed radix complex transform path. JTransforms allocates a new + * 2N float scratch array on every mixed radix transform, which creates significant garbage collection pressure when + * the transform runs at block rate (e.g. the polyphase channelizer IFFT). This implementation pre-allocates the + * scratch buffer once and reuses it for each transform. + * + * Results match JTransforms FloatFFT_1D complexForward() and complexInverse(): the forward transform uses the + * exp(-j...) kernel, the inverse uses exp(+j...), and the inverse optionally scales by 1/N. Twiddle factors are + * calculated in double precision, so results can differ from JTransforms in the least significant bits. + * + * Samples are interleaved complex: a[offset + 2k] = Re[k] and a[offset + 2k + 1] = Im[k] for 0 <= k < N. + * + * Note: this class is NOT thread-safe since it reuses an internal scratch buffer. Use one instance per thread. + */ +public class MixedRadixComplexFFT +{ + private static final int[] FACTORS = {4, 2, 3, 5}; + private static final int FORWARD = -1; + private static final int INVERSE = 1; + + private final int mSize; + private final int[] mFactors; + private final float[] mTwiddles; + private final float[] mScratch; + private final float mScale; + + /** + * Constructs an instance + * @param size of the transform in complex samples. Must be a positive value with no prime factors other than + * 2, 3 and 5. + * @throws IllegalArgumentException if the size is not supported + */ + public MixedRadixComplexFFT(int size) + { + if(!isSupported(size)) + { + throw new IllegalArgumentException("Unsupported FFT size [" + size + + "] - size must be positive and have no prime factors other than 2, 3 and 5"); + } + + mSize = size; + mFactors = factorize(size); + mTwiddles = new float[2 * size + 2]; + mScratch = new float[2 * size]; + mScale = 1.0f / size; + calculateTwiddles(); + } + + /** + * Indicates if the size can be transformed by this class. + * @param size in complex samples + * @return true if the size is positive and has no prime factors other than 2, 3 and 5. + */ + public static boolean isSupported(int size) + { + if(size < 1) + { + return false; + } + + int remainder = size; + + for(int factor: new int[]{2, 3, 5}) + { + while(remainder % factor == 0) + { + remainder /= factor; + } + } + + return remainder == 1; + } + + /** + * Size of this transform + * @return size in complex samples + */ + public int getSize() + { + return mSize; + } + + /** + * Computes the forward DFT of the interleaved complex samples, in place. Equivalent to JTransforms + * FloatFFT_1D.complexForward(samples, offset). + * @param samples interleaved complex, with at least 2N values from the offset + * @param offset to the first sample + */ + public void forward(float[] samples, int offset) + { + transform(samples, offset, FORWARD); + } + + /** + * Computes the forward DFT of the interleaved complex samples, in place. + * @param samples interleaved complex, with at least 2N values + */ + public void forward(float[] samples) + { + forward(samples, 0); + } + + /** + * Computes the inverse DFT of the interleaved complex samples, in place. Equivalent to JTransforms + * FloatFFT_1D.complexInverse(samples, offset, scale). + * @param samples interleaved complex, with at least 2N values from the offset + * @param offset to the first sample + * @param scale true to scale the results by 1/N + */ + public void inverse(float[] samples, int offset, boolean scale) + { + transform(samples, offset, INVERSE); + + if(scale && mSize > 1) + { + int end = offset + 2 * mSize; + + for(int x = offset; x < end; x++) + { + samples[x] *= mScale; + } + } + } + + /** + * Computes the inverse DFT of the interleaved complex samples, in place. + * @param samples interleaved complex, with at least 2N values + * @param scale true to scale the results by 1/N + */ + public void inverse(float[] samples, boolean scale) + { + inverse(samples, 0, scale); + } + + /** + * Factors the size into the radix sequence used by the transform, ported from JTransforms cffti(). Radix 4 is + * preferred, then 2, 3 and 5, and a radix 2 factor is moved to the front of the sequence. + */ + private static int[] factorize(int size) + { + if(size == 1) + { + return new int[0]; + } + + int[] factors = new int[32]; + int count = 0; + int remainder = size; + + for(int factor: FACTORS) + { + while(remainder % factor == 0) + { + factors[count++] = factor; + remainder /= factor; + + if(factor == 2 && count != 1) + { + System.arraycopy(factors, 0, factors, 1, count - 1); + factors[0] = 2; + } + } + } + + if(remainder != 1) + { + throw new IllegalArgumentException("Unsupported FFT size [" + size + "]"); + } + + int[] result = new int[count]; + System.arraycopy(factors, 0, result, 0, count); + return result; + } + + /** + * Calculates the twiddle factor table, ported from JTransforms cffti() with the table offset removed and the + * angles calculated in double precision. + */ + private void calculateTwiddles() + { + if(mSize == 1) + { + return; + } + + double angleIncrement = 2.0 * Math.PI / mSize; + int i = 1; + int l1 = 1; + + for(int factor: mFactors) + { + int ld = 0; + int l2 = l1 * factor; + int ido = mSize / l2; + int idot = ido + ido + 2; + + for(int j = 1; j < factor; j++) + { + mTwiddles[i - 1] = 1; + mTwiddles[i] = 0; + ld += l1; + int fi = 0; + double argld = ld * angleIncrement; + + for(int ii = 4; ii <= idot; ii += 2) + { + i += 2; + fi++; + double arg = fi * argld; + mTwiddles[i - 1] = (float)Math.cos(arg); + mTwiddles[i] = (float)Math.sin(arg); + } + } + + l1 = l2; + } + } + + /** + * Performs the transform, ported from JTransforms cfftf() using the reusable scratch buffer. + * @param a samples + * @param offa offset to first sample + * @param isign -1 for forward and +1 for inverse + */ + private void transform(float[] a, int offa, int isign) + { + if(mSize == 1) + { + return; + } + + float[] ch = mScratch; + int twon = 2 * mSize; + boolean inScratch = false; + int l1 = 1; + int iw = 0; + + for(int factor: mFactors) + { + int l2 = factor * l1; + int idot = 2 * (mSize / l2); + + switch(factor) + { + case 4: + if(inScratch) + { + passf4(idot, l1, ch, 0, a, offa, iw, isign); + } + else + { + passf4(idot, l1, a, offa, ch, 0, iw, isign); + } + break; + case 2: + if(inScratch) + { + passf2(idot, l1, ch, 0, a, offa, iw, isign); + } + else + { + passf2(idot, l1, a, offa, ch, 0, iw, isign); + } + break; + case 3: + if(inScratch) + { + passf3(idot, l1, ch, 0, a, offa, iw, isign); + } + else + { + passf3(idot, l1, a, offa, ch, 0, iw, isign); + } + break; + case 5: + if(inScratch) + { + passf5(idot, l1, ch, 0, a, offa, iw, isign); + } + else + { + passf5(idot, l1, a, offa, ch, 0, iw, isign); + } + break; + default: + throw new IllegalStateException("Unsupported radix: " + factor); + } + + inScratch = !inScratch; + l1 = l2; + iw += (factor - 1) * idot; + } + + if(inScratch) + { + System.arraycopy(ch, 0, a, offa, twon); + } + } + + /** + * Radix 2 butterfly, ported from JTransforms passf2(). + */ + private void passf2(final int ido, final int l1, final float[] in, final int in_off, final float[] out, + final int out_off, final int offset, final int isign) + { + final float[] wtable = mTwiddles; + int iw1 = offset; + int idx = ido * l1; + + if(ido <= 2) + { + for(int k = 0; k < l1; k++) + { + int idx0 = k * ido; + int iidx1 = in_off + 2 * idx0; + int iidx2 = iidx1 + ido; + float a1r = in[iidx1]; + float a1i = in[iidx1 + 1]; + float a2r = in[iidx2]; + float a2i = in[iidx2 + 1]; + + int oidx1 = out_off + idx0; + int oidx2 = oidx1 + idx; + out[oidx1] = a1r + a2r; + out[oidx1 + 1] = a1i + a2i; + out[oidx2] = a1r - a2r; + out[oidx2 + 1] = a1i - a2i; + } + } + else + { + for(int k = 0; k < l1; k++) + { + for(int i = 0; i < ido - 1; i += 2) + { + int idx0 = k * ido; + int iidx1 = in_off + i + 2 * idx0; + int iidx2 = iidx1 + ido; + float i1r = in[iidx1]; + float i1i = in[iidx1 + 1]; + float i2r = in[iidx2]; + float i2i = in[iidx2 + 1]; + + int widx1 = i + iw1; + float w1r = wtable[widx1]; + float w1i = isign * wtable[widx1 + 1]; + + float t1r = i1r - i2r; + float t1i = i1i - i2i; + + int oidx1 = out_off + i + idx0; + int oidx2 = oidx1 + idx; + out[oidx1] = i1r + i2r; + out[oidx1 + 1] = i1i + i2i; + out[oidx2] = w1r * t1r - w1i * t1i; + out[oidx2 + 1] = w1r * t1i + w1i * t1r; + } + } + } + } + + /** + * Radix 3 butterfly, ported from JTransforms passf3(). + */ + private void passf3(final int ido, final int l1, final float[] in, final int in_off, final float[] out, + final int out_off, final int offset, final int isign) + { + final float taur = -0.5f; + final float taui = 0.866025403784438707610604524234076962f; + final float[] wtable = mTwiddles; + int iw1 = offset; + int iw2 = iw1 + ido; + final int idxt = l1 * ido; + + if(ido == 2) + { + for(int k = 1; k <= l1; k++) + { + int iidx1 = in_off + (3 * k - 2) * ido; + int iidx2 = iidx1 + ido; + int iidx3 = iidx1 - ido; + float i1r = in[iidx1]; + float i1i = in[iidx1 + 1]; + float i2r = in[iidx2]; + float i2i = in[iidx2 + 1]; + float i3r = in[iidx3]; + float i3i = in[iidx3 + 1]; + + float tr2 = i1r + i2r; + float cr2 = i3r + taur * tr2; + float ti2 = i1i + i2i; + float ci2 = i3i + taur * ti2; + float cr3 = isign * taui * (i1r - i2r); + float ci3 = isign * taui * (i1i - i2i); + + int oidx1 = out_off + (k - 1) * ido; + int oidx2 = oidx1 + idxt; + int oidx3 = oidx2 + idxt; + out[oidx1] = i3r + tr2; + out[oidx1 + 1] = i3i + ti2; + out[oidx2] = cr2 - ci3; + out[oidx2 + 1] = ci2 + cr3; + out[oidx3] = cr2 + ci3; + out[oidx3 + 1] = ci2 - cr3; + } + } + else + { + for(int k = 1; k <= l1; k++) + { + int idx1 = in_off + (3 * k - 2) * ido; + int idx2 = out_off + (k - 1) * ido; + + for(int i = 0; i < ido - 1; i += 2) + { + int iidx1 = i + idx1; + int iidx2 = iidx1 + ido; + int iidx3 = iidx1 - ido; + float a1r = in[iidx1]; + float a1i = in[iidx1 + 1]; + float a2r = in[iidx2]; + float a2i = in[iidx2 + 1]; + float a3r = in[iidx3]; + float a3i = in[iidx3 + 1]; + + float tr2 = a1r + a2r; + float cr2 = a3r + taur * tr2; + float ti2 = a1i + a2i; + float ci2 = a3i + taur * ti2; + float cr3 = isign * taui * (a1r - a2r); + float ci3 = isign * taui * (a1i - a2i); + float dr2 = cr2 - ci3; + float dr3 = cr2 + ci3; + float di2 = ci2 + cr3; + float di3 = ci2 - cr3; + + int widx1 = i + iw1; + int widx2 = i + iw2; + float w1r = wtable[widx1]; + float w1i = isign * wtable[widx1 + 1]; + float w2r = wtable[widx2]; + float w2i = isign * wtable[widx2 + 1]; + + int oidx1 = i + idx2; + int oidx2 = oidx1 + idxt; + int oidx3 = oidx2 + idxt; + out[oidx1] = a3r + tr2; + out[oidx1 + 1] = a3i + ti2; + out[oidx2] = w1r * dr2 - w1i * di2; + out[oidx2 + 1] = w1r * di2 + w1i * dr2; + out[oidx3] = w2r * dr3 - w2i * di3; + out[oidx3 + 1] = w2r * di3 + w2i * dr3; + } + } + } + } + + /** + * Radix 4 butterfly, ported from JTransforms passf4(). + */ + private void passf4(final int ido, final int l1, final float[] in, final int in_off, final float[] out, + final int out_off, final int offset, final int isign) + { + final float[] wtable = mTwiddles; + int iw1 = offset; + int iw2 = iw1 + ido; + int iw3 = iw2 + ido; + int idx0 = l1 * ido; + + if(ido == 2) + { + for(int k = 0; k < l1; k++) + { + int idxt1 = k * ido; + int iidx1 = in_off + 4 * idxt1 + 1; + int iidx2 = iidx1 + ido; + int iidx3 = iidx2 + ido; + int iidx4 = iidx3 + ido; + + float i1i = in[iidx1 - 1]; + float i1r = in[iidx1]; + float i2i = in[iidx2 - 1]; + float i2r = in[iidx2]; + float i3i = in[iidx3 - 1]; + float i3r = in[iidx3]; + float i4i = in[iidx4 - 1]; + float i4r = in[iidx4]; + + float ti1 = i1r - i3r; + float ti2 = i1r + i3r; + float tr4 = i4r - i2r; + float ti3 = i2r + i4r; + float tr1 = i1i - i3i; + float tr2 = i1i + i3i; + float ti4 = i2i - i4i; + float tr3 = i2i + i4i; + + int oidx1 = out_off + idxt1; + int oidx2 = oidx1 + idx0; + int oidx3 = oidx2 + idx0; + int oidx4 = oidx3 + idx0; + out[oidx1] = tr2 + tr3; + out[oidx1 + 1] = ti2 + ti3; + out[oidx2] = tr1 + isign * tr4; + out[oidx2 + 1] = ti1 + isign * ti4; + out[oidx3] = tr2 - tr3; + out[oidx3 + 1] = ti2 - ti3; + out[oidx4] = tr1 - isign * tr4; + out[oidx4 + 1] = ti1 - isign * ti4; + } + } + else + { + for(int k = 0; k < l1; k++) + { + int idx1 = k * ido; + int idx2 = in_off + 1 + 4 * idx1; + + for(int i = 0; i < ido - 1; i += 2) + { + int iidx1 = i + idx2; + int iidx2 = iidx1 + ido; + int iidx3 = iidx2 + ido; + int iidx4 = iidx3 + ido; + float i1i = in[iidx1 - 1]; + float i1r = in[iidx1]; + float i2i = in[iidx2 - 1]; + float i2r = in[iidx2]; + float i3i = in[iidx3 - 1]; + float i3r = in[iidx3]; + float i4i = in[iidx4 - 1]; + float i4r = in[iidx4]; + + float ti1 = i1r - i3r; + float ti2 = i1r + i3r; + float ti3 = i2r + i4r; + float tr4 = i4r - i2r; + float tr1 = i1i - i3i; + float tr2 = i1i + i3i; + float ti4 = i2i - i4i; + float tr3 = i2i + i4i; + float cr3 = tr2 - tr3; + float ci3 = ti2 - ti3; + float cr2 = tr1 + isign * tr4; + float cr4 = tr1 - isign * tr4; + float ci2 = ti1 + isign * ti4; + float ci4 = ti1 - isign * ti4; + + int widx1 = i + iw1; + int widx2 = i + iw2; + int widx3 = i + iw3; + float w1r = wtable[widx1]; + float w1i = isign * wtable[widx1 + 1]; + float w2r = wtable[widx2]; + float w2i = isign * wtable[widx2 + 1]; + float w3r = wtable[widx3]; + float w3i = isign * wtable[widx3 + 1]; + + int oidx1 = out_off + i + idx1; + int oidx2 = oidx1 + idx0; + int oidx3 = oidx2 + idx0; + int oidx4 = oidx3 + idx0; + out[oidx1] = tr2 + tr3; + out[oidx1 + 1] = ti2 + ti3; + out[oidx2] = w1r * cr2 - w1i * ci2; + out[oidx2 + 1] = w1r * ci2 + w1i * cr2; + out[oidx3] = w2r * cr3 - w2i * ci3; + out[oidx3 + 1] = w2r * ci3 + w2i * cr3; + out[oidx4] = w3r * cr4 - w3i * ci4; + out[oidx4 + 1] = w3r * ci4 + w3i * cr4; + } + } + } + } + + /** + * Radix 5 butterfly, ported from JTransforms passf5(). + */ + private void passf5(final int ido, final int l1, final float[] in, final int in_off, final float[] out, + final int out_off, final int offset, final int isign) + { + final float tr11 = 0.309016994374947451262869435595348477f; + final float ti11 = 0.951056516295153531181938433292089030f; + final float tr12 = -0.809016994374947340240566973079694435f; + final float ti12 = 0.587785252292473248125759255344746634f; + final float[] wtable = mTwiddles; + int iw1 = offset; + int iw2 = iw1 + ido; + int iw3 = iw2 + ido; + int iw4 = iw3 + ido; + int idx0 = l1 * ido; + + if(ido == 2) + { + for(int k = 1; k <= l1; ++k) + { + int iidx1 = in_off + (5 * k - 4) * ido + 1; + int iidx2 = iidx1 + ido; + int iidx3 = iidx1 - ido; + int iidx4 = iidx2 + ido; + int iidx5 = iidx4 + ido; + + float i1i = in[iidx1 - 1]; + float i1r = in[iidx1]; + float i2i = in[iidx2 - 1]; + float i2r = in[iidx2]; + float i3i = in[iidx3 - 1]; + float i3r = in[iidx3]; + float i4i = in[iidx4 - 1]; + float i4r = in[iidx4]; + float i5i = in[iidx5 - 1]; + float i5r = in[iidx5]; + + float ti5 = i1r - i5r; + float ti2 = i1r + i5r; + float ti4 = i2r - i4r; + float ti3 = i2r + i4r; + float tr5 = i1i - i5i; + float tr2 = i1i + i5i; + float tr4 = i2i - i4i; + float tr3 = i2i + i4i; + float cr2 = i3i + tr11 * tr2 + tr12 * tr3; + float ci2 = i3r + tr11 * ti2 + tr12 * ti3; + float cr3 = i3i + tr12 * tr2 + tr11 * tr3; + float ci3 = i3r + tr12 * ti2 + tr11 * ti3; + float cr5 = isign * (ti11 * tr5 + ti12 * tr4); + float ci5 = isign * (ti11 * ti5 + ti12 * ti4); + float cr4 = isign * (ti12 * tr5 - ti11 * tr4); + float ci4 = isign * (ti12 * ti5 - ti11 * ti4); + + int oidx1 = out_off + (k - 1) * ido; + int oidx2 = oidx1 + idx0; + int oidx3 = oidx2 + idx0; + int oidx4 = oidx3 + idx0; + int oidx5 = oidx4 + idx0; + out[oidx1] = i3i + tr2 + tr3; + out[oidx1 + 1] = i3r + ti2 + ti3; + out[oidx2] = cr2 - ci5; + out[oidx2 + 1] = ci2 + cr5; + out[oidx3] = cr3 - ci4; + out[oidx3 + 1] = ci3 + cr4; + out[oidx4] = cr3 + ci4; + out[oidx4 + 1] = ci3 - cr4; + out[oidx5] = cr2 + ci5; + out[oidx5 + 1] = ci2 - cr5; + } + } + else + { + for(int k = 1; k <= l1; k++) + { + int idx1 = in_off + 1 + (k * 5 - 4) * ido; + int idx2 = out_off + (k - 1) * ido; + + for(int i = 0; i < ido - 1; i += 2) + { + int iidx1 = i + idx1; + int iidx2 = iidx1 + ido; + int iidx3 = iidx1 - ido; + int iidx4 = iidx2 + ido; + int iidx5 = iidx4 + ido; + float i1i = in[iidx1 - 1]; + float i1r = in[iidx1]; + float i2i = in[iidx2 - 1]; + float i2r = in[iidx2]; + float i3i = in[iidx3 - 1]; + float i3r = in[iidx3]; + float i4i = in[iidx4 - 1]; + float i4r = in[iidx4]; + float i5i = in[iidx5 - 1]; + float i5r = in[iidx5]; + + float ti5 = i1r - i5r; + float ti2 = i1r + i5r; + float ti4 = i2r - i4r; + float ti3 = i2r + i4r; + float tr5 = i1i - i5i; + float tr2 = i1i + i5i; + float tr4 = i2i - i4i; + float tr3 = i2i + i4i; + float cr2 = i3i + tr11 * tr2 + tr12 * tr3; + float ci2 = i3r + tr11 * ti2 + tr12 * ti3; + float cr3 = i3i + tr12 * tr2 + tr11 * tr3; + float ci3 = i3r + tr12 * ti2 + tr11 * ti3; + float cr5 = isign * (ti11 * tr5 + ti12 * tr4); + float ci5 = isign * (ti11 * ti5 + ti12 * ti4); + float cr4 = isign * (ti12 * tr5 - ti11 * tr4); + float ci4 = isign * (ti12 * ti5 - ti11 * ti4); + float dr3 = cr3 - ci4; + float dr4 = cr3 + ci4; + float di3 = ci3 + cr4; + float di4 = ci3 - cr4; + float dr5 = cr2 + ci5; + float dr2 = cr2 - ci5; + float di5 = ci2 - cr5; + float di2 = ci2 + cr5; + + int widx1 = i + iw1; + int widx2 = i + iw2; + int widx3 = i + iw3; + int widx4 = i + iw4; + float w1r = wtable[widx1]; + float w1i = isign * wtable[widx1 + 1]; + float w2r = wtable[widx2]; + float w2i = isign * wtable[widx2 + 1]; + float w3r = wtable[widx3]; + float w3i = isign * wtable[widx3 + 1]; + float w4r = wtable[widx4]; + float w4i = isign * wtable[widx4 + 1]; + + int oidx1 = i + idx2; + int oidx2 = oidx1 + idx0; + int oidx3 = oidx2 + idx0; + int oidx4 = oidx3 + idx0; + int oidx5 = oidx4 + idx0; + out[oidx1] = i3i + tr2 + tr3; + out[oidx1 + 1] = i3r + ti2 + ti3; + out[oidx2] = w1r * dr2 - w1i * di2; + out[oidx2 + 1] = w1r * di2 + w1i * dr2; + out[oidx3] = w2r * dr3 - w2i * di3; + out[oidx3 + 1] = w2r * di3 + w2i * dr3; + out[oidx4] = w3r * dr4 - w3i * di4; + out[oidx4 + 1] = w3r * di4 + w3i * dr4; + out[oidx5] = w4r * dr5 - w4i * di5; + out[oidx5 + 1] = w4r * di5 + w4i * dr5; + } + } + } + } +} 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..f82b45dc 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 @@ -18,6 +18,7 @@ */ package io.github.dsheirer.dsp.filter.channelizer; +import io.github.dsheirer.dsp.fft.MixedRadixComplexFFT; import io.github.dsheirer.dsp.filter.FilterFactory; import io.github.dsheirer.dsp.filter.design.FilterDesignException; import io.github.dsheirer.sample.complex.InterleavedComplexSamples; @@ -27,7 +28,6 @@ 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.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -85,7 +85,8 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann //Sized to process 40 times per second private IFFTProcessorDispatcher mIFFTProcessorDispatcher = new IFFTProcessorDispatcher(25); - private FloatFFT_1D mFFT; + //Allocation-free IFFT. Only used by the IFFT processor dispatcher thread since it reuses a scratch buffer. + private MixedRadixComplexFFT mFFT; private ChannelDFT mChannelDFT; private float[] mInlineSamples; private float[] mInlineFilter; @@ -453,7 +454,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann private void init(float[] coefficients) { int channelCount = getChannelCount(); - mFFT = new FloatFFT_1D(channelCount); + mFFT = new MixedRadixComplexFFT(channelCount); mChannelDFT = ChannelDFTFactory.getChannelDFT(channelCount); int bufferLength = getSubChannelCount() * mTapsPerChannel; mSamplesPerBlock = channelCount; //Same as subChannelCount / 2 @@ -505,7 +506,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.inverse(samples, offset, true); } } diff --git a/src/test/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFTTest.java b/src/test/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFTTest.java new file mode 100644 index 00000000..99317d57 --- /dev/null +++ b/src/test/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFTTest.java @@ -0,0 +1,286 @@ +/* + * ***************************************************************************** + * 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 io.github.dsheirer.dsp.fft; + +import io.github.dsheirer.dsp.filter.channelizer.ComplexPolyphaseChannelizerM2; +import java.lang.management.ManagementFactory; +import java.util.Random; +import org.jtransforms.fft.FloatFFT_1D; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertArrayEquals; +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; +import static org.junit.jupiter.api.Assumptions.assumeTrue; + +/** + * Verifies that the MixedRadixComplexFFT produces the same results as JTransforms FloatFFT_1D and does not allocate + * memory per transform. + */ +public class MixedRadixComplexFFTTest +{ + /** + * Maximum error relative to the largest output magnitude. + */ + private static final double RELATIVE_TOLERANCE = 2e-6; + + /** + * Sizes covering each radix (2, 3, 4, 5), mixed factorizations, powers of two and the channelizer channel counts. + */ + private static final int[] SIZES = {1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 25, 27, 30, 32, 40, 45, 50, 60, 64, 72, + 80, 90, 96, 100, 120, 125, 128, 150, 160, 180, 192, 200, 240, 250, 256, 270, 300, 320, 360, 384, 400, 450, 480, + 500, 512, 600, 640, 720, 750, 768, 800, 900, 960, 1000, 1024, 1080, 1200, 1500, 2000, 2400, 3000, 4096, 6000}; + + private static float[] randomSamples(int length, Random random) + { + float[] samples = new float[length]; + + for(int x = 0; x < samples.length; x++) + { + samples[x] = random.nextFloat() * 2.0f - 1.0f; + } + + return samples; + } + + /** + * Asserts that the actual values match the expected values within a tolerance relative to the largest expected + * magnitude. + */ + private static void assertClose(float[] expected, float[] actual, String message) + { + assertEquals(expected.length, actual.length, message + " length"); + + double maxMagnitude = 1e-30; + double maxError = 0; + + for(int x = 0; x < expected.length; x++) + { + maxMagnitude = Math.max(maxMagnitude, Math.abs(expected[x])); + maxError = Math.max(maxError, Math.abs(expected[x] - actual[x])); + } + + double relativeError = maxError / maxMagnitude; + assertTrue(relativeError <= RELATIVE_TOLERANCE, message + " relative error " + relativeError + + " exceeds " + RELATIVE_TOLERANCE); + } + + /** + * Double precision brute force DFT reference. + * @param samples interleaved complex + * @param sign -1 forward, +1 inverse + */ + private static float[] referenceDFT(float[] samples, int sign) + { + int n = samples.length / 2; + float[] result = new float[samples.length]; + + for(int k = 0; k < n; k++) + { + double re = 0; + double im = 0; + + for(int t = 0; t < n; t++) + { + double angle = sign * 2.0 * Math.PI * ((long)k * t % n) / n; + double cos = Math.cos(angle); + double sin = Math.sin(angle); + re += samples[2 * t] * cos - samples[2 * t + 1] * sin; + im += samples[2 * t] * sin + samples[2 * t + 1] * cos; + } + + result[2 * k] = (float)re; + result[2 * k + 1] = (float)im; + } + + return result; + } + + @Test + public void forwardMatchesJTransforms() + { + Random random = new Random(1); + + for(int size : SIZES) + { + float[] input = randomSamples(2 * size, random); + float[] expected = input.clone(); + new FloatFFT_1D(size).complexForward(expected); + float[] actual = input.clone(); + new MixedRadixComplexFFT(size).forward(actual); + assertClose(expected, actual, "Forward size " + size); + } + } + + @Test + public void inverseMatchesJTransforms() + { + Random random = new Random(2); + + for(int size : SIZES) + { + for(boolean scale : new boolean[]{true, false}) + { + float[] input = randomSamples(2 * size, random); + float[] expected = input.clone(); + new FloatFFT_1D(size).complexInverse(expected, scale); + float[] actual = input.clone(); + new MixedRadixComplexFFT(size).inverse(actual, scale); + assertClose(expected, actual, "Inverse size " + size + " scale " + scale); + } + } + } + + @Test + public void matchesDoublePrecisionReference() + { + Random random = new Random(3); + + for(int size : SIZES) + { + if(size > 1200) + { + continue; //Brute force reference is O(N^2) + } + + MixedRadixComplexFFT fft = new MixedRadixComplexFFT(size); + float[] input = randomSamples(2 * size, random); + + float[] forward = input.clone(); + fft.forward(forward); + assertClose(referenceDFT(input, -1), forward, "Forward reference size " + size); + + float[] inverse = input.clone(); + fft.inverse(inverse, false); + assertClose(referenceDFT(input, 1), inverse, "Inverse reference size " + size); + } + } + + /** + * Mirrors the channelizer usage: repeated in-place inverse transforms of consecutive blocks in a larger array, + * reusing a single instance and its scratch buffer. Also verifies that data outside each block is not modified. + */ + @Test + public void offsetBlocksMatchJTransforms() + { + Random random = new Random(4); + int blocks = 7; + + for(int size : new int[]{2, 72, 96, 100, 128, 240, 400, 480, 600}) + { + int blockLength = 2 * size; + int guard = 3; + float[] input = randomSamples(blocks * blockLength + 2 * guard, random); + + float[] expected = input.clone(); + FloatFFT_1D jtransforms = new FloatFFT_1D(size); + float[] actual = input.clone(); + MixedRadixComplexFFT fft = new MixedRadixComplexFFT(size); + + for(int offset = guard; offset < guard + blocks * blockLength; offset += blockLength) + { + jtransforms.complexInverse(expected, offset, true); + fft.inverse(actual, offset, true); + } + + assertClose(expected, actual, "Offset blocks size " + size); + + for(int x = 0; x < guard; x++) + { + assertEquals(input[x], actual[x], "Leading guard modified size " + size); + assertEquals(input[input.length - 1 - x], actual[actual.length - 1 - x], + "Trailing guard modified size " + size); + } + } + } + + @Test + public void roundTripRestoresInput() + { + Random random = new Random(5); + + for(int size : SIZES) + { + MixedRadixComplexFFT fft = new MixedRadixComplexFFT(size); + float[] input = randomSamples(2 * size, random); + float[] samples = input.clone(); + fft.forward(samples); + fft.inverse(samples, true); + assertArrayEquals(input, samples, 1e-5f, "Round trip size " + size); + } + } + + @Test + public void supportsAllChannelizerChannelCounts() + { + for(double sampleRate = 50_000; sampleRate <= 61_440_000; sampleRate += 25_000) + { + int channelCount = ComplexPolyphaseChannelizerM2.getChannelCount(sampleRate); + assertTrue(MixedRadixComplexFFT.isSupported(channelCount), "Channel count " + channelCount + + " for sample rate " + sampleRate); + } + } + + @Test + public void rejectsUnsupportedSizes() + { + for(int size : new int[]{0, -4, 7, 14, 49, 244, 1001}) + { + assertFalse(MixedRadixComplexFFT.isSupported(size), "Size " + size); + assertThrows(IllegalArgumentException.class, () -> new MixedRadixComplexFFT(size), "Size " + size); + } + } + + /** + * Verifies that repeated transforms do not allocate memory, which is the reason this class exists. + */ + @Test + public void doesNotAllocatePerTransform() + { + assumeTrue(ManagementFactory.getThreadMXBean() instanceof com.sun.management.ThreadMXBean); + com.sun.management.ThreadMXBean threadMXBean = (com.sun.management.ThreadMXBean)ManagementFactory.getThreadMXBean(); + assumeTrue(threadMXBean.isThreadAllocatedMemorySupported() && threadMXBean.isThreadAllocatedMemoryEnabled()); + + int size = 400; + int iterations = 10_000; + MixedRadixComplexFFT fft = new MixedRadixComplexFFT(size); + float[] samples = randomSamples(2 * size, new Random(6)); + + //Warm up so that class loading and JIT compilation are complete before measuring + for(int x = 0; x < iterations; x++) + { + fft.inverse(samples, 0, true); + } + + long threadId = Thread.currentThread().threadId(); + long before = threadMXBean.getThreadAllocatedBytes(threadId); + + for(int x = 0; x < iterations; x++) + { + fft.inverse(samples, 0, true); + } + + long allocated = threadMXBean.getThreadAllocatedBytes(threadId) - before; + + //JTransforms allocates 2N floats (3200 bytes at N=400) per call, about 32 MB for these iterations + assertTrue(allocated < 64 * 1024, "Allocated " + allocated + " bytes over " + iterations + " transforms"); + } +} diff --git a/src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java b/src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java index 92472897..2d4977ec 100644 --- a/src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java +++ b/src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java @@ -19,12 +19,12 @@ package io.github.dsheirer.dsp.filter.channelizer; +import io.github.dsheirer.dsp.fft.MixedRadixComplexFFT; import io.github.dsheirer.sample.complex.InterleavedComplexSamples; import java.util.Arrays; import java.util.List; import java.util.Random; import java.util.concurrent.CopyOnWriteArrayList; -import org.jtransforms.fft.FloatFFT_1D; import org.junit.jupiter.api.Test; import static org.junit.jupiter.api.Assertions.assertEquals; @@ -142,7 +142,9 @@ public class ComplexPolyphaseChannelizerM2Test channelizer.stop(); assertEquals(batchCount, batches.size(), "Batches (channelCount=" + channelCount + ")"); - FloatFFT_1D fft = new FloatFFT_1D(channelCount); + //Use the same IFFT as the channelizer so the comparison remains bit-exact. Equivalence of this IFFT with + //JTransforms is verified separately in MixedRadixComplexFFTTest. + MixedRadixComplexFFT fft = new MixedRadixComplexFFT(channelCount); int blockLength = channelCount * 2; boolean topBlock = TAPS_PER_CHANNEL % 2 == 0; @@ -152,7 +154,7 @@ public class ComplexPolyphaseChannelizerM2Test float[] expected = harness.processFused(topBlock); topBlock = !topBlock; harness.shift(channelCount); - fft.complexInverse(expected, true); + fft.inverse(expected, true); float[] batch = batches.get(block / blocksPerBatch); int offset = (block % blocksPerBatch) * blockLength;