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

src/main/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFT.java

@@ -1,779 +0,0 @@
-/*
- * *****************************************************************************
- * 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/>
- * ****************************************************************************
- *
- * 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;
-                }
-            }
-        }
-    }
-}

src/main/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2.java

@@ -18,7 +18,6 @@
  */
 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;
@@ -28,6 +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.ReusableFloatInverseFFT;
 import org.slf4j.Logger;
 import org.slf4j.LoggerFactory;
 
@@ -85,8 +85,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
 
     //Sized to process 40 times per second
     private IFFTProcessorDispatcher mIFFTProcessorDispatcher = new IFFTProcessorDispatcher(25);
-    //Allocation-free IFFT.  Only used by the IFFT processor dispatcher thread since it reuses a scratch buffer.
+    private ReusableFloatInverseFFT mFFT;
-    private MixedRadixComplexFFT mFFT;
     private ChannelDFT mChannelDFT;
     private float[] mInlineSamples;
     private float[] mInlineFilter;
@@ -454,7 +453,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
     private void init(float[] coefficients)
     {
         int channelCount = getChannelCount();
-        mFFT = new MixedRadixComplexFFT(channelCount);
+        mFFT = new ReusableFloatInverseFFT(channelCount);
         mChannelDFT = ChannelDFTFactory.getChannelDFT(channelCount);
         int bufferLength = getSubChannelCount() * mTapsPerChannel;
         mSamplesPerBlock = channelCount; //Same as subChannelCount / 2
@@ -506,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.inverse(samples, offset, true);
+                            mFFT.complexInverse(samples, offset);
                         }
                     }
 

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 <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;
+        }
+    }
+}

src/test/java/io/github/dsheirer/dsp/fft/MixedRadixComplexFFTTest.java

@@ -1,286 +0,0 @@
-/*
- * *****************************************************************************
- * 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 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");
-    }
-}

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,9 +142,7 @@ public class ComplexPolyphaseChannelizerM2Test
             channelizer.stop();
             assertEquals(batchCount, batches.size(), "Batches (channelCount=" + channelCount + ")");
 
-            //Use the same IFFT as the channelizer so the comparison remains bit-exact.  Equivalence of this IFFT with
+            FloatFFT_1D fft = new FloatFFT_1D(channelCount);
-            //JTransforms is verified separately in MixedRadixComplexFFTTest.
-            MixedRadixComplexFFT fft = new MixedRadixComplexFFT(channelCount);
             int blockLength = channelCount * 2;
             boolean topBlock = TAPS_PER_CHANNEL % 2 == 0;
 
@@ -154,7 +152,7 @@ public class ComplexPolyphaseChannelizerM2Test
                 float[] expected = harness.processFused(topBlock);
                 topBlock = !topBlock;
                 harness.shift(channelCount);
-                fft.inverse(expected, true);
+                fft.complexInverse(expected, true);
 
                 float[] batch = batches.get(block / blocksPerBatch);
                 int offset = (block % blocksPerBatch) * blockLength;

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 <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");
+        }
+    }
+}