From 2292af6c07d7c01450b908fe6c4aa2a0d25a17b3 Mon Sep 17 00:00:00 2001 From: John Stockdale Date: Sun, 12 Jul 2026 04:38:20 -0700 Subject: [PATCH] #2428 Polyphase Channelizer efficiency enhancements Closes #2428 Fuse polyphase channelizer multiply-accumulate and reuse scratch accumulator to reduce per-block allocation (#2429) --- .../ComplexPolyphaseChannelizerM2.java | 77 ++--- .../ComplexPolyphaseChannelizerM2Test.java | 279 ++++++++++++++++++ 2 files changed, 323 insertions(+), 33 deletions(-) create mode 100644 src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java 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 3daa92ad..4262a330 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 @@ -24,6 +24,7 @@ import io.github.dsheirer.sample.complex.InterleavedComplexSamples; import io.github.dsheirer.util.Dispatcher; import java.text.DecimalFormat; import java.util.ArrayList; +import java.util.Arrays; import java.util.List; import org.apache.commons.math3.util.FastMath; import org.jtransforms.fft.FloatFFT_1D; @@ -46,10 +47,9 @@ import org.slf4j.LoggerFactory; * * Instead of using an array of channel filters as described in the Harris text, this filter and the sample buffer * are arranged as a contiguous array to maximize Java's ability to leverage native processor Single Instruction - * Multiple Data (SIMD) intrinsics (since Java 8). The filter process is broken into four steps: + * Multiple Data (SIMD) intrinsics (since Java 8). The filter process is broken into three steps: * - * -Multiply the inline array of samples and filter coefficients - * -Accumulate the results for each sub-channel + * -Fused multiply and accumulate of the inline sample/filter arrays into each I/Q sub-channel accumulator * -Rearrange the sub-channel results to correctly order the sub-channels * -Perform IFFT * @@ -73,6 +73,15 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann private FloatFFT_1D mFFT; private float[] mInlineSamples; private float[] mInlineFilter; + + /** + * Reusable sub-channel accumulator. This is scratch space that is fully overwritten on each process() invocation + * and never escapes the method, so it is allocated once here (sized to the sub-channel count) and reused rather + * than being reallocated on every block. process() is only ever invoked from the single-threaded polyphase buffer + * dispatcher (see PolyphaseChannelManager), so reusing this buffer across calls is thread-safe. It is (re)sized in + * init() alongside the other rate-dependent buffers. + */ + private float[] mFilterAccumulator; private boolean mTopBlockIndicator = true; private int[] mTopBlockMap; private int[] mMiddleBlockMap; @@ -333,48 +342,49 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann /** * Processes the sample buffer for each new block of sample data that is loaded and distributes the results to any * registered channel listeners. + * + * Note: mInlineSamples and mInlineFilter are both exactly (subChannelCount * tapsPerChannel) in length (see init() + * and getAlignedFilter()), so the fused multiply-accumulate loop below indexes every element of each array exactly + * once. The multiply and per-sub-channel accumulation are fused into a single pass, which removes the previous + * per-call interim product array entirely and preserves bit-for-bit identical output: for each sub-channel the + * tap products are summed in the same ascending tap order as before, and the JVM does not contract the separate + * float multiply/add into an FMA without an explicit Math.fma() call. */ private float[] process() { - int bufferLength = getSubChannelCount() * mTapsPerChannel; - float[] inlineInterimOutput = new float[bufferLength]; + final int subChannelCount = getSubChannelCount(); + final int tapsPerChannel = mTapsPerChannel; + final float[] inlineSamples = mInlineSamples; + final float[] inlineFilter = mInlineFilter; + final float[] filterAccumulator = mFilterAccumulator; - //Multiply each of the samples by the corresponding filter tap - for(int x = 0; x < mInlineSamples.length; x++) + //Reset the reused accumulator to zero before summing this block's tap products + Arrays.fill(filterAccumulator, 0.0f); + + //Fused multiply-accumulate: multiply each sample by its corresponding filter tap and accumulate the product + //directly into the appropriate I/Q sub-channel accumulator. + for(int tap = 0; tap < tapsPerChannel; tap++) { - inlineInterimOutput[x] = mInlineSamples[x] * mInlineFilter[x]; - } + final int tapOffset = tap * subChannelCount; - float[] filterAccumulator = new float[getSubChannelCount()]; - - int tapOffset = 0; - - //Accumulate the sample/filter product results into each of the I/Q sub-channels - for(int tap = 0; tap < mTapsPerChannel; tap++) - { - tapOffset = tap * getSubChannelCount(); - - for(int channel = 0; channel < getSubChannelCount(); channel++) + for(int channel = 0; channel < subChannelCount; channel++) { - filterAccumulator[channel] += inlineInterimOutput[tapOffset + channel]; + final int index = tapOffset + channel; + filterAccumulator[channel] += inlineSamples[index] * inlineFilter[index]; } } - float[] processed = new float[getSubChannelCount()]; + //This output array is handed off to the IFFT dispatcher thread and therefore must be freshly allocated + //per block (it escapes this method and cannot be reused). + final float[] processed = new float[subChannelCount]; - if(mTopBlockIndicator) + //Rearrange the accumulated sub-channels into the order required by the IFFT, alternating between the + //top-block and middle-block orderings on successive blocks. + final int[] blockMap = mTopBlockIndicator ? mTopBlockMap : mMiddleBlockMap; + + for(int x = 0; x < subChannelCount; x++) { - for(int x = 0; x < getSubChannelCount(); x++) - { - processed[x] = filterAccumulator[mTopBlockMap[x]]; - } - } - else - { - for(int x = 0; x < getSubChannelCount(); x++) - { - processed[x] = filterAccumulator[mMiddleBlockMap[x]]; - } + processed[x] = filterAccumulator[blockMap[x]]; } mTopBlockIndicator = !mTopBlockIndicator; @@ -397,6 +407,7 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann mMiddleBlockMap = getMiddleBlockMap(channelCount); mInlineFilter = getAlignedFilter(coefficients, channelCount, mTapsPerChannel); mInlineSamples = new float[bufferLength]; + mFilterAccumulator = new float[getSubChannelCount()]; } /** 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 new file mode 100644 index 00000000..a92b15de --- /dev/null +++ b/src/test/java/io/github/dsheirer/dsp/filter/channelizer/ComplexPolyphaseChannelizerM2Test.java @@ -0,0 +1,279 @@ +/* + * ***************************************************************************** + * 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.filter.channelizer; + +import java.util.Arrays; +import java.util.Random; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * Verifies that the fused multiply-accumulate used by {@link ComplexPolyphaseChannelizerM2} produces bit-for-bit + * identical output to the prior multiply-then-accumulate implementation. + * + * The channelizer's per-block filter step was changed from (a) multiplying the entire inline sample/filter arrays into + * a per-call interim product array and then summing that array across taps, to (b) fusing the multiply directly into + * the per-sub-channel accumulation in a single pass (which also removes the per-call interim array and reuses the + * accumulator). Both forms are expected to be numerically identical: the tap products for each sub-channel are summed + * in the same ascending tap order, and the JVM does not contract a separate float multiply/add into an FMA without an + * explicit {@code Math.fma()} call. + * + * This test reproduces both the reference and the fused filter+reorder computation (including the same aligned-filter + * arrangement and top/middle block reordering the channelizer uses) and asserts the outputs match exactly, using + * {@link Float#floatToRawIntBits(float)} so that differences in NaN bit patterns and the sign of zero would also be + * detected. Multiple channel counts (spanning the 25 kHz-per-channel sizing from ~50 kHz to 10 MHz sample rates) and + * several stressed floating point input distributions are exercised. + */ +public class ComplexPolyphaseChannelizerM2Test +{ + private static final int TAPS_PER_CHANNEL = 9; //Matches PolyphaseChannelManager.POLYPHASE_CHANNELIZER_TAPS_PER_CHANNEL + + /** + * Test: for every channel count and input distribution, the fused multiply-accumulate output must equal the + * reference multiply-then-accumulate output bit-for-bit. + * + * Success: zero differing output values across all configurations, input distributions, and block phases. + */ + @Test + void fusedMultiplyAccumulateMatchesReference() + { + Random rng = new Random(0xC0FFEEL); + + //Even channel counts corresponding to sampleRate/25000 for ~50 kHz, 100 kHz, 1 MHz, 2.4 MHz, 3.2 MHz, + //6.1 MHz and 10 MHz sample rates, plus small/odd sizes. + int[] channelCounts = {2, 4, 8, 40, 96, 128, 244, 400}; + + for(int channelCount : channelCounts) + { + Harness harness = new Harness(channelCount, TAPS_PER_CHANNEL, rng); + + boolean topBlockReference = true; + boolean topBlockFused = true; + + for(int mode = 0; mode < 4; mode++) + { + for(int iteration = 0; iteration < 1000; iteration++) + { + harness.randomizeSamples(rng, mode); + + float[] reference = harness.processReference(topBlockReference); + topBlockReference = !topBlockReference; + + float[] fused = harness.processFused(topBlockFused); + topBlockFused = !topBlockFused; + + assertBitIdentical(reference, fused, channelCount, mode, iteration); + } + } + } + } + + private static void assertBitIdentical(float[] reference, float[] fused, int channelCount, int mode, int iteration) + { + assertEquals(reference.length, fused.length, + "Output length mismatch (channelCount=" + channelCount + ")"); + + for(int x = 0; x < reference.length; x++) + { + assertEquals(Float.floatToRawIntBits(reference[x]), Float.floatToRawIntBits(fused[x]), + "Fused output differs from reference at index " + x + " (channelCount=" + channelCount + + ", mode=" + mode + ", iteration=" + iteration + "): reference=" + reference[x] + + " fused=" + fused[x]); + } + } + + /** + * Reproduces the channelizer's per-block filter state and both the reference and fused implementations of the + * filter+reorder step so their outputs can be compared directly, without needing the (private) process() method + * or the downstream IFFT. The aligned-filter arrangement and block-reorder maps mirror the production helpers so + * the inputs have realistic structure. + */ + private static final class Harness + { + private final int mSubChannelCount; // channelCount * 2 + private final int mTapsPerChannel; + private final int mBufferLength; // mSubChannelCount * mTapsPerChannel + private final float[] mInlineSamples; // length mBufferLength + private final float[] mInlineFilter; // length mBufferLength + private final int[] mTopBlockMap; // length mSubChannelCount + private final int[] mMiddleBlockMap; // length mSubChannelCount + private final float[] mReusableAccumulator; + + private Harness(int channelCount, int tapsPerChannel, Random rng) + { + mSubChannelCount = channelCount * 2; + mTapsPerChannel = tapsPerChannel; + mBufferLength = mSubChannelCount * tapsPerChannel; + + float[] coefficients = new float[channelCount * tapsPerChannel]; + for(int i = 0; i < coefficients.length; i++) + { + coefficients[i] = rng.nextFloat() * 2f - 1f; + } + + mInlineFilter = getAlignedFilter(coefficients, channelCount, tapsPerChannel); + mInlineSamples = new float[mBufferLength]; + mTopBlockMap = getTopBlockMap(channelCount); + mMiddleBlockMap = getMiddleBlockMap(channelCount); + mReusableAccumulator = new float[mSubChannelCount]; + } + + //Reference: multiply the whole array into an interim buffer, then accumulate across taps (prior behavior). + private float[] processReference(boolean topBlock) + { + float[] interim = new float[mBufferLength]; + for(int x = 0; x < mInlineSamples.length; x++) + { + interim[x] = mInlineSamples[x] * mInlineFilter[x]; + } + + float[] accumulator = new float[mSubChannelCount]; + for(int tap = 0; tap < mTapsPerChannel; tap++) + { + int tapOffset = tap * mSubChannelCount; + for(int channel = 0; channel < mSubChannelCount; channel++) + { + accumulator[channel] += interim[tapOffset + channel]; + } + } + + return reorder(accumulator, topBlock); + } + + //Fused: multiply-accumulate in a single pass into a reused accumulator (new behavior). + private float[] processFused(boolean topBlock) + { + float[] accumulator = mReusableAccumulator; + Arrays.fill(accumulator, 0.0f); + + for(int tap = 0; tap < mTapsPerChannel; tap++) + { + int tapOffset = tap * mSubChannelCount; + for(int channel = 0; channel < mSubChannelCount; channel++) + { + int index = tapOffset + channel; + accumulator[channel] += mInlineSamples[index] * mInlineFilter[index]; + } + } + + return reorder(accumulator, topBlock); + } + + private float[] reorder(float[] accumulator, boolean topBlock) + { + float[] processed = new float[mSubChannelCount]; + int[] map = topBlock ? mTopBlockMap : mMiddleBlockMap; + for(int x = 0; x < mSubChannelCount; x++) + { + processed[x] = accumulator[map[x]]; + } + return processed; + } + + private void randomizeSamples(Random rng, int mode) + { + for(int i = 0; i < mInlineSamples.length; i++) + { + switch(mode) + { + case 0 -> mInlineSamples[i] = rng.nextFloat() * 2f - 1f; //typical baseband + case 1 -> mInlineSamples[i] = (rng.nextInt(5) == 0) ? 0f : rng.nextFloat() * 2f - 1f; //sprinkle zeros + case 2 -> mInlineSamples[i] = (rng.nextFloat() * 2f - 1f) + * (float)Math.pow(10, rng.nextInt(30) - 15); //wide dynamic range + default -> mInlineSamples[i] = switch(rng.nextInt(6)) //edge values + { + case 0 -> -0.0f; + case 1 -> 0.0f; + case 2 -> Float.MIN_VALUE; + case 3 -> -Float.MIN_VALUE; + case 4 -> rng.nextFloat() * 1e6f; + default -> rng.nextFloat() * 2f - 1f; + }; + } + } + } + + //---- Local copies of the production coefficient-arrangement and reorder-map helpers (values only matter in + // that they produce the same realistic structure; the equivalence being tested holds regardless). ---- + + private static float[] getAlignedFilter(float[] coefficients, int channelCount, int tapsPerChannel) + { + float[] filter = new float[channelCount * tapsPerChannel * 2]; + int blockSize = channelCount; + + int coefficientPointer = 0; + int filterPointer = 0; + while(coefficientPointer < coefficients.length) + { + filter[filterPointer++] = coefficients[coefficientPointer]; + filter[filterPointer++] = coefficients[coefficientPointer++]; + } + + for(int x = 0; x < filter.length; x += blockSize) + { + for(int y = 0; y < blockSize / 2; y++) + { + int index1 = x + y; + int index2 = x + (blockSize - y - 1); + float temp = filter[index2]; + filter[index2] = filter[index1]; + filter[index1] = temp; + } + } + + return filter; + } + + private static int[] getTopBlockMap(int channelCount) + { + int[] newMap = new int[channelCount * 2]; + int blockSize = channelCount / 2; + for(int channel = 0; channel < blockSize; channel++) + { + int newIndex = 2 * channel; + int originalIndex = 2 * (blockSize - channel - 1); + int offset = 2 * blockSize; + newMap[originalIndex] = newIndex; + newMap[originalIndex + 1] = newIndex + 1; + newMap[offset + originalIndex] = offset + newIndex; + newMap[offset + originalIndex + 1] = offset + newIndex + 1; + } + return newMap; + } + + private static int[] getMiddleBlockMap(int channelCount) + { + int[] newMap = new int[channelCount * 2]; + int blockSize = channelCount / 2; + for(int channel = 0; channel < blockSize; channel++) + { + int newIndex = 2 * channel; + int originalIndex = 2 * (blockSize - channel - 1); + int offset = 2 * blockSize; + newMap[offset + originalIndex] = newIndex; + newMap[offset + originalIndex + 1] = newIndex + 1; + newMap[originalIndex] = offset + newIndex; + newMap[originalIndex + 1] = offset + newIndex + 1; + } + return newMap; + } + } +}