GitHub/sdrtrunk
1
0
Fork 0
mirror of https://github.com/DSheirer/sdrtrunk.git synced 2026-11-12 16:47:31 -05:00
Code Issues Packages Projects Releases 2 Wiki Activity

#2518 P25P1 LSM Decoder enhancements. (#2523)

Browse Source 
* Uses same shared architecture and symbol processor pattern as other C4FM decoders.
* Adds soft symbol processing and sync detection.
* Removes RRC filter and replaces with improved channel filter.
* Viterbi soft symbol decoding with demodulator reliability weighting
* Improved sync reacquisition in low SNR
* Improved handling for phase wrapping and interpolation
* NAC Tracker updated to track by NAC count
* Native buffers returning partial buffers
* Message processor uses only valid LDU1 words
* Reduced phantom TDU messages
* Tuner/Channe frequency error corrections don't reset the symbol trackers and fix possible divide by zero error
* LSM timing error detector now updates per symbol block size instead of sample buffer size - faster convergence
* New unit tests
* Resolves several issues, including a lagging sync detect offset issue in the recently updated P25P1 C4FM decoder.
* Consolidate filter designers for each protocol back into the FilterFactory class.
* DMR, P25 Phase 1 C4FM and LSM, and P25 Phase 2 decoders design the baseband filter once and share the taps across the I and Q filters.

Co-authored-by: denny <dennis.sheirer@gmail.com>
This commit is contained in:
Denny Sheirer
2026-10-29 04:17:47 -04:00
committed by GitHub
parent f2c98cdb40
commit 6b9cd2d317
42 changed files with 3586 additions and 2767 deletions
Download Patch File Download Diff File Expand all files Collapse all files

14
src/main/java/io/github/dsheirer/dsp/c4fm/C4FMMessageFramer.java
View File

@@ -47,6 +47,20 @@ public abstract class C4FMMessageFramer
return process(symbol);
}
/**
* Process the symbol decision along with the soft symbol it was decided from and the soft symbol reliability (ie
* inverse of the soft symbol noise variance). Framers that can weight soft decision decoding override this method.
* Default implementation discards the reliability.
* @param symbol decision
* @param softSymbol equalized soft symbol in radians
* @param reliability of the soft symbol
* @return true if the framer emits a valid message that includes this symbol.
*/
public boolean process(Dibit symbol, float softSymbol, float reliability)
{
return process(symbol, softSymbol);
}
/**
* A sync pattern is detected
*/

109
src/main/java/io/github/dsheirer/dsp/c4fm/C4FMSymbolProcessor.java
View File

@@ -75,10 +75,17 @@ public abstract class C4FMSymbolProcessor
private static final double IDEAL_OUTER_SYMBOL_PHASE = 3.0 * Math.PI / 4.0;
//Allowable gain deviation from the nominal gain, as a fraction of the nominal gain.
private static final double SAMPLE_GAIN_TOLERANCE = 0.1;
//Reliability delay line length, which must exceed the symbol equalizer delay.
private static final int RELIABILITY_DELAY_LINE_LENGTH = 32;
private double mNoiseStandardDeviationThreshold;
private double mSamplePoint;
protected double mSamplesPerSymbol;
protected float[] mBuffer;
//Optional soft symbol reliabilities for each sample, parallel to the sample buffer (see process(float[], float[]))
private float[] mReliabilityBuffer;
//Delays the sampled reliabilities to align with the soft symbols emerging from the symbol equalizer delay line
private final float[] mReliabilityDelayLine = new float[RELIABILITY_DELAY_LINE_LENGTH];
private int mReliabilityDelayLinePointer = 0;
private float mLaggingSyncOffset;
protected float mOptimizeFineIncrement;
private int mBufferPointer;
@@ -216,6 +223,16 @@ public abstract class C4FMSymbolProcessor
*/
public abstract float getSyncThresholdEqualized();
/**
* Minimum correlation against the reference sync waveform, as a fraction of the ideal correlation score, for the
* first two sync detections following a sync loss (see validateSyncWaveform()). Default is 85%.
* @return validation threshold as a fraction of the ideal correlation score.
*/
public double getSyncWaveformValidationThreshold()
{
return 0.85;
}
/**
* Notification that a sync detection was validated and the correction applied. Subclasses can override this to
* capture protocol specific details (ie the detected sync pattern) from the sync detector that triggered the
@@ -269,6 +286,19 @@ public abstract class C4FMSymbolProcessor
* @param samples to process
*/
public void process(float[] samples)
{
process(samples, null);
}
/**
* Primary input method for receiving a stream of demodulated samples, with an optional parallel stream of soft
* symbol reliabilities (ie inverse of the soft symbol noise variance) for each sample, to process into symbols.
* When provided, each soft symbol's reliability is sampled at the symbol sampling instant and is provided to the
* message framer with the soft symbol to weight soft decision decoding.
* @param samples to process
* @param reliabilities for each sample, or null.
*/
public void process(float[] samples, float[] reliabilities)
{
//Establish buffer size based on incoming sample array length which should be consistent across decoder lifecycle
if(mBuffer == null)
@@ -278,6 +308,11 @@ public abstract class C4FMSymbolProcessor
mBufferPointer = mBufferReloadThreshold;
}
if(reliabilities != null && mReliabilityBuffer == null)
{
mReliabilityBuffer = new float[mBuffer.length];
}
//Create shadow copy of heap variables on stack
double samplePoint = mSamplePoint;
double samplesPerSymbol = mSamplesPerSymbol;
@@ -319,6 +354,22 @@ public abstract class C4FMSymbolProcessor
int copyLength = samples.length;
System.arraycopy(mBuffer, copyLength, mBuffer, 0, mBuffer.length - copyLength);
System.arraycopy(samples, samplesPointer, mBuffer, mBuffer.length - copyLength, copyLength);
if(mReliabilityBuffer != null)
{
System.arraycopy(mReliabilityBuffer, copyLength, mReliabilityBuffer, 0, mBuffer.length - copyLength);
if(reliabilities != null && reliabilities.length == samples.length)
{
System.arraycopy(reliabilities, samplesPointer, mReliabilityBuffer, mBuffer.length - copyLength,
copyLength);
}
else
{
Arrays.fill(mReliabilityBuffer, mBuffer.length - copyLength, mBuffer.length, Float.NaN);
}
}
samplesPointer += copyLength;
bufferPointer -= copyLength;
}
@@ -352,6 +403,13 @@ public abstract class C4FMSymbolProcessor
softSymbol = mSampleEqualizer.getEqualizedSymbol(mBuffer[bufferPointer], mBuffer[bufferPointer + 1], samplePoint);
float timingSymbol = softSymbol;
if(mReliabilityBuffer != null)
{
mReliabilityDelayLinePointer = (mReliabilityDelayLinePointer + 1) % RELIABILITY_DELAY_LINE_LENGTH;
mReliabilityDelayLine[mReliabilityDelayLinePointer] = (float)(mReliabilityBuffer[bufferPointer] +
(mReliabilityBuffer[bufferPointer + 1] - mReliabilityBuffer[bufferPointer]) * samplePoint);
}
//Process dynamic updates to the PLL when supported by the decoder
softSymbol = mSampleEqualizer.correct(softSymbol);
@@ -413,7 +471,21 @@ public abstract class C4FMSymbolProcessor
// 2 consecutive sync detections, send the current sample equalizer's balance value as a pseudo-PLL
// error measurement so the tuner can self adjust to correct PPM. The balance is applied in the
// equalized symbol domain, so divide by gain to convert to radians per sample for the decoder.
if(mMessageFramer.process(symbol, softSymbol) && mConsecutiveSyncPeriods > 1)
boolean validMessage;
if(mReliabilityBuffer != null)
{
int delay = Math.min(mSymbolEqualizer.getDelaySymbolCount(), RELIABILITY_DELAY_LINE_LENGTH - 1);
float reliability = mReliabilityDelayLine[(mReliabilityDelayLinePointer - delay +
RELIABILITY_DELAY_LINE_LENGTH) % RELIABILITY_DELAY_LINE_LENGTH];
validMessage = mMessageFramer.process(symbol, softSymbol, reliability);
}
else
{
validMessage = mMessageFramer.process(symbol, softSymbol);
}
if(validMessage && mConsecutiveSyncPeriods > 1)
{
mFeedbackDecoder.processPLLError(mSampleEqualizer.balanceOptimized() / mSampleEqualizer.gain());
}
@@ -542,14 +614,14 @@ public abstract class C4FMSymbolProcessor
for(int x = referenceSyncWaveform.samples().length - 1; x >= 0; x--)
{
index = integral + x;
float sample = LinearInterpolator.calculate(mBuffer[index], mBuffer[index + 1], mu);
float sample = interpolate(mBuffer[index], mBuffer[index + 1], mu);
sample *= gain;
sample += balance;
sample = Math.clamp(sample, -referenceSyncWaveform.maxSample(), referenceSyncWaveform.maxSample());
accumulator += referenceSyncWaveform.samples()[x] * sample;
}
double threshold = referenceSyncWaveform.idealCorrelationScore() * 0.85; //85% of the ideal correlation score
double threshold = referenceSyncWaveform.idealCorrelationScore() * getSyncWaveformValidationThreshold();
//Return the correction as valid if the correlation score exceeds the threshold, otherwise return the static
// invalid sync detection to nullify the correction.
@@ -1145,7 +1217,7 @@ public abstract class C4FMSymbolProcessor
}
/**
* Equalizes both samples with optional phase wrap detection and returns the interpolated value at the mu offset.
* Interpolates between the samples at the mu offset and equalizes the interpolated value.
*
* @param sample1 of the sample stream
* @param sample2 of the sample stream
@@ -1154,10 +1226,7 @@ public abstract class C4FMSymbolProcessor
*/
public float getEqualizedSymbol(float sample1, float sample2, double mu)
{
sample1 = equalize(sample1);
sample2 = equalize(sample2);
//Calculate interpolated sample and pass through the dynamic balance correction method.
return LinearInterpolator.calculate(sample1, sample2, mu);
return equalize(interpolate(sample1, sample2, mu));
}
/**
@@ -1179,6 +1248,21 @@ public abstract class C4FMSymbolProcessor
}
/**
* Interpolates the (unequalized) sample value between two consecutive samples of the demodulated sample stream.
* Default implementation is linear interpolation. Protocols with a phase wrapped sample stream (ie P25 LSM) can
* override this to interpolate across the phase wrap.
*
* @param sample1 of the sample stream
* @param sample2 of the sample stream
* @param mu interpolation point between sample 1 and sample 2, in the range 0.0 to 1.0
* @return interpolated sample
*/
protected float interpolate(float sample1, float sample2, double mu)
{
return LinearInterpolator.calculate(sample1, sample2, mu);
}
/**
* Calculates the sync correlation score at the specified offset and samples per symbol interval using the
* current equalizer gain and supplied balance value.
@@ -1209,16 +1293,13 @@ public abstract class C4FMSymbolProcessor
double fractional = pointer - bufferPointer;
float score = 0;
float sample1, sample2;
float idealSymbol;
for(int x = 0; x < patternSymbols.length; x++)
{
if(bufferPointer < maxPointer)
{
sample1 = mBuffer[bufferPointer] * gain + balance;
sample2 = mBuffer[bufferPointer + 1] * gain + balance;
softSymbol = LinearInterpolator.calculate(sample1, sample2, fractional);
softSymbol = interpolate(mBuffer[bufferPointer], mBuffer[bufferPointer + 1], fractional) * gain + balance;
}
else
{
@@ -1309,8 +1390,8 @@ public abstract class C4FMSymbolProcessor
if(resampleStartIntegral >= 0)
{
symbol = patternSymbols[x];
resampledSoftSymbol = LinearInterpolator.calculate(mBuffer[resampleStartIntegral],
mBuffer[resampleStartIntegral + 1], resampleStart - resampleStartIntegral);
resampledSoftSymbol = interpolate(mBuffer[resampleStartIntegral], mBuffer[resampleStartIntegral + 1],
resampleStart - resampleStartIntegral);
resampledSoftSymbol = resampledSoftSymbol * mSampleEqualizer.mGain + measuredBalance;
resampledDibit = Dibit.fromSample(resampledSoftSymbol);
bitErrorCount += patternDibits[x].getBitErrorFrom(resampledDibit);

675
src/main/java/io/github/dsheirer/dsp/filter/FilterFactory.java
View File

@@ -57,8 +57,11 @@ import io.github.dsheirer.vector.calibrate.CalibrationManager;
import io.github.dsheirer.vector.calibrate.CalibrationType;
import io.github.dsheirer.vector.calibrate.Implementation;
import java.text.DecimalFormat;
import java.util.Arrays;
import java.util.Map;
import java.util.Set;
import java.util.TreeSet;
import java.util.concurrent.ConcurrentHashMap;
import org.apache.commons.math3.util.FastMath;
import org.jtransforms.fft.FloatFFT_1D;
import org.slf4j.Logger;
@@ -72,6 +75,7 @@ public class FilterFactory
private final static Logger mLog = LoggerFactory.getLogger(FilterFactory.class);
private static final double PERFECT_RECONSTRUCTION_GAIN_AT_BAND_EDGE = -6.020599842071533; //decibel(0.5, 0.0)
private static final double MARGIN_OF_ERROR = 0.0003;
private static final Map<LowPassSpecification,float[]> REMEZ_LOW_PASS_FILTERS = new ConcurrentHashMap<>();
/**
* Generates coefficients for a unity-gain, windowed low-pass filter
@@ -664,111 +668,6 @@ public class FilterFactory
return taps;
}
/**
* Creates a Root Raised Cosine filter.
* @param samplesPerSymbol
* @param symbolCount
* @param alpha
* @return
*/
public static float[] getRRC(double samplesPerSymbol, int symbolCount, float alpha)
{
//TODO: this wraps the existing method which is reducing the samples per symbol value x2 internally and is likely
// a mistake. This method is used by NXDN to correct that error until the DMR and P25 decoders can be updated to
// ensure that the correction doesn't impact those decoders. Once the cleanup is complete, this method can be
// deleted and the NXDN implementation can revert to using the future corrected version.
return getRootRaisedCosine(samplesPerSymbol * 2, symbolCount, alpha);
}
/**
* Creates root raised cosine filter coefficients with a tap count equal to the symbols x samplesPerSymbol + 1.
*
* Symbol count should be an even integer.
* Ported to java from gnuradio/filter/firdes.cc
*
* For 40db attenuation, calculate the number of symbols based on the following formula:
* Symbols = -44 * alpha + 33
*
* Polyphase Channelizer notes:
* -Set samples Per Symbol at 2 (or more) * number of channels
* -Set symbolCount to sufficient size to produce required attenuation
* -Set symbol rate in hertz
* -Set alpha = (desired channel bandwidth / (symbol rate * samples per symbol) - 1.0
* -Channel bandwidth must be larger than symbol rate * samples per symbol, in order for there to be a positive alpha
*
* @param samplesPerSymbol - number of samples per symbol
* @param symbolCount - number of symbol intervals for the filter. This directly impacts the filter size
* @param alpha - roll-off factor.
* @param gain of the filter
* @return - filter coefficients
*/
public static float[] getRootRaisedCosine(double samplesPerSymbol, int symbolCount, float alpha)
{
samplesPerSymbol /= 2;
int taps = (int)Math.round(samplesPerSymbol * symbolCount);
float scale = 0;
float[] coefficients = new float[taps];
for(int x = 0; x < taps; x++)
{
double index = x - (taps / 2.0f);
double x1 = FastMath.PI * index / samplesPerSymbol;
double x2 = 4.0 * alpha * index / samplesPerSymbol;
double x3 = x2 * x2 - 1.0;
double numerator, denominator;
if(FastMath.abs(x3) >= 0.000001)
{
if(x != taps / 2)
{
numerator = FastMath.cos((1.0 + alpha) * x1) +
FastMath.sin((1.0 - alpha) * x1) / (4.0 * alpha * index / samplesPerSymbol);
}
else
{
numerator = FastMath.cos((1.0 + alpha) * x1) + (1.0 - alpha) * FastMath.PI / (4.0 * alpha);
}
denominator = x3 * FastMath.PI;
}
else
{
if(alpha == 1.0)
{
coefficients[x] = -1.0f;
continue;
}
x3 = (1.0 - alpha) * x1;
x2 = (1.0 + alpha) * x1;
numerator = (FastMath.sin(x2) * (1.0 + alpha) * FastMath.PI - FastMath.cos(x3) * ((1.0 - alpha) * FastMath.PI *
samplesPerSymbol) / (4.0 * alpha * index) + FastMath.sin(x3) * samplesPerSymbol *
samplesPerSymbol / (4.0 * alpha * index * index));
denominator = -32.0 * FastMath.PI * alpha * alpha * index / samplesPerSymbol;
}
coefficients[x] = (float)(4.0 * alpha * numerator / denominator);
scale += coefficients[x];
}
/**
* Normalize (scale) coefficients to 1.0 sum and apply gain
*/
for(int x = 0; x < taps; x++)
{
coefficients[x] = coefficients[x] / scale;
}
return coefficients;
}
/**
* Creates a filter from the filter specification using the remez exchange design algorithm
*
@@ -788,6 +687,572 @@ public class FilterFactory
return null;
}
/**
* Creates a low pass filter using the remez exchange design algorithm with 0.01 dB pass band and stop band ripple.
* See getRemezLowPass(double, double, double, double, double).
*
* @param sampleRate in Hertz
* @param passBandCutoff in Hertz
* @param stopBandStart in Hertz
* @return filter coefficients
* @throws IllegalStateException if the filter cannot be designed
*/
public static float[] getRemezLowPass(double sampleRate, double passBandCutoff, double stopBandStart)
{
return getRemezLowPass(sampleRate, passBandCutoff, 0.01, stopBandStart, 0.01);
}
/**
* Creates a unity gain low pass filter using the remez exchange design algorithm. Designed filters are cached and
* reused, since the design is computationally expensive and decoders request the same filters each time they start.
*
* @param sampleRate in Hertz
* @param passBandCutoff in Hertz
* @param passBandRipple in decibels
* @param stopBandStart in Hertz
* @param stopBandRipple in decibels
* @return filter coefficients, a copy that the caller can modify (ie RealFIRFilter reverses its coefficients)
* @throws IllegalStateException if the filter cannot be designed
*/
public static float[] getRemezLowPass(double sampleRate, double passBandCutoff, double passBandRipple,
double stopBandStart, double stopBandRipple)
{
LowPassSpecification key = new LowPassSpecification(sampleRate, passBandCutoff, passBandRipple, stopBandStart,
stopBandRipple);
//Filter cache is shared across decoder instances that may be started from different threads
float[] coefficients = REMEZ_LOW_PASS_FILTERS.get(key);
if(coefficients == null)
{
FIRFilterSpecification specification = FIRFilterSpecification
.lowPassBuilder()
.sampleRate(sampleRate)
.passBandCutoff(passBandCutoff)
.passBandAmplitude(1.0).passBandRipple(passBandRipple)
.stopBandAmplitude(0.0).stopBandStart(stopBandStart)
.stopBandRipple(stopBandRipple).build();
try
{
coefficients = getTaps(specification);
if(coefficients != null)
{
REMEZ_LOW_PASS_FILTERS.put(key, coefficients);
}
}
catch(Exception e) //FilterDesignException or runtime failures from the Remez designer
{
mLog.error("Error designing low pass filter for " + key, e);
}
}
if(coefficients == null)
{
throw new IllegalStateException("Unable to design low pass filter for " + key);
}
return Arrays.copyOf(coefficients, coefficients.length);
}
/**
* Remez low pass filter design parameters, used as the filter cache key.
*/
private record LowPassSpecification(double sampleRate, double passBandCutoff, double passBandRipple,
double stopBandStart, double stopBandRipple)
{
@Override
public String toString()
{
return "sample rate [" + sampleRate + "] pass band [" + passBandCutoff + "] ripple [" + passBandRipple +
"] stop band [" + stopBandStart + "] ripple [" + stopBandRipple + "]";
}
}
/**
* Designs an inverse sinc (x / sin(x)) filter that compensates for the sinc (integrate and dump) response of a
* transmitter's shaping filter (ie NXDN). The response is inverse sinc up to the band edge of the companion root
* raised cosine filter at (1 + alpha) * symbolRate / 2, and is held at the band edge gain above it.
*
* @param sampleRate in Hertz
* @param symbolRate in symbols per second
* @param alpha roll-off factor of the companion RRC filter, defining the band edge (1 + alpha) * symbolRate / 2
* @param symbolSpan filter length in symbol periods
* @return filter taps, odd length (minimum 9), linear phase, unity gain at DC
*/
public static float[] getInverseSincFilter(double sampleRate, double symbolRate, double alpha, int symbolSpan)
{
if(sampleRate <= 0 || symbolRate <= 0)
{
throw new IllegalArgumentException("Sample rate and symbol rate must be positive");
}
if(alpha < 0 || alpha > 1)
{
throw new IllegalArgumentException("Alpha must be in the range 0 - 1");
}
if(symbolSpan < 1)
{
throw new IllegalArgumentException("Symbol span must be at least 1");
}
double bandEdge = (1.0 + alpha) * symbolRate / 2.0;
if(bandEdge >= sampleRate / 2.0)
{
throw new IllegalArgumentException("Sample rate [" + sampleRate +
"] is too low for symbol rate [" + symbolRate + "] - band edge exceeds Nyquist");
}
int length = Math.max(getOddLength(symbolSpan * sampleRate / symbolRate), 9);
double symbolPeriod = 1.0 / symbolRate;
double edgeGain = inverseSinc(bandEdge * symbolPeriod);
//Least-squares weight applied beyond the band edge, relative to a weight of 1.0 in band.
double stopBandWeight = 0.01;
double[] taps = getLeastSquaresSymmetricTaps(sampleRate, length, frequency -> frequency <= bandEdge ?
new double[]{inverseSinc(frequency * symbolPeriod), 1.0} : new double[]{edgeGain, stopBandWeight});
return normalizeToUnityDC(taps);
}
/**
* Designs a root raised cosine filter from the closed form impulse response, normalized to unity gain at DC.
*
* @param sampleRate in Hertz
* @param symbolRate in symbols per second
* @param alpha roll-off factor (0 - 1]
* @param symbolSpan filter length in symbol periods
* @return filter taps, odd length, linear phase
*/
public static float[] getRootRaisedCosineFilter(double sampleRate, double symbolRate, double alpha, int symbolSpan)
{
double[] taps = getRootRaisedCosineTaps(sampleRate, symbolRate, alpha, symbolSpan);
float[] output = new float[taps.length];
for(int n = 0; n < taps.length; n++)
{
output[n] = (float)taps[n];
}
return output;
}
/**
* Designs a receive filter that is a root raised cosine filter convolved with an inverse sinc filter (see
* getInverseSincFilter()). Cascaded with a transmit filter that is a root raised cosine and sinc filter (ie NXDN),
* the overall response is a raised cosine Nyquist filter with no intersymbol interference at the symbol sampling
* instants.
*
* @param sampleRate in Hertz
* @param symbolRate in symbols per second
* @param alpha roll-off factor (0 - 1]
* @param rrcSymbolSpan root raised cosine filter length in symbol periods
* @param inverseSincSymbolSpan inverse sinc filter length in symbol periods
* @return filter taps, odd length (sum of both filter lengths - 1), linear phase, unity gain at DC
*/
public static float[] getRootRaisedCosineInverseSincFilter(double sampleRate, double symbolRate, double alpha,
int rrcSymbolSpan, int inverseSincSymbolSpan)
{
double[] rrc = getRootRaisedCosineTaps(sampleRate, symbolRate, alpha, rrcSymbolSpan);
float[] inverseSinc = getInverseSincFilter(sampleRate, symbolRate, alpha, inverseSincSymbolSpan);
//Full linear convolution - both filters are odd length and symmetric, so the result is too
double[] combined = new double[rrc.length + inverseSinc.length - 1];
for(int i = 0; i < rrc.length; i++)
{
for(int j = 0; j < inverseSinc.length; j++)
{
combined[i + j] += rrc[i] * inverseSinc[j];
}
}
//Normalize to unity gain at DC so the discriminator output keeps its deviation scaling
return normalizeToUnityDC(combined);
}
/**
* Designs a band limited integrate and dump receive filter for the FM discriminator output (ie P25 Phase 1 C4FM).
*
* TIA-102.BAAA specifies a raised cosine (alpha 0.2) Nyquist filter cascaded with an inverse sinc (x / sin(x))
* shaping filter at the transmitter, and an integrate and dump filter (sinc response) at the receiver, so that the
* overall response is the raised cosine with no intersymbol interference at the symbol sampling instants.
*
* The transmitted signal is band limited to (1 + alpha) * symbolRate / 2, but an integrate and dump filter has
* slowly decaying sidelobes that pass the FM discriminator noise above the signal band, where the discriminator
* noise power spectral density rises with the square of frequency. This produces a linear phase FIR filter with
* the integrate and dump (sinc) response across the signal band and a stopband above it, which preserves the zero
* ISI property while rejecting the out of band noise.
*
* @param sampleRate in Hertz
* @param symbolRate in symbols per second
* @param alpha roll-off factor of the transmitter's raised cosine filter
* @param stopBandStart frequency in Hertz where the stopband starts. The passband extends to the signal band
* edge of (1 + alpha) * symbolRate / 2 and the region between the passband edge and the stopband start is a
* transition band that is excluded from the fit.
* @param symbolSpan filter length in symbol periods
* @return filter taps, odd length, linear phase, unity gain at DC
*/
public static float[] getBandLimitedIntegrateAndDump(double sampleRate, double symbolRate, double alpha,
double stopBandStart, int symbolSpan)
{
if(sampleRate <= 0 || symbolRate <= 0)
{
throw new IllegalArgumentException("Sample rate and symbol rate must be positive");
}
if(symbolSpan < 1)
{
throw new IllegalArgumentException("Symbol span must be at least 1");
}
double nyquist = sampleRate / 2.0;
double passBandEdge = (1.0 + alpha) * symbolRate / 2.0;
if(stopBandStart <= passBandEdge || stopBandStart >= nyquist)
{
throw new IllegalArgumentException("Stop band start [" + stopBandStart + "] must be between the pass band " +
"edge [" + passBandEdge + "] and Nyquist [" + nyquist + "]");
}
int length = getOddLength(symbolSpan * sampleRate / symbolRate);
double symbolPeriod = 1.0 / symbolRate;
//Least-squares weight applied in the stopband, relative to a weight of 1.0 in the passband.
double stopBandWeight = 1.0;
double[] taps = getLeastSquaresSymmetricTaps(sampleRate, length, frequency -> {
if(frequency <= passBandEdge)
{
return new double[]{sinc(frequency * symbolPeriod), 1.0};
}
else if(frequency >= stopBandStart)
{
return new double[]{0.0, stopBandWeight};
}
return null; //Transition band
});
return normalizeToUnityDC(taps);
}
/**
* Designs a one symbol period integrator (moving sum) for the FM discriminator output of a differentially encoded
* QPSK (ie P25 LSM / CQPSK) signal. The FM discriminator produces the phase change between consecutive samples,
* so summing one symbol period of the discriminator output produces the (unwrapped) phase change across the symbol
* period, which is the differential demodulation of the symbol: arg(x[n] * conj(x[n - samplesPerSymbol])), modulo
* 2 * PI.
*
* A fractional samples per symbol value is handled by weighting the final tap with the fractional remainder, which
* is equivalent to linearly interpolating the delayed (unwrapped) phase.
*
* @param samplesPerSymbol must be at least 1.0
* @return filter taps scaled to unity gain at DC so the output is in radians per sample, like the C4FM receive
* filter output.
*/
public static float[] getDifferentialIntegrator(double samplesPerSymbol)
{
if(samplesPerSymbol < 1.0)
{
throw new IllegalArgumentException("Samples per symbol [" + samplesPerSymbol + "] must be at least 1.0");
}
int integral = (int)Math.floor(samplesPerSymbol);
double fractional = samplesPerSymbol - integral;
float[] taps = new float[fractional > 0.0 ? integral + 1 : integral];
Arrays.fill(taps, 0, integral, (float)(1.0 / samplesPerSymbol));
if(fractional > 0.0)
{
taps[integral] = (float)(fractional / samplesPerSymbol);
}
return taps;
}
/**
* Root raised cosine filter taps in double precision, normalized to unity gain at DC. See
* getRootRaisedCosineFilter().
*/
private static double[] getRootRaisedCosineTaps(double sampleRate, double symbolRate, double alpha, int symbolSpan)
{
if(sampleRate <= 0 || symbolRate <= 0)
{
throw new IllegalArgumentException("Sample rate and symbol rate must be positive");
}
if(alpha <= 0 || alpha > 1)
{
throw new IllegalArgumentException("Alpha must be in the range (0 - 1]");
}
if(symbolSpan < 1)
{
throw new IllegalArgumentException("Symbol span must be at least 1");
}
int length = getOddLength(symbolSpan * sampleRate / symbolRate);
int middle = (length - 1) / 2;
double[] taps = new double[length];
double sum = 0.0;
for(int n = 0; n < length; n++)
{
//Time in symbol periods relative to the filter center
double t = (n - middle) * symbolRate / sampleRate;
double value;
if(Math.abs(t) < 1e-9)
{
value = 1.0 - alpha + 4.0 * alpha / Math.PI;
}
else if(Math.abs(Math.abs(t) - 1.0 / (4.0 * alpha)) < 1e-9)
{
//Limit at the singularity t = +/- T / (4 * alpha)
value = (alpha / Math.sqrt(2.0)) *
((1.0 + 2.0 / Math.PI) * Math.sin(Math.PI / (4.0 * alpha)) +
(1.0 - 2.0 / Math.PI) * Math.cos(Math.PI / (4.0 * alpha)));
}
else
{
double numerator = Math.sin(Math.PI * t * (1.0 - alpha)) +
4.0 * alpha * t * Math.cos(Math.PI * t * (1.0 + alpha));
double denominator = Math.PI * t * (1.0 - Math.pow(4.0 * alpha * t, 2.0));
value = numerator / denominator;
}
taps[n] = value;
sum += value;
}
for(int n = 0; n < length; n++)
{
taps[n] /= sum;
}
return taps;
}
/**
* Odd tap count that is at least the specified length, so the filter has an integer group delay of
* (length - 1) / 2 samples.
*/
private static int getOddLength(double length)
{
int oddLength = (int)Math.ceil(length);
if(oddLength % 2 == 0)
{
oddLength++;
}
return oddLength;
}
/**
* Desired response and least-squares weight for each frequency grid point of getLeastSquaresSymmetricTaps().
*/
private interface LeastSquaresObjective
{
/**
* @param frequency in Hertz
* @return {desired response, weight} at the frequency, or null to exclude it from the fit (ie transition band)
*/
double[] evaluate(double frequency);
}
/**
* Designs a symmetric (linear phase) FIR filter whose response is a weighted least-squares fit to the objective.
*
* @param sampleRate in Hertz
* @param length of the filter, odd
* @param objective desired response and weight at each frequency
* @return filter taps, not normalized
*/
private static double[] getLeastSquaresSymmetricTaps(double sampleRate, int length, LeastSquaresObjective objective)
{
double nyquist = sampleRate / 2.0;
int middle = (length - 1) / 2;
int coefficientCount = middle + 1;
//A symmetric filter has the zero-phase response A(f) = c[0] + sum(m=1..middle) c[m] * 2cos(2*pi*f*m/fs)
//where c[m] = h[middle - m] = h[middle + m]. Accumulate the weighted normal equations (B'WB)c = B'Wd over a
//dense frequency grid (16 points per coefficient) from DC to Nyquist.
double[][] normal = new double[coefficientCount][coefficientCount];
double[] rhs = new double[coefficientCount];
double[] basis = new double[coefficientCount];
int gridSize = 16 * coefficientCount;
for(int g = 0; g <= gridSize; g++)
{
double frequency = nyquist * g / gridSize;
double[] desiredAndWeight = objective.evaluate(frequency);
if(desiredAndWeight == null)
{
continue;
}
double desired = desiredAndWeight[0];
double weight = desiredAndWeight[1];
basis[0] = 1.0;
for(int m = 1; m < coefficientCount; m++)
{
basis[m] = 2.0 * Math.cos(2.0 * Math.PI * frequency * m / sampleRate);
}
for(int i = 0; i < coefficientCount; i++)
{
double weighted = weight * basis[i];
rhs[i] += weighted * desired;
for(int j = i; j < coefficientCount; j++)
{
normal[i][j] += weighted * basis[j];
}
}
}
//Fill in the lower triangle of the symmetric normal matrix
for(int i = 0; i < coefficientCount; i++)
{
for(int j = 0; j < i; j++)
{
normal[i][j] = normal[j][i];
}
}
double[] coefficients = solve(normal, rhs);
//Expand to the full symmetric impulse response
double[] taps = new double[length];
for(int m = 0; m < coefficientCount; m++)
{
taps[middle - m] = coefficients[m];
taps[middle + m] = coefficients[m];
}
return taps;
}
/**
* Normalizes the filter taps to unity gain at DC.
* @param taps to normalize
* @return normalized taps
*/
private static float[] normalizeToUnityDC(double[] taps)
{
double sum = 0.0;
for(double tap : taps)
{
sum += tap;
}
float[] output = new float[taps.length];
for(int n = 0; n < taps.length; n++)
{
output[n] = (float)(taps[n] / sum);
}
return output;
}
/**
* Sinc: sin(x) / x where x = pi * f * T, the frequency response of an integrate and dump filter of length T.
*
* @param normalizedFrequency f * T (frequency divided by the symbol rate)
*/
private static double sinc(double normalizedFrequency)
{
if(normalizedFrequency == 0.0)
{
return 1.0;
}
double x = Math.PI * normalizedFrequency;
return Math.sin(x) / x;
}
/**
* Inverse sinc: x / sin(x) where x = pi * f * T.
*
* @param normalizedFrequency f * T (frequency divided by the symbol rate)
*/
private static double inverseSinc(double normalizedFrequency)
{
if(normalizedFrequency == 0.0)
{
return 1.0;
}
double x = Math.PI * normalizedFrequency;
return x / Math.sin(x);
}
/**
* Solves Ax = b using Gaussian elimination with partial pivoting. A and b are modified.
*/
private static double[] solve(double[][] a, double[] b)
{
int n = b.length;
for(int column = 0; column < n; column++)
{
int pivot = column;
for(int row = column + 1; row < n; row++)
{
if(Math.abs(a[row][column]) > Math.abs(a[pivot][column]))
{
pivot = row;
}
}
double[] swapRow = a[column];
a[column] = a[pivot];
a[pivot] = swapRow;
double swapValue = b[column];
b[column] = b[pivot];
b[pivot] = swapValue;
for(int row = column + 1; row < n; row++)
{
double factor = a[row][column] / a[column][column];
b[row] -= factor * b[column];
for(int k = column; k < n; k++)
{
a[row][k] -= factor * a[column][k];
}
}
}
double[] x = new double[n];
for(int row = n - 1; row >= 0; row--)
{
double sum = b[row];
for(int k = row + 1; k < n; k++)
{
sum -= a[row][k] * x[k];
}
x[row] = sum / a[row][row];
}
return x;
}
/**
* Evaluates the filter to determine the magnitude response (dB) at the specified frequency.
*

1
src/main/java/io/github/dsheirer/dsp/filter/design/FilterViewer.java
View File

@@ -82,7 +82,6 @@ public class FilterViewer extends Application
}
// float[] taps = FilterFactory.getHighPass(8000, 3000, filterLength, windowType);
// float[] taps = FilterFactory.getRootRaisedCosine(50000.0 / 4800.0, 20, 0.2f);
// taps = HilbertTransform.HALF_BAND_FILTER_47_TAP;

100
src/main/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliability.java Normal file
View File

@@ -0,0 +1,100 @@
/*
* *****************************************************************************
* 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.fm;
import io.github.dsheirer.dsp.filter.fir.real.IRealFilter;
import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter;
/**
* Calculates the soft symbol reliability (inverse of the relative noise variance) of each sample produced by an FM
* discriminator followed by a post-discriminator FIR filter (ie a C4FM receive filter), from the complex samples that
* are input to the discriminator.
*
* The discriminator output is the phase difference between consecutive samples: d[n] = phase(x[n]) - phase(x[n - 1]),
* so the filtered output y[n] = sum(h[k] * d[n - k]) = sum(g[k] * phase(x[n - k])) where g[k] = h[k] - h[k - 1] (with
* h[-1] = h[L] = 0). The phase noise variance of each sample is inversely proportional to the sample power, so the
* noise variance of the filtered output is proportional to sum(g[k]^2 / |x[n - k]|^2), and the reliability is the
* inverse of that sum. The sum is calculated as an FIR filter of the inverse sample power with taps g[k]^2.
*
* Testing against a P25 LSM simulcast recording found that most of the symbol errors occur in the small fraction of
* symbols with the lowest envelope, so weighting the soft decision decoding by the reliability recovers many of the
* blocks that fail unweighted soft decision decoding.
*/
public class FMDiscriminatorReliability
{
//Minimum sample power, avoiding division by zero for samples with zero power
private static final float MINIMUM_POWER = 1e-30f;
private final IRealFilter mInversePowerFilter;
/**
* Constructs an instance
* @param postDiscriminatorFilterTaps for the FIR filter applied to the FM discriminator output.
*/
public FMDiscriminatorReliability(float[] postDiscriminatorFilterTaps)
{
mInversePowerFilter = new RealFIRFilter(getInversePowerFilterTaps(postDiscriminatorFilterTaps));
}
/**
* Calculates the inverse sample power filter taps: g[k]^2 where g[k] = h[k] - h[k - 1] for k = 0 to L, with
* h[-1] = h[L] = 0.
* @param taps of the post-discriminator filter (h)
* @return inverse power filter taps, one longer than the post-discriminator filter.
*/
public static float[] getInversePowerFilterTaps(float[] taps)
{
float[] squared = new float[taps.length + 1];
for(int k = 0; k <= taps.length; k++)
{
float current = k < taps.length ? taps[k] : 0.0f;
float previous = k > 0 ? taps[k - 1] : 0.0f;
float difference = current - previous;
squared[k] = difference * difference;
}
return squared;
}
/**
* Calculates the reliability of each filtered discriminator output sample from the discriminator input samples.
* @param i inphase samples input to the FM discriminator
* @param q quadrature samples input to the FM discriminator
* @return reliabilities, parallel to the filtered discriminator output samples.
*/
public float[] calculate(float[] i, float[] q)
{
float[] inversePower = new float[i.length];
for(int x = 0; x < i.length; x++)
{
inversePower[x] = 1.0f / Math.max(i[x] * i[x] + q[x] * q[x], MINIMUM_POWER);
}
float[] reliabilities = mInversePowerFilter.filter(inversePower);
for(int x = 0; x < reliabilities.length; x++)
{
reliabilities[x] = reliabilities[x] > 0.0f ? 1.0f / reliabilities[x] : 0.0f;
}
return reliabilities;
}
}