GitHub/sdrtrunk
1
0
Fork 0
mirror of https://github.com/DSheirer/sdrtrunk.git synced 2026-11-09 15:17:24 -05:00
Code Issues Packages Projects Releases 2 Wiki Activity

#2526 Polyphase Channelizer gain imbalance across synthesis of two adjacent channels. (#2527)

Browse Source 
- Fix two-channel synthesis asymmetry across polyphase channel boundaries.
- TwoChannelSynthesizerM2: replace incorrect 2-path polyphase split (which applied
  different responses to each channel) with channel sum + prototype filter
  modulated to the channel boundary (+FS/4).
- ComplexPolyphaseChannelizerM2: start on the middle block map when taps per
  channel is odd so adjacent channel outputs are phase coherent.
- TwoChannelOutputProcessor: translate channel 2 by FS/2 using the absolute
  block index parity derived from the even-sized dispatch batches.
- SynthesizerViewer: translate channel 2 before synthesis.
- Add TwoChannelSynthesizerM2Test for flat, symmetric boundary response.

Co-authored-by: denny <dennis.sheirer@gmail.com>
This commit is contained in:
Denny Sheirer
2026-10-29 08:51:37 -04:00
committed by GitHub
parent 9a6dc2b0a0
commit 248fe3e4f7
6 changed files with 275 additions and 122 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -64,7 +64,11 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
private static final int DEFAULT_MINIMUM_CHANNEL_BANDWIDTH = 25000;
/**
* Determines how many processed channel results to dispatch for threaded IFFT processing per batch
* Determines how many processed channel results to dispatch for threaded IFFT processing per batch.
*
* Note: this value must be even. Each dispatched batch then starts on an even block index, which lets downstream
* channel output processors (e.g. TwoChannelOutputProcessor) derive the absolute block index parity from each
* result's position within the batch.
*/
private static final int PROCESSED_CHANNEL_RESULTS_THRESHOLD = 1024;
@@ -404,6 +408,11 @@ public class ComplexPolyphaseChannelizerM2 extends AbstractComplexPolyphaseChann
int bufferLength = getSubChannelCount() * mTapsPerChannel;
mSamplesPerBlock = channelCount; //Same as subChannelCount / 2
mTopBlockMap = getTopBlockMap(channelCount);
//The prototype filter group delay (taps per channel * channel count / 2 samples) rotates each channel k by
//PI * k * tapsPerChannel, so with an odd taps per channel count adjacent channel outputs are 180 degrees out of
//phase. Starting with the middle block map applies an M/2 circular shift ahead of the IFFT, which multiplies
//channel k by (-1)^k and makes adjacent channels phase coherent for synthesis.
mTopBlockIndicator = (mTapsPerChannel % 2 == 0);
mMiddleBlockMap = getMiddleBlockMap(channelCount);
mInlineFilter = getAlignedFilter(coefficients, channelCount, mTapsPerChannel);
mInlineSamples = new float[bufferLength];

201
src/main/java/io/github/dsheirer/dsp/filter/channelizer/TwoChannelSynthesizerM2.java
View File

@@ -20,46 +20,36 @@ package io.github.dsheirer.dsp.filter.channelizer;
import io.github.dsheirer.sample.complex.ComplexSamples;
import java.util.Arrays;
import org.apache.commons.math3.util.FastMath;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
/**
* Implements a two-channel polyphase filer synthesizer. This class is intended to be used with an M2 polyphase
* channelizer in order to recover a signal that overlaps the boundary between two channels. This synthesizer will
* rejoin the two M2 neighboring channels, mix the signal of interest to baseband, and filter the output stream to
* half of the sample rate, producing an overall M2 output sample rate.
* Implements a two-channel synthesizer. This class is intended to be used with an M2 polyphase channelizer in order
* to recover a signal that overlaps the boundary between two adjacent channels.
*
* The synthesizer concept is modeled on the structure described by Fred Harris et al, in 'Interleaving Different
* Bandwidth Narrowband Channels in Perfect Reconstruction Cascade Polyphase Filter Banks for Efficient Flexible
* Variable Bandwidth Filters in Wideband Digital Transceivers'.
*
* The data buffers and sub-filter kernels bear some resemblance to the structures defined in the paper. The data and
* filter structures in this class are organized to take advantage of Java's ability to leverage SIMD processor
* intrinsics for optimal efficiency in processing. The paper specifies splitting the filter kernel into N polyphase
* partitions and creating N data buffers. This class organizes the data buffer as a contiguous sample
* array and creates an I/Q interleaved filter kernel. Instead of calculating the dot-product for each sub-filter,
* we calculate the product of the full data array against the filter and then accumulate each sub-filter. This allows
* java to use SIMD for the array product and then normal processing for the accumulation. Since this is a two-channel
* processor and the results of each filter accumulation are added, we use a single accumulator across both filters.
* Each M2 channel is sampled at twice the channel bandwidth, so the composite of two adjacent channels can be
* represented at the same sample rate as a single channel: channel 1 remains at baseband and channel 2 is translated
* to FS/2, placing the boundary between the two channels at +FS/4. The caller is responsible for translating channel 2
* by FS/2 (multiplying each sample by (-1)^n, where n is the absolute channelizer output index) since only the caller
* knows the absolute sample index parity. With a perfect reconstruction (-6 dB band edge) channelizer filter, the sum
* of the two channels recovers the original spectrum across the channel boundary.
*
* The composite is then filtered with the supplied low-pass prototype, modulated up by +FS/4 so that it is centered
* on the boundary between the two channels, to suppress the out-of-band channelizer transition regions that alias
* into the outer portions of the composite spectrum. The output remains centered at +FS/4 and the caller is
* responsible for translating it to baseband.
*/
public class TwoChannelSynthesizerM2
{
private final static Logger mLog = LoggerFactory.getLogger(TwoChannelSynthesizerM2.class);
private float[] mSerpentineDataBuffer;
private float[] mIQInterleavedFilter;
private float[] mFilterVectorProduct;
private float mIAccumulator;
private float mQAccumulator;
private FloatFFT_1D mFFT = new FloatFFT_1D(2);
private boolean mTopBlockFlag = true;
private float[] mFilterI;
private float[] mFilterQ;
private float[] mBufferI;
private float[] mBufferQ;
private int mBufferPointer;
/**
* Polyphase synthesizer for mixing two M2 oversampled channels into a composite channel using perfect
* reconstruction filters (-6db at band edge). Output sample rate is the same as the input sample rate of one of
* the channels.
* Synthesizer for combining two M2 oversampled channels into a composite channel. Output sample rate is the same
* as the input sample rate of one of the channels.
*
* @param filter low-pass prototype filter designed at the channel sample rate.
*/
public TwoChannelSynthesizerM2(float[] filter)
{
@@ -74,36 +64,64 @@ public class TwoChannelSynthesizerM2
{
StringBuilder sb = new StringBuilder();
sb.append("Two-Channel Synthesizer");
sb.append("\n\t\t Interleaved Filter: " + Arrays.toString(mIQInterleavedFilter));
sb.append("\n\t\t Serpentine Buffer Length: ").append(mSerpentineDataBuffer.length);
sb.append("\n\t\t Filter Vector Product Length: ").append(mFilterVectorProduct.length);
sb.append("\n\t\t Filter I: " + Arrays.toString(mFilterI));
sb.append("\n\t\t Filter Q: " + Arrays.toString(mFilterQ));
return sb.toString();
}
/**
* Initializes the synthesizer filter and data buffers for operation.
*
* @param filter to use for polyphase synthesis of two channels.
* The prototype low-pass filter is modulated by +FS/4 (multiplied by j^n) to create a complex band-pass filter
* centered on the boundary between the two channels. The filter is stored in time-reversed order so that it
* can be applied directly against the oldest-to-newest ordered sample history.
*
* @param filter low-pass prototype filter.
*/
private void init(float[] filter)
{
int tapsPerChannel = (int) FastMath.ceil(filter.length / 2);
mIQInterleavedFilter = getInterleavedFilter(filter, tapsPerChannel);
mSerpentineDataBuffer = new float[mIQInterleavedFilter.length];
mFilterVectorProduct = new float[mIQInterleavedFilter.length];
int length = filter.length;
mFilterI = new float[length];
mFilterQ = new float[length];
for(int x = 0; x < length; x++)
{
int reversed = length - 1 - x;
switch(x % 4)
{
case 0:
mFilterI[reversed] = filter[x];
break;
case 1:
mFilterQ[reversed] = filter[x];
break;
case 2:
mFilterI[reversed] = -filter[x];
break;
case 3:
mFilterQ[reversed] = -filter[x];
break;
}
}
//Double-length buffers so that the newest filter-length span of samples is always contiguous
mBufferI = new float[length * 2];
mBufferQ = new float[length * 2];
mBufferPointer = 0;
}
/**
* Synthesizes a new channel from the channel 1 and channel 2 sample arrays that are arranged as
* i0, q0, i1, q1 ... iN-1, qN-1
* Synthesizes a new channel from the channel 1 and channel 2 sample buffers. Channel 2 must already be
* translated to FS/2.
*
* @param channelBuffer1 input channel
* @param channelBuffer2 input channel
* @return synthesized channel results of the same length as both channel 1/2 input arrays.
* @param channelBuffer1 input channel at baseband
* @param channelBuffer2 input channel translated to FS/2
* @return synthesized channel results, centered at +FS/4, of the same length as the input buffers.
*/
public ComplexSamples process(ComplexSamples channelBuffer1, ComplexSamples channelBuffer2)
{
if(channelBuffer1.i().length != channelBuffer1.i().length)
if(channelBuffer1.i().length != channelBuffer2.i().length)
{
throw new IllegalArgumentException("Channel 1 and 2 array length must be equal");
}
@@ -116,90 +134,41 @@ public class TwoChannelSynthesizerM2
float[] i = new float[iCh1.length];
float[] q = new float[qCh1.length];
float[] IFFTBuffer = new float[4];
int length = mFilterI.length;
float iAccumulator, qAccumulator;
for(int x = 0; x < iCh1.length; x++)
{
//Load samples from each channel into buffer for IFFT
IFFTBuffer[0] = iCh1[x]; //i
IFFTBuffer[1] = qCh1[x]; //q
IFFTBuffer[2] = iCh2[x]; //i
IFFTBuffer[3] = qCh2[x]; //q
//Combine the channels and load the composite sample into both halves of the double-length buffer
float iSample = iCh1[x] + iCh2[x];
float qSample = qCh1[x] + qCh2[x];
mBufferI[mBufferPointer] = iSample;
mBufferQ[mBufferPointer] = qSample;
mBufferI[mBufferPointer + length] = iSample;
mBufferQ[mBufferPointer + length] = qSample;
mBufferPointer++;
//Perform Inverse FFT (IFFT)
mFFT.complexInverse(IFFTBuffer, true);
//Perform serpentine shift of data blocks in the data buffer - make room for 4 new samples
System.arraycopy(mSerpentineDataBuffer, 0, mSerpentineDataBuffer, 4, mSerpentineDataBuffer.length - 4);
//Top Block - load samples into data buffer in normal order
if(mTopBlockFlag)
if(mBufferPointer >= length)
{
System.arraycopy(IFFTBuffer, 0, mSerpentineDataBuffer, 0, IFFTBuffer.length);
}
//Bottom Block - swap samples via data loading to account for phase shift
else
{
System.arraycopy(IFFTBuffer, 0, mSerpentineDataBuffer, 2, 2);
System.arraycopy(IFFTBuffer, 2, mSerpentineDataBuffer, 0, 2);
mBufferPointer = 0;
}
//Note: in order to use Java's ability to leverage SIMD intrinsics, we perform filtering in two steps
//(multiply then accumulate) since our filter and data are structured with I and Q vectors interleaved.
//This approach allows the Hotspot compiler to more easily recognize the scalor operations.
//Oldest sample is at mBufferPointer and newest sample is at mBufferPointer + length - 1
iAccumulator = 0.0f;
qAccumulator = 0.0f;
//Multiply data samples by the I/Q interleaved filter to form the vector product
for(int y = 0; y < mSerpentineDataBuffer.length; y++)
for(int y = 0; y < length; y++)
{
mFilterVectorProduct[y] = mSerpentineDataBuffer[y] * mIQInterleavedFilter[y];
float iData = mBufferI[mBufferPointer + y];
float qData = mBufferQ[mBufferPointer + y];
iAccumulator += (iData * mFilterI[y]) - (qData * mFilterQ[y]);
qAccumulator += (qData * mFilterI[y]) + (iData * mFilterQ[y]);
}
//Accumulate output I/Q samples from vector product
mIAccumulator = 0.0f;
mQAccumulator = 0.0f;
for(int y = 0; y < mFilterVectorProduct.length; y += 2)
{
mIAccumulator += mFilterVectorProduct[y];
mQAccumulator += mFilterVectorProduct[y + 1];
}
i[x] = mIAccumulator;
q[x] = mQAccumulator;
mTopBlockFlag = !mTopBlockFlag;
i[x] = iAccumulator;
q[x] = qAccumulator;
}
return new ComplexSamples(i, q, channelBuffer1.timestamp());
}
/**
* Creates an interleaved I/Q filter where each coefficient from the filter argument is duplicated and the returned
* filter is twice the length of the original filter.
*
* Note: the returned filter array is sized to: 2 * channel count * taps per channel, which may be slightly more
* than twice the length of the original filter. Any Added filter array elements will contain zero values.
*
* @param coefficients to create an interleaved filter
* @param tapsPerChannel count
* @return
*/
private static float[] getInterleavedFilter(float[] coefficients, int tapsPerChannel)
{
int channelCount = 2;
float[] filter = new float[channelCount * tapsPerChannel * 2];
int coefficientPointer = 0;
int filterPointer = 0;
//Create a new filter that duplicates each tap to produce an interleaved I/Q filter
while(coefficientPointer < coefficients.length)
{
filter[filterPointer++] = coefficients[coefficientPointer];
filter[filterPointer++] = coefficients[coefficientPointer++];
}
return filter;
}
}

2
src/main/java/io/github/dsheirer/dsp/filter/channelizer/output/TwoChannelMixerAssembler.java
View File

@@ -98,7 +98,7 @@ public class TwoChannelMixerAssembler extends MixerAssembler
ComplexSamples channel1 = mChannel1SampleAssembler.getBufferAndReset(timestamp);
ComplexSamples channel2 = mChannel2SampleAssembler.getBufferAndReset(timestamp);
//Synthesize two channels into one, reducing sample rate by 2 so it's same as single channel rate
//Synthesize two channels into one at the single channel sample rate, centered on the channel boundary at +FS/4
ComplexSamples buffer = mTwoChannelSynthesizer.process(channel1, channel2);
//Mix to baseband

18
src/main/java/io/github/dsheirer/dsp/filter/channelizer/output/TwoChannelOutputProcessor.java
View File

@@ -110,10 +110,22 @@ public class TwoChannelOutputProcessor extends ChannelOutputProcessor
@Override
public void process(List<float[]> channelResultsList)
{
for(float[] channelResults : channelResultsList)
//Translate channel 2 up by FS/2 (multiply by (-1)^n) so that it sits adjacent to channel 1 for synthesis. The
//channelizer dispatches even-sized batches, so the index within the batch carries the absolute block parity.
for(int x = 0; x < channelResultsList.size(); x++)
{
mMixerAssembler.receive(channelResults[mChannelOffset1], channelResults[mChannelOffset1 + 1],
channelResults[mChannelOffset2], channelResults[mChannelOffset2 + 1]);
float[] channelResults = channelResultsList.get(x);
if((x & 1) == 0)
{
mMixerAssembler.receive(channelResults[mChannelOffset1], channelResults[mChannelOffset1 + 1],
channelResults[mChannelOffset2], channelResults[mChannelOffset2 + 1]);
}
else
{
mMixerAssembler.receive(channelResults[mChannelOffset1], channelResults[mChannelOffset1 + 1],
-channelResults[mChannelOffset2], -channelResults[mChannelOffset2 + 1]);
}
if(mMixerAssembler.hasBuffer())
{

19
src/main/java/io/github/dsheirer/gui/channelizer/SynthesizerViewer.java
View File

@@ -264,6 +264,7 @@ public class SynthesizerViewer extends JFrame
private TwoChannelSynthesizerM2 mSynthesizer;
private FS4DownConverter mFS4DownConverter = new FS4DownConverter();
private int mSamplesPerCycle = CHANNEL_SAMPLE_RATE / DATA_GENERATOR_FRAME_RATE;
private boolean mOddSample = false;
public DataGenerationManager()
{
@@ -284,7 +285,23 @@ public class SynthesizerViewer extends JFrame
ComplexSamples channel1Buffer = getChannel1ControlPanel().getOscillator().generateComplexSamples(mSamplesPerCycle, 0l);
ComplexSamples channel2Buffer = getChannel2ControlPanel().getOscillator().generateComplexSamples(mSamplesPerCycle, 0l);
ComplexSamples synthesizedBuffer = mSynthesizer.process(channel1Buffer, channel2Buffer);
//Translate channel 2 by FS/2 so that it sits adjacent to channel 1 for synthesis
float[] i2 = channel2Buffer.i().clone();
float[] q2 = channel2Buffer.q().clone();
for(int x = 0; x < i2.length; x++)
{
if(mOddSample)
{
i2[x] = -i2[x];
q2[x] = -q2[x];
}
mOddSample = !mOddSample;
}
ComplexSamples synthesizedBuffer = mSynthesizer.process(channel1Buffer,
new ComplexSamples(i2, q2, channel2Buffer.timestamp()));
getChannel1Panel().receive(new FloatNativeBuffer(channel1Buffer));
getChannel2Panel().receive(new FloatNativeBuffer(channel2Buffer));

146
src/test/java/io/github/dsheirer/dsp/filter/channelizer/TwoChannelSynthesizerM2Test.java Normal file
View File

@@ -0,0 +1,146 @@
/*
* *****************************************************************************
* 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.filter.channelizer;
import io.github.dsheirer.dsp.filter.FilterFactory;
import io.github.dsheirer.sample.complex.ComplexSamples;
import io.github.dsheirer.sample.complex.InterleavedComplexSamples;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
import org.junit.jupiter.api.Test;
import static org.junit.jupiter.api.Assertions.assertEquals;
import static org.junit.jupiter.api.Assertions.assertTrue;
/**
* Verifies that two adjacent M2 channelizer outputs synthesized by {@link TwoChannelSynthesizerM2} recover a tone
* anywhere across the channel boundary with a flat, symmetric response. This requires that adjacent channelizer
* outputs are phase coherent (for both odd and even taps per channel) and that the synthesizer applies the same
* response to each channel.
*/
public class TwoChannelSynthesizerM2Test
{
private static final int CHANNEL_COUNT = 8;
private static final double CHANNEL_BANDWIDTH = 25000.0;
private static final double SAMPLE_RATE = CHANNEL_COUNT * CHANNEL_BANDWIDTH;
private static final int LOWER_CHANNEL = 2;
private static final int SAMPLES = 2048 * CHANNEL_COUNT; //Produces 4096 M2 channel results (4 dispatch batches)
private static final int SETTLE = 100;
/**
* Test: tones at mirrored offsets from the boundary between two channels produce equal output levels, and a tone
* exactly on the boundary produces the same level as one near the center of the composite.
*/
@Test
void boundaryResponseIsFlatAndSymmetric() throws Exception
{
for(int tapsPerChannel : new int[]{8, 9})
{
double center = gain(0.0, tapsPerChannel);
assertEquals(1.0, center, 0.01, "Unity gain at channel boundary (taps per channel=" + tapsPerChannel + ")");
for(double offset : new double[]{0.1, 0.2, 0.3, 0.4})
{
double below = gain(-offset, tapsPerChannel);
double above = gain(offset, tapsPerChannel);
assertEquals(below, above, 0.01, "Asymmetric response at +/-" + offset +
" channel bandwidths from boundary (taps per channel=" + tapsPerChannel + ")");
//Synthesis filter cutoff is 0.55 channel bandwidths, so only check flatness for the inner region
if(offset <= 0.25)
{
assertTrue(below > 0.98, "Passband droop at " + offset + " channel bandwidths: " + below);
}
}
}
}
/**
* Channelizes a tone at the specified offset from the boundary between two channels, synthesizes the two channels
* and returns the output tone amplitude, normalized for the channelizer gain.
* @param offset from the channel boundary as a fraction of the channel bandwidth
* @param tapsPerChannel for the channelizer filter
*/
private static double gain(double offset, int tapsPerChannel) throws Exception
{
List<float[]> results = Collections.synchronizedList(new ArrayList<>());
float[] taps = FilterFactory.getSincM2Channelizer(CHANNEL_BANDWIDTH, CHANNEL_COUNT, tapsPerChannel, false);
ComplexPolyphaseChannelizerM2 channelizer = new ComplexPolyphaseChannelizerM2(taps, (int)SAMPLE_RATE,
CHANNEL_COUNT)
{
@Override
protected void dispatch(List<float[]> channelResultsList)
{
results.addAll(channelResultsList);
}
};
double frequency = (LOWER_CHANNEL + 0.5 + offset) * CHANNEL_BANDWIDTH;
float[] samples = new float[SAMPLES * 2];
for(int x = 0; x < SAMPLES; x++)
{
double phase = 2.0 * Math.PI * frequency * x / SAMPLE_RATE;
samples[2 * x] = (float)Math.cos(phase);
samples[2 * x + 1] = (float)Math.sin(phase);
}
channelizer.start();
channelizer.receive(new InterleavedComplexSamples(samples, 0));
int expected = SAMPLES * 2 / CHANNEL_COUNT;
long timeout = System.currentTimeMillis() + 5000;
while(results.size() < expected && System.currentTimeMillis() < timeout)
{
Thread.sleep(10);
}
channelizer.stop();
assertEquals(expected, results.size(), "Channelizer results");
//Extract both channels and translate the upper channel by FS/2 using the absolute block index parity
float[] i1 = new float[expected], q1 = new float[expected], i2 = new float[expected], q2 = new float[expected];
for(int x = 0; x < expected; x++)
{
float[] result = results.get(x);
float sign = (x & 1) == 0 ? 1.0f : -1.0f;
i1[x] = result[2 * LOWER_CHANNEL];
q1[x] = result[2 * LOWER_CHANNEL + 1];
i2[x] = sign * result[2 * (LOWER_CHANNEL + 1)];
q2[x] = sign * result[2 * (LOWER_CHANNEL + 1) + 1];
}
float[] filter = FilterFactory.getSincM2Synthesizer(2.0 * CHANNEL_BANDWIDTH, CHANNEL_BANDWIDTH, 2, 9);
ComplexSamples output = new TwoChannelSynthesizerM2(filter).process(new ComplexSamples(i1, q1, 0),
new ComplexSamples(i2, q2, 0));
//Correlate against the expected tone, centered at +FS/4 + offset/2 in the composite output
double outputFrequency = 0.25 + offset / 2.0;
double re = 0, im = 0;
int count = 0;
for(int x = SETTLE; x < expected; x++, count++)
{
double angle = -2.0 * Math.PI * outputFrequency * x;
re += output.i()[x] * Math.cos(angle) - output.q()[x] * Math.sin(angle);
im += output.i()[x] * Math.sin(angle) + output.q()[x] * Math.cos(angle);
}
return Math.sqrt(re * re + im * im) / count * CHANNEL_COUNT;
}
}