From 6b9cd2d3173986d90e9f6e13c0ebe3e08169bcfd Mon Sep 17 00:00:00 2001 From: Denny Sheirer Date: Thu, 29 Oct 2026 04:17:47 -0400 Subject: [PATCH] #2518 P25P1 LSM Decoder enhancements. (#2523) * 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 --- .../dsheirer/buffer/FloatNativeBuffer.java | 12 +- .../dsheirer/dsp/c4fm/C4FMMessageFramer.java | 14 + .../dsp/c4fm/C4FMSymbolProcessor.java | 109 +- .../dsheirer/dsp/filter/FilterFactory.java | 675 +++++++-- .../dsp/filter/design/FilterViewer.java | 1 - .../dsp/fm/FMDiscriminatorReliability.java | 100 ++ .../edac/trellis/ViterbiDecoder_1_2_P25.java | 125 +- .../module/decode/dmr/DMRDecoder.java | 75 +- .../module/decode/nxdn/NXDNDecoder.java | 71 +- .../decode/nxdn/NXDNFilterDesigner.java | 358 ----- .../module/decode/p25/phase1/NACTracker.java | 41 +- .../decode/p25/phase1/P25P1DecoderC4FM.java | 99 +- .../decode/p25/phase1/P25P1DecoderLSM.java | 143 +- .../p25/phase1/P25P1DemodulatorC4FM.java | 1278 ----------------- .../p25/phase1/P25P1DemodulatorLSM.java | 391 ----- .../P25P1DifferentialDemodulatorLSM.java | 371 +++++ .../p25/phase1/P25P1FilterDesigner.java | 244 ---- .../decode/p25/phase1/P25P1Interleave.java | 1 + .../p25/phase1/P25P1MessageAssembler.java | 111 +- .../decode/p25/phase1/P25P1MessageFramer.java | 51 +- .../p25/phase1/P25P1MessageProcessor.java | 13 +- .../p25/phase1/P25P1SymbolProcessorC4FM.java | 12 +- .../p25/phase1/P25P1SymbolProcessorLSM.java | 131 ++ .../p25/phase1/message/hdu/HDUMessage.java | 11 + .../phase1/message/pdu/PDUMessageFactory.java | 34 +- .../pdu/block/UnconfirmedDataBlock.java | 15 +- .../p25/phase1/message/tdu/TDULCMessage.java | 10 + .../message/tsbk/TSBKMessageFactory.java | 37 +- .../decode/p25/phase2/P25P2DecoderHDQPSK.java | 53 +- .../dsheirer/source/mixer/MixerReader.java | 2 +- .../ChannelFrequencyErrorManager.java | 40 +- .../frequency/TunerFrequencyErrorManager.java | 7 + .../buffer/FloatNativeBufferTest.java | 109 ++ .../fm/FMDiscriminatorReliabilityTest.java | 95 ++ .../trellis/ViterbiDecoder_1_2_P25Test.java | 238 +++ .../decode/p25/phase1/NACTrackerTest.java | 75 + .../p25/phase1/P25P1DecoderLSMTest.java | 409 ++++++ .../P25P1DifferentialDemodulatorLSMTest.java | 353 +++++ .../p25/phase1/P25P1MessageAssemblerTest.java | 138 ++ .../phase1/P25P1SymbolProcessorLSMTest.java | 80 ++ .../message/P25P1MessageValidityTest.java | 117 ++ .../ChannelFrequencyErrorManagerTest.java | 104 ++ 42 files changed, 3586 insertions(+), 2767 deletions(-) create mode 100644 src/main/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliability.java delete mode 100644 src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNFilterDesigner.java delete mode 100644 src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorC4FM.java delete mode 100644 src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorLSM.java create mode 100644 src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSM.java delete mode 100644 src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1FilterDesigner.java create mode 100644 src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSM.java create mode 100644 src/test/java/io/github/dsheirer/buffer/FloatNativeBufferTest.java create mode 100644 src/test/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliabilityTest.java create mode 100644 src/test/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25Test.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/NACTrackerTest.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSMTest.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSMTest.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssemblerTest.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSMTest.java create mode 100644 src/test/java/io/github/dsheirer/module/decode/p25/phase1/message/P25P1MessageValidityTest.java create mode 100644 src/test/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManagerTest.java diff --git a/src/main/java/io/github/dsheirer/buffer/FloatNativeBuffer.java b/src/main/java/io/github/dsheirer/buffer/FloatNativeBuffer.java index 11cedf25..9371565b 100644 --- a/src/main/java/io/github/dsheirer/buffer/FloatNativeBuffer.java +++ b/src/main/java/io/github/dsheirer/buffer/FloatNativeBuffer.java @@ -100,8 +100,10 @@ public class FloatNativeBuffer extends AbstractNativeBuffer throw new IllegalStateException("No more samples"); } - float[] chunk = Arrays.copyOfRange(mInterleavedComplexSamples, mBufferPointer, mBufferPointer + BUFFER_SIZE * 2); - mBufferPointer += BUFFER_SIZE * 2; + //The final chunk holds the remaining samples, so that it isn't padded with zero valued samples + int end = Math.min(mBufferPointer + BUFFER_SIZE * 2, mInterleavedComplexSamples.length); + float[] chunk = Arrays.copyOfRange(mInterleavedComplexSamples, mBufferPointer, end); + mBufferPointer = end; return SampleUtils.deinterleave(chunk, getTimestamp()); } @@ -125,8 +127,10 @@ public class FloatNativeBuffer extends AbstractNativeBuffer throw new IllegalStateException("No more samples"); } - float[] chunk = Arrays.copyOfRange(mInterleavedComplexSamples, mBufferPointer, mBufferPointer + BUFFER_SIZE * 2); - mBufferPointer += BUFFER_SIZE * 2; + //The final chunk holds the remaining samples, so that it isn't padded with zero valued samples + int end = Math.min(mBufferPointer + BUFFER_SIZE * 2, mInterleavedComplexSamples.length); + float[] chunk = Arrays.copyOfRange(mInterleavedComplexSamples, mBufferPointer, end); + mBufferPointer = end; return new InterleavedComplexSamples(chunk, getTimestamp()); } diff --git a/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMMessageFramer.java b/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMMessageFramer.java index 57fd5bf9..e09790f5 100644 --- a/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMMessageFramer.java +++ b/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMMessageFramer.java @@ -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 */ diff --git a/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMSymbolProcessor.java b/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMSymbolProcessor.java index 5a7a8744..052a0ac2 100644 --- a/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMSymbolProcessor.java +++ b/src/main/java/io/github/dsheirer/dsp/c4fm/C4FMSymbolProcessor.java @@ -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); diff --git a/src/main/java/io/github/dsheirer/dsp/filter/FilterFactory.java b/src/main/java/io/github/dsheirer/dsp/filter/FilterFactory.java index f7cc527e..33434423 100644 --- a/src/main/java/io/github/dsheirer/dsp/filter/FilterFactory.java +++ b/src/main/java/io/github/dsheirer/dsp/filter/FilterFactory.java @@ -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 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. * diff --git a/src/main/java/io/github/dsheirer/dsp/filter/design/FilterViewer.java b/src/main/java/io/github/dsheirer/dsp/filter/design/FilterViewer.java index abc4a495..6303564a 100644 --- a/src/main/java/io/github/dsheirer/dsp/filter/design/FilterViewer.java +++ b/src/main/java/io/github/dsheirer/dsp/filter/design/FilterViewer.java @@ -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; diff --git a/src/main/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliability.java b/src/main/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliability.java new file mode 100644 index 00000000..d4b9c53a --- /dev/null +++ b/src/main/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliability.java @@ -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 + * **************************************************************************** + */ + +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; + } +} diff --git a/src/main/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25.java b/src/main/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25.java index 8ef4256b..ff8f2ad8 100644 --- a/src/main/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25.java +++ b/src/main/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25.java @@ -32,9 +32,21 @@ public class ViterbiDecoder_1_2_P25 extends ViterbiDecoder /** * Maximum per-symbol error magnitude (radians) applied to the soft decision branch metric. FM discriminator noise * is impulsive (clicks), and an unbounded squared error lets a single corrupted symbol dominate the path metric. - * Noise injection testing against a P25 control channel recording found a peak at 1.75 - 2.0 radians. + * Noise injection testing against a P25 control channel recording found a peak at 1.75 - 2.0 radians. Noise + * injection testing against P25 LSM recordings (phase wrapped soft symbols) also found 1.75 performed best. */ private static final double SOFT_SYMBOL_ERROR_CLIP = 1.75; + /** + * Per-symbol error magnitude clip applied when the branch metric of phase wrapped soft symbols (ie P25 LSM) is + * weighted by symbol reliabilities. The reliability weighting de-emphasizes the corrupted (low envelope) symbols + * that the clip otherwise protects against, so the weighted metric doesn't clip (PI is the maximum phase wrapped + * error). Noise injection testing against P25 LSM recordings found that no clip produced 1.6-7.0% more valid + * messages than the 1.75 clip. Soft symbols that aren't phase wrapped (ie C4FM) can have errors beyond PI from FM + * discriminator clicks, so the weighted metric for those uses the SOFT_SYMBOL_ERROR_CLIP, which testing against a + * P25 C4FM recording found 4.5% better than no clip. + */ + private static final double WEIGHTED_PHASE_WRAPPED_SOFT_SYMBOL_ERROR_CLIP = Math.PI; + private static final double TWO_PI = 2.0 * Math.PI; public ViterbiDecoder_1_2_P25() { super(2,4); @@ -68,12 +80,48 @@ public class ViterbiDecoder_1_2_P25 extends ViterbiDecoder * @return decoded 96-bit message */ public CorrectedBinaryMessage decodeSoft(float[] softSymbols) + { + return decodeSoft(softSymbols, false); + } + + /** + * Soft decision Viterbi decoding of a deinterleaved 1/2 rate TCM encoded P25 block. See decodeSoft(float[]). + * + * @param softSymbols 98 deinterleaved soft symbols (radians, ideal values +/- PI/4 and +/- 3*PI/4). + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie differentially demodulated P25 LSM) + * where the soft symbol error is measured around the circle so that a soft symbol near +/- PI is near both the + * +3 and -3 symbols. False for C4FM where soft symbols are not wrapped and a soft symbol beyond +/- PI favors the + * outer symbol of the same polarity. + * @return decoded 96-bit message + */ + public CorrectedBinaryMessage decodeSoft(float[] softSymbols, boolean phaseWrapped) + { + return decodeSoft(softSymbols, null, phaseWrapped); + } + + /** + * Soft decision Viterbi decoding of a deinterleaved 1/2 rate TCM encoded P25 block with each symbol's branch metric + * weighted by the symbol's reliability (ie the inverse of the soft symbol noise variance). See decodeSoft(float[]). + * + * The weights only need to be proportional to the symbol reliabilities, since they're normalized to unity mean + * across the block. + * + * @param softSymbols 98 deinterleaved soft symbols (radians, ideal values +/- PI/4 and +/- 3*PI/4). + * @param reliabilities for each deinterleaved soft symbol, or null to decode without weighting. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie differentially demodulated P25 LSM) + * @return decoded 96-bit message + */ + public CorrectedBinaryMessage decodeSoft(float[] softSymbols, float[] reliabilities, boolean phaseWrapped) { if(softSymbols == null || softSymbols.length % 2 != 0 || softSymbols.length < 4) { throw new IllegalArgumentException("Soft symbols must contain an integral number of symbol pairs"); } + double[] weights = getWeights(reliabilities, softSymbols.length); + double clip = weights != null && phaseWrapped ? WEIGHTED_PHASE_WRAPPED_SOFT_SYMBOL_ERROR_CLIP : + SOFT_SYMBOL_ERROR_CLIP; + int steps = softSymbols.length / 2; int states = getInputValueCount(); double[] metrics = new double[states]; @@ -92,9 +140,18 @@ public class ViterbiDecoder_1_2_P25 extends ViterbiDecoder for(int output = 0; output < branchMetrics.length; output++) { - double errorFirst = Math.min(Math.abs(first - getIdealPhase(output >> 2)), SOFT_SYMBOL_ERROR_CLIP); - double errorSecond = Math.min(Math.abs(second - getIdealPhase(output & 3)), SOFT_SYMBOL_ERROR_CLIP); - branchMetrics[output] = errorFirst * errorFirst + errorSecond * errorSecond; + double errorFirst = Math.min(getError(first, getIdealPhase(output >> 2), phaseWrapped), clip); + double errorSecond = Math.min(getError(second, getIdealPhase(output & 3), phaseWrapped), clip); + + if(weights != null) + { + branchMetrics[output] = weights[2 * step] * errorFirst * errorFirst + + weights[2 * step + 1] * errorSecond * errorSecond; + } + else + { + branchMetrics[output] = errorFirst * errorFirst + errorSecond * errorSecond; + } } java.util.Arrays.fill(nextMetrics, Double.MAX_VALUE); @@ -170,6 +227,66 @@ public class ViterbiDecoder_1_2_P25 extends ViterbiDecoder return message; } + /** + * Absolute error between the soft symbol and the ideal symbol phase. + * @param softSymbol received (radians) + * @param idealPhase of the candidate symbol (radians) + * @param phaseWrapped to measure the error around the circle (+/- PI) + * @return absolute error in radians + */ + private static double getError(float softSymbol, float idealPhase, boolean phaseWrapped) + { + double error = softSymbol - idealPhase; + + if(phaseWrapped) + { + error -= TWO_PI * Math.round(error / TWO_PI); + } + + return Math.abs(error); + } + + /** + * Normalizes the symbol reliabilities to unity mean for use as branch metric weights. + * @param reliabilities for each soft symbol, or null. + * @param length of the soft symbol array. + * @return weights or null when the reliabilities are unavailable or unusable (ie all zero). + */ + private static double[] getWeights(float[] reliabilities, int length) + { + if(reliabilities == null || reliabilities.length != length) + { + return null; + } + + double sum = 0.0; + + for(float reliability : reliabilities) + { + if(!(reliability >= 0.0f) || Float.isInfinite(reliability)) //Rejects NaN + { + return null; + } + + sum += reliability; + } + + if(sum <= 0.0) + { + return null; + } + + double scale = length / sum; + double[] weights = new double[length]; + + for(int x = 0; x < length; x++) + { + weights[x] = reliabilities[x] * scale; + } + + return weights; + } + /** * Ideal soft symbol phase for a 2-bit dibit value where the MSB is the first transmitted bit. */ diff --git a/src/main/java/io/github/dsheirer/module/decode/dmr/DMRDecoder.java b/src/main/java/io/github/dsheirer/module/decode/dmr/DMRDecoder.java index 754979c7..7d07ccbe 100644 --- a/src/main/java/io/github/dsheirer/module/decode/dmr/DMRDecoder.java +++ b/src/main/java/io/github/dsheirer/module/decode/dmr/DMRDecoder.java @@ -23,7 +23,6 @@ import io.github.dsheirer.alias.AliasList; import io.github.dsheirer.dsp.filter.FilterFactory; import io.github.dsheirer.dsp.filter.decimate.DecimationFilterFactory; import io.github.dsheirer.dsp.filter.decimate.IRealDecimationFilter; -import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification; import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter; import io.github.dsheirer.dsp.fm.FmDemodulatorFactory; @@ -40,7 +39,6 @@ import io.github.dsheirer.module.decode.dmr.message.data.lc.LCMessage; import io.github.dsheirer.module.decode.dmr.message.data.lc.full.FullLCMessage; import io.github.dsheirer.module.decode.dmr.message.data.lc.shorty.ShortLCMessage; import io.github.dsheirer.module.decode.dmr.message.data.terminator.Terminator; -import io.github.dsheirer.module.decode.nxdn.NXDNFilterDesigner; import io.github.dsheirer.preference.UserPreferences; import io.github.dsheirer.sample.Listener; import io.github.dsheirer.sample.buffer.IByteBufferProvider; @@ -54,8 +52,6 @@ import java.io.IOException; import java.nio.ByteBuffer; import java.text.DecimalFormat; import java.util.Iterator; -import java.util.Map; -import java.util.concurrent.ConcurrentHashMap; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -88,7 +84,8 @@ public class DMRDecoder extends FeedbackDecoder implements IByteBufferProvider, //DMR transmit filter is a root raised cosine with 0.2 roll-off, so the receive filter is a matched RRC filter. private static final double RECEIVE_FILTER_ALPHA = 0.2; private static final int RECEIVE_FILTER_SYMBOL_SPAN = 12; - private static final Map BASEBAND_FILTERS = new ConcurrentHashMap<>(); + private static final double BASEBAND_FILTER_PASS_BAND = 5100.0; //Hertz + private static final double BASEBAND_FILTER_STOP_BAND = 6500.0; //Hertz private IDemodulator mDemodulator; private final DMRMessageFramer mMessageFramer; private final DMRSymbolProcessor mSymbolProcessor; @@ -155,17 +152,13 @@ public class DMRDecoder extends FeedbackDecoder implements IByteBufferProvider, setDecimatedSampleRate(decimatedSampleRate); //Receive filter: root raised cosine matched to the DMR transmit filter, normalized to unity gain at DC. - double[] receiveFilterTaps = NXDNFilterDesigner.designRootRaisedCosine(decimatedSampleRate, SYMBOL_RATE, - RECEIVE_FILTER_ALPHA, RECEIVE_FILTER_SYMBOL_SPAN); - float[] taps = new float[receiveFilterTaps.length]; - for(int x = 0; x < receiveFilterTaps.length; x++) - { - taps[x] = (float)receiveFilterTaps[x]; - } - mReceiveFilter = new RealFIRFilter(taps); + mReceiveFilter = new RealFIRFilter(FilterFactory.getRootRaisedCosineFilter(decimatedSampleRate, SYMBOL_RATE, + RECEIVE_FILTER_ALPHA, RECEIVE_FILTER_SYMBOL_SPAN)); - mBasebandFilterI = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); - mBasebandFilterQ = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); + float[] basebandTaps = FilterFactory.getRemezLowPass(decimatedSampleRate, BASEBAND_FILTER_PASS_BAND, + BASEBAND_FILTER_STOP_BAND); + mBasebandFilterI = FilterFactory.getRealFilter(basebandTaps); + mBasebandFilterQ = FilterFactory.getRealFilter(basebandTaps); mDemodulator = FmDemodulatorFactory.getFmDemodulator(); mSymbolProcessor.setSamplesPerSymbol(samplesPerSymbol); @@ -197,51 +190,6 @@ public class DMRDecoder extends FeedbackDecoder implements IByteBufferProvider, mSymbolProcessor.process(demodulated); } - /** - * Constructs a baseband filter for this decoder using the current sample rate - */ - private float[] getBasebandFilter(double sampleRate) - { - //Filter cache is shared across decoder instances that may be started from different threads - float[] cached = BASEBAND_FILTERS.get(sampleRate); - - if(cached != null) - { - return cached; - } - - FIRFilterSpecification specification = FIRFilterSpecification - .lowPassBuilder() - .sampleRate(sampleRate) - .passBandCutoff(5100) - .passBandAmplitude(1.0).passBandRipple(0.01) //.01 - .stopBandAmplitude(0.0).stopBandStart(6500) //6500 - .stopBandRipple(0.01).build(); - - float[] coefficients = null; - - try - { - coefficients = FilterFactory.getTaps(specification); - - if(coefficients != null) - { - BASEBAND_FILTERS.put(sampleRate, coefficients); - } - } - catch(Exception e) //FilterDesignException or runtime failures from the Remez designer - { - LOGGER.error("Error designing DMR baseband filter for sample rate [" + sampleRate + "]", e); - } - - if(coefficients == null) - { - throw new IllegalStateException("Unable to design low pass filter for sample rate [" + sampleRate + "]"); - } - - return coefficients; - } - /** * Implements the IByteBufferProvider interface - delegates to the symbol processor */ @@ -306,9 +254,14 @@ public class DMRDecoder extends FeedbackDecoder implements IByteBufferProvider, switch(sourceEvent.getEvent()) { case NOTIFICATION_FREQUENCY_CHANGE: - case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: mSymbolProcessor.resetBalance(); break; + case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: + //Tuner PPM changes are small (auto-PPM steps are <= 10 Hz) and the symbol processor balance tracks + //the change. Resetting drops the symbol processor synchronization, and the symbol timing tracking + //with it, until the next sync detection, which costs messages on channels with long sync intervals + //(ie P25 voice). + break; case NOTIFICATION_SAMPLE_RATE_CHANGE: setSampleRate(sourceEvent.getValue().doubleValue()); break; diff --git a/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNDecoder.java b/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNDecoder.java index 6c2cf0d4..7c0e66a8 100644 --- a/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNDecoder.java +++ b/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNDecoder.java @@ -23,7 +23,6 @@ import com.google.common.eventbus.Subscribe; import io.github.dsheirer.dsp.filter.FilterFactory; import io.github.dsheirer.dsp.filter.decimate.DecimationFilterFactory; import io.github.dsheirer.dsp.filter.decimate.IRealDecimationFilter; -import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification; import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter; import io.github.dsheirer.dsp.fm.FmDemodulatorFactory; @@ -45,8 +44,6 @@ import java.io.File; import java.io.IOException; import java.nio.ByteBuffer; import java.util.Iterator; -import java.util.Map; -import java.util.concurrent.ConcurrentHashMap; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -57,7 +54,8 @@ public class NXDNDecoder extends FeedbackDecoder implements IByteBufferProvider, ISourceEventListener, Listener { private static final Logger LOGGER = LoggerFactory.getLogger(NXDNDecoder.class); - private static final Map BASEBAND_FILTERS = new ConcurrentHashMap<>(); + //NXDN transmit filter root raised cosine roll-off + private static final double RECEIVE_FILTER_ALPHA = 0.2; private final DecodeConfigNXDN mConfig; private final NXDNSymbolProcessor mSymbolProcessor; private final NXDNMessageFramer mMessageFramer; @@ -141,13 +139,13 @@ public class NXDNDecoder extends FeedbackDecoder implements IByteBufferProvider, int rrcSymbolSpan = 12; int inverseSyncSymbolSpan = 4; //Create an RRC pulse shaping filter convolved with an inverse sync (ie restoration) filter. - float[] taps = NXDNFilterDesigner.designNXDN(decimatedSampleRate, mConfig.getTransmissionMode().getSymbolRate(), - rrcSymbolSpan, inverseSyncSymbolSpan); -// int symbolLength = 26; -// float rrcAlpha = 0.2f; -// taps = FilterFactory.getRRC(samplesPerSymbol, symbolLength, rrcAlpha); + float[] taps = FilterFactory.getRootRaisedCosineInverseSincFilter(decimatedSampleRate, + mConfig.getTransmissionMode().getSymbolRate(), RECEIVE_FILTER_ALPHA, rrcSymbolSpan, + inverseSyncSymbolSpan); mPulseShapingFilter = new RealFIRFilter(taps); - float[] basebandTaps = getBasebandFilter(decimatedSampleRate, mConfig); + ChannelSpecification channel = mConfig.getChannelSpecification(); + float[] basebandTaps = FilterFactory.getRemezLowPass(decimatedSampleRate, channel.getPassFrequency(), + channel.getStopFrequency()); mBasebandFilterI = FilterFactory.getRealFilter(basebandTaps); mBasebandFilterQ = FilterFactory.getRealFilter(basebandTaps); mDemodulator = FmDemodulatorFactory.getFmDemodulator(); @@ -178,52 +176,6 @@ public class NXDNDecoder extends FeedbackDecoder implements IByteBufferProvider, mSymbolProcessor.process(demodulated); } - /** - * Constructs a baseband filter for this decoder using the current sample rate - */ - private float[] getBasebandFilter(double sampleRate, DecodeConfigNXDN config) - { - String mapKey = config.getTransmissionMode().name() + sampleRate; - - //Filter cache is shared across decoder instances that may be started from different threads - float[] coefficients = BASEBAND_FILTERS.get(mapKey); - - if(coefficients == null) - { - ChannelSpecification channel = config.getChannelSpecification(); - - FIRFilterSpecification specification = FIRFilterSpecification - .lowPassBuilder() - .sampleRate(sampleRate) - .passBandCutoff(channel.getPassFrequency()) - .passBandAmplitude(1.0).passBandRipple(0.01) - .stopBandAmplitude(0.0).stopBandStart(channel.getStopFrequency()) - .stopBandRipple(0.01).build(); - - try - { - coefficients = FilterFactory.getTaps(specification); - - if(coefficients != null) - { - BASEBAND_FILTERS.put(mapKey, coefficients); - } - } - catch(Exception e) //FilterDesignException or runtime failures from the Remez designer - { - LOGGER.error("Error designing NXDN " + config.getTransmissionMode() + " baseband filter for sample rate [" + - sampleRate + "]", e); - } - } - - if(coefficients == null) - { - throw new IllegalStateException("Unable to design low pass filter for sample rate [" + sampleRate + "]"); - } - - return coefficients; - } - /** * Implements the IByteBufferProvider interface - delegates to the symbol processor */ @@ -272,9 +224,14 @@ public class NXDNDecoder extends FeedbackDecoder implements IByteBufferProvider, switch(sourceEvent.getEvent()) { case NOTIFICATION_FREQUENCY_CHANGE: - case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: mSymbolProcessor.resetBalance(); break; + case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: + //Tuner PPM changes are small (auto-PPM steps are <= 10 Hz) and the symbol processor balance tracks the + //change. Resetting drops the symbol processor synchronization, and the symbol timing tracking with it, + //until the next sync detection, which costs messages on channels with long sync intervals (ie P25 + //voice). + break; case NOTIFICATION_SAMPLE_RATE_CHANGE: setSampleRate(sourceEvent.getValue().doubleValue()); break; diff --git a/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNFilterDesigner.java b/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNFilterDesigner.java deleted file mode 100644 index d8987db2..00000000 --- a/src/main/java/io/github/dsheirer/module/decode/nxdn/NXDNFilterDesigner.java +++ /dev/null @@ -1,358 +0,0 @@ -package io.github.dsheirer.module.decode.nxdn; - -/** - * Designs an inverse sinc (x / sin(x)) FIR filter for NXDN 4FSK receivers. - * - * The NXDN transmitter applies a sinc response, sin(pi*f*T)/(pi*f*T), to the symbol stream along with the root - * raised cosine filter. The receiver cascades this inverse sinc filter, H(f) = (pi*f*T)/sin(pi*f*T), with its own - * root raised cosine filter to restore an overall raised cosine (Nyquist) response. - * - * The inverse sinc grows without bound as f approaches the symbol rate, so the response is only matched over the - * occupied bandwidth of the signal, |f| <= (1 + alpha) * symbolRate / 2. Beyond that edge the cascaded RRC filter - * suppresses the signal, so that region is lightly weighted toward a flat gain equal to the band-edge value, which - * keeps the filter from boosting out-of-band noise. - * - * The taps are a weighted least-squares fit of a linear phase (symmetric, odd length) FIR filter to that target - * response on a dense frequency grid. - */ -public class NXDNFilterDesigner -{ - /** NXDN root raised cosine roll-off factor. */ - public static final double NXDN_ALPHA = 0.2; - - /** Default filter length in symbol periods. */ - private static final int DEFAULT_SYMBOL_SPAN = 4; - - /** Least-squares weight applied beyond the band edge, relative to a weight of 1.0 in band. */ - private static final double STOPBAND_WEIGHT = 0.01; - - /** Frequency grid points per filter coefficient used for the least-squares fit. */ - private static final int GRID_DENSITY = 16; - - /** - * Designs an inverse sinc filter using the NXDN roll-off factor (0.2) to define the matched bandwidth. - * - * @param sampleRate in Hertz - * @param symbolRate in symbols per second (2400 for NXDN48, 4800 for NXDN96) - * @return filter taps, odd length, linear phase, unity gain at DC - */ - public static float[] design(double sampleRate, double symbolRate) - { - return design(sampleRate, symbolRate, NXDN_ALPHA, DEFAULT_SYMBOL_SPAN); - } - - /** - * Designs an inverse sinc filter. - * - * @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, linear phase, unity gain at DC - */ - public static float[] design(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 nyquist = sampleRate / 2.0; - double bandEdge = (1.0 + alpha) * symbolRate / 2.0; - - if(bandEdge >= nyquist) - { - throw new IllegalArgumentException("Sample rate [" + sampleRate + - "] is too low for symbol rate [" + symbolRate + "] - band edge exceeds Nyquist"); - } - - //Odd tap count so the filter has an integer group delay of (length - 1) / 2 samples - int length = (int)Math.ceil(symbolSpan * sampleRate / symbolRate); - if(length % 2 == 0) - { - length++; - } - length = Math.max(length, 9); - - int middle = (length - 1) / 2; - int coefficientCount = middle + 1; - double symbolPeriod = 1.0 / symbolRate; - double edgeGain = inverseSinc(bandEdge * symbolPeriod); - - //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 from DC to Nyquist. - double[][] normal = new double[coefficientCount][coefficientCount]; - double[] rhs = new double[coefficientCount]; - double[] basis = new double[coefficientCount]; - int gridSize = GRID_DENSITY * coefficientCount; - - for(int g = 0; g <= gridSize; g++) - { - double frequency = nyquist * g / gridSize; - boolean inBand = frequency <= bandEdge; - double desired = inBand ? inverseSinc(frequency * symbolPeriod) : edgeGain; - double weight = inBand ? 1.0 : STOPBAND_WEIGHT; - - 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 and normalize to unity gain at DC - double[] taps = new double[length]; - double sum = 0.0; - - for(int m = 0; m < coefficientCount; m++) - { - taps[middle - m] = coefficients[m]; - taps[middle + m] = coefficients[m]; - } - - for(double tap : taps) - { - sum += tap; - } - - float[] output = new float[length]; - for(int n = 0; n < length; n++) - { - output[n] = (float)(taps[n] / sum); - } - - return output; - } - - /** Default root raised cosine filter length in symbol periods for the combined NXDN receive filter. */ - private static final int DEFAULT_RRC_SYMBOL_SPAN = 8; - - /** - * Designs the complete NXDN receive filter: a root raised cosine filter (alpha 0.2) convolved with the inverse - * sinc filter. Cascaded with the NXDN transmit filter (RRC and sinc), 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 (2400 for NXDN48, 4800 for NXDN96) - * @return filter taps, odd length, linear phase, unity gain at DC - */ - public static float[] designNXDN(double sampleRate, double symbolRate) - { - return designNXDN(sampleRate, symbolRate, DEFAULT_RRC_SYMBOL_SPAN, DEFAULT_SYMBOL_SPAN); - } - - /** - * Designs the complete NXDN receive filter: a root raised cosine filter (alpha 0.2) convolved with the inverse - * sinc filter. - * - * @param sampleRate in Hertz - * @param symbolRate in symbols per second - * @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[] designNXDN(double sampleRate, double symbolRate, int rrcSymbolSpan, int inverseSincSymbolSpan) - { - double[] rrc = designRootRaisedCosine(sampleRate, symbolRate, NXDN_ALPHA, rrcSymbolSpan); - float[] inverseSinc = design(sampleRate, symbolRate, NXDN_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 - double sum = 0.0; - for(double tap : combined) - { - sum += tap; - } - - float[] output = new float[combined.length]; - for(int n = 0; n < combined.length; n++) - { - output[n] = (float)(combined[n] / sum); - } - - return output; - } - - /** - * 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 double[] designRootRaisedCosine(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 = (int)Math.ceil(symbolSpan * sampleRate / symbolRate); - if(length % 2 == 0) - { - length++; - } - - 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; - } - - /** - * 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; - } -} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/NACTracker.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/NACTracker.java index be27ed68..1cd3b414 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/NACTracker.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/NACTracker.java @@ -21,7 +21,6 @@ package io.github.dsheirer.module.decode.p25.phase1; import java.util.ArrayList; import java.util.Comparator; -import java.util.Date; import java.util.HashMap; import java.util.List; import java.util.Map; @@ -32,10 +31,14 @@ import java.util.Map; public class NACTracker { private static final Comparator COUNT_DESCENDING_COMPARATOR = new CountDescendingSorter(); - private static final Comparator TIME_ASCENDING_COMPARATOR = new TimeAscendingSorter(); + private static final Comparator OBSERVATION_ASCENDING_COMPARATOR = new ObservationAscendingSorter(); private final Map mTrackerMap = new HashMap<>(); private static final int MAX_TRACKER_COUNT = 3; private static final int MIN_OBSERVATION_THRESHOLD = 3; + //Observation sequence number, used instead of the system clock to order the trackers by their last observation. + //Recordings decode faster than real time, so many observations share the same millisecond, and the clock would + //make the eviction choice (and the decoded output) depend on the processing speed. + private long mObservationCount = 0; /** * Removes all tracked NAC values. Invoke this method after an extended loss of sync. @@ -47,13 +50,14 @@ public class NACTracker /** * Tracks the NAC value. Each time the NAC value is observed, the observation count is incremented so that the - * dominant NAC value will have the highest observation count. Additionally, the last updated timestamp for the - * tracker is updated to now. When the tracker count exceeds the max tracker threshold, the tracker with the - * oldest updated timestamp is thrown out. + * dominant NAC value will have the highest observation count. Additionally, the tracker's last observation is + * updated to this observation. When the tracker count exceeds the max tracker threshold, the tracker with the + * oldest last observation is thrown out. * @param nac to track */ public void track(int nac) { + mObservationCount++; Tracker tracker = mTrackerMap.get(nac); if(tracker == null) @@ -63,7 +67,7 @@ public class NACTracker if(mTrackerMap.size() > MAX_TRACKER_COUNT) { List trackerList = new ArrayList<>(mTrackerMap.values()); - trackerList.sort(TIME_ASCENDING_COMPARATOR); + trackerList.sort(OBSERVATION_ASCENDING_COMPARATOR); Tracker oldest = trackerList.getFirst(); mTrackerMap.remove(oldest.nac()); } @@ -104,11 +108,12 @@ public class NACTracker private class Tracker { private int mNAC; - private long mTimestamp = System.currentTimeMillis(); + private long mLastObservation = mObservationCount; private int mCount = 1; /** - * Constructs an instance to track the NAC. Note: count is set to 1 and timestamp is set to now. + * Constructs an instance to track the NAC. Note: count is set to 1 and the last observation is set to the + * current observation. * @param nac to track. */ private Tracker(int nac) @@ -126,7 +131,7 @@ public class NACTracker } /** - * Increments the observation count for this tracker and updates the timestamp. + * Increments the observation count for this tracker and updates the last observation. */ public void increment() { @@ -138,16 +143,16 @@ public class NACTracker mCount = Integer.MAX_VALUE; } - mTimestamp = System.currentTimeMillis(); + mLastObservation = mObservationCount; } /** - * Last observation timestamp. - * @return timestamp. + * Observation sequence number of the last observation of this tracker's NAC. + * @return last observation. */ - public long timestamp() + public long lastObservation() { - return mTimestamp; + return mLastObservation; } /** @@ -165,7 +170,7 @@ public class NACTracker @Override public String toString() { - return "NAC [" + mNAC + "] COUNT [" + mCount + "] " + new Date(mTimestamp); + return "NAC [" + mNAC + "] COUNT [" + mCount + "] LAST OBSERVATION [" + mLastObservation + "]"; } } @@ -183,14 +188,14 @@ public class NACTracker } /** - * Sorts the tracked NACs in time order where index 0 is the oldest tracked value. + * Sorts the tracked NACs in last observation order where index 0 is the least recently observed value. */ - public static class TimeAscendingSorter implements Comparator + public static class ObservationAscendingSorter implements Comparator { @Override public int compare(Tracker o1, Tracker o2) { - return Long.compare(o1.timestamp(), o2.timestamp()); + return Long.compare(o1.lastObservation(), o2.lastObservation()); } } } diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderC4FM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderC4FM.java index d143bfde..6e96451d 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderC4FM.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderC4FM.java @@ -23,9 +23,9 @@ import io.github.dsheirer.alias.AliasList; import io.github.dsheirer.dsp.filter.FilterFactory; import io.github.dsheirer.dsp.filter.decimate.DecimationFilterFactory; import io.github.dsheirer.dsp.filter.decimate.IRealDecimationFilter; -import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification; import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter; +import io.github.dsheirer.dsp.fm.FMDiscriminatorReliability; import io.github.dsheirer.dsp.fm.FmDemodulatorFactory; import io.github.dsheirer.dsp.fm.IDemodulator; import io.github.dsheirer.message.IMessage; @@ -48,8 +48,6 @@ import java.io.IOException; import java.nio.ByteBuffer; import java.text.DecimalFormat; import java.util.Iterator; -import java.util.Map; -import java.util.concurrent.ConcurrentHashMap; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -62,6 +60,9 @@ import org.slf4j.LoggerFactory; * * As a child of the FeedbackDecoder, this decoder provides periodic PLL measurements to the tuner for automatic PPM * correction. It also provides a stream of demodulated soft symbols (in radians) for display to the user. + * + * Soft decision decoding is weighted by each soft symbol's reliability (see FMDiscriminatorReliability), calculated + * from the power of the complex samples spanned by the receive filter. */ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProvider, IComplexSamplesListener, ISourceEventListener, Listener @@ -69,11 +70,14 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv private static final Logger LOGGER = LoggerFactory.getLogger(P25P1DecoderC4FM.class); private static final DecimalFormat DECIMAL_FORMAT = new DecimalFormat("#.##"); private static final int SYMBOL_RATE = 4800; + //TIA-102.BAAA C4FM transmit raised cosine filter roll-off + private static final double RECEIVE_FILTER_ALPHA = 0.2; //Band limited integrate and dump receive filter settings. Noise injection testing against P25 C4FM recordings: //stopband starts of 3300 - 4000 Hz perform similarly, a 5-symbol span is too short, and 12 symbols adds nothing. private static final double RECEIVE_FILTER_STOP_BAND_START = 3600.0; //Hertz private static final int RECEIVE_FILTER_SYMBOL_SPAN = 8; - private static final Map BASEBAND_FILTERS = new ConcurrentHashMap<>(); + private static final double BASEBAND_FILTER_PASS_BAND = 5200.0; //Hertz + private static final double BASEBAND_FILTER_STOP_BAND = 6500.0; //Hertz private final P25P1SymbolProcessorC4FM mSymbolProcessor; private final P25P1MessageFramer mMessageFramer = new P25P1MessageFramer(); @@ -84,6 +88,7 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv private IRealFilter mBasebandFilterI; private IRealFilter mBasebandFilterQ; private IRealFilter mReceiveFilter; + private FMDiscriminatorReliability mReliability; @Override public DecoderType getDecoderType() @@ -135,12 +140,15 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv //Receive filter: integrate and dump (sinc) response across the signal band with a stopband above the signal band //to reject FM discriminator noise that an integrate and dump filter's sidelobes would otherwise pass. - float[] receiveFilterTaps = P25P1FilterDesigner.designReceiveFilter(decimatedSampleRate, SYMBOL_RATE, - RECEIVE_FILTER_STOP_BAND_START, RECEIVE_FILTER_SYMBOL_SPAN); + float[] receiveFilterTaps = FilterFactory.getBandLimitedIntegrateAndDump(decimatedSampleRate, SYMBOL_RATE, + RECEIVE_FILTER_ALPHA, RECEIVE_FILTER_STOP_BAND_START, RECEIVE_FILTER_SYMBOL_SPAN); mReceiveFilter = new RealFIRFilter(receiveFilterTaps); + mReliability = new FMDiscriminatorReliability(receiveFilterTaps); - mBasebandFilterI = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); - mBasebandFilterQ = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); + float[] basebandTaps = FilterFactory.getRemezLowPass(decimatedSampleRate, BASEBAND_FILTER_PASS_BAND, + BASEBAND_FILTER_STOP_BAND); + mBasebandFilterI = FilterFactory.getRealFilter(basebandTaps); + mBasebandFilterQ = FilterFactory.getRealFilter(basebandTaps); mDemodulator = FmDemodulatorFactory.getFmDemodulator(); mSymbolProcessor.setSamplesPerSymbol(samplesPerSymbol); @@ -151,6 +159,24 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv mMessageProcessor.setMessageListener(getMessageListener()); } + /** + * Converts the soft symbol reliabilities to soft decision decoding weights. Noise injection testing against a P25 + * C4FM control channel recording found that the square root of the reliability produced 1.6% more valid messages + * at 9 dB carrier to noise (in-band) than the reliability, and 6.3% more than unweighted decoding. The P25 LSM + * decoder uses the reliability, which tested 5% better than the square root against LSM recordings. + * @param reliabilities to convert (modified in place) + * @return weights + */ + private static float[] getWeights(float[] reliabilities) + { + for(int x = 0; x < reliabilities.length; x++) + { + reliabilities[x] = (float)Math.sqrt(reliabilities[x]); + } + + return reliabilities; + } + /** * Primary method for processing incoming complex sample buffers * @param samples containing channelized complex samples @@ -170,53 +196,9 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv float[] demodulated = mDemodulator.demodulate(i, q); demodulated = mReceiveFilter.filter(demodulated); - //Process demodulated samples into symbols and apply message sync detection and framing. - mSymbolProcessor.process(demodulated); - } - - /** - * Constructs a baseband filter for this decoder using the current sample rate - */ - private float[] getBasebandFilter(double sampleRate) - { - //Filter cache is shared across decoder instances that may be started from different threads - float[] cached = BASEBAND_FILTERS.get(sampleRate); - - if(cached != null) - { - return cached; - } - - FIRFilterSpecification specification = FIRFilterSpecification - .lowPassBuilder() - .sampleRate(sampleRate) - .passBandCutoff(5200) - .passBandAmplitude(1.0).passBandRipple(0.01) //.01 - .stopBandAmplitude(0.0).stopBandStart(6500) //6500 - .stopBandRipple(0.01).build(); - - float[] coefficients = null; - - try - { - coefficients = FilterFactory.getTaps(specification); - - if(coefficients != null) - { - BASEBAND_FILTERS.put(sampleRate, coefficients); - } - } - catch(Exception e) //FilterDesignException or runtime failures from the Remez designer - { - LOGGER.error("Error designing P25 Phase 1 C4FM baseband filter for sample rate [" + sampleRate + "]", e); - } - - if(coefficients == null) - { - throw new IllegalStateException("Unable to design low pass filter for sample rate [" + sampleRate + "]"); - } - - return coefficients; + //Process demodulated samples into symbols and apply message sync detection and framing. The soft symbol + //reliabilities weight the soft decision decoding. + mSymbolProcessor.process(demodulated, getWeights(mReliability.calculate(i, q))); } /** @@ -275,9 +257,14 @@ public class P25P1DecoderC4FM extends FeedbackDecoder implements IByteBufferProv switch(sourceEvent.getEvent()) { case NOTIFICATION_FREQUENCY_CHANGE: - case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: mSymbolProcessor.resetBalance(); break; + case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: + //Tuner PPM changes are small (auto-PPM steps are <= 10 Hz) and the symbol processor balance tracks the + //change. Resetting drops the symbol processor synchronization, and the symbol timing tracking with it, + //until the next sync detection, which costs messages on channels with long sync intervals (ie P25 + //voice). + break; case NOTIFICATION_SAMPLE_RATE_CHANGE: setSampleRate(sourceEvent.getValue().doubleValue()); break; diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSM.java index e6d99a9b..290fb33c 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSM.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSM.java @@ -23,9 +23,7 @@ import io.github.dsheirer.alias.AliasList; import io.github.dsheirer.dsp.filter.FilterFactory; import io.github.dsheirer.dsp.filter.decimate.DecimationFilterFactory; import io.github.dsheirer.dsp.filter.decimate.IRealDecimationFilter; -import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification; import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; -import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter; import io.github.dsheirer.message.IMessage; import io.github.dsheirer.message.SyncLossMessage; import io.github.dsheirer.module.decode.DecoderType; @@ -44,17 +42,21 @@ import java.io.File; import java.io.IOException; import java.nio.ByteBuffer; import java.text.DecimalFormat; -import java.util.HashMap; import java.util.Iterator; -import java.util.Map; import org.slf4j.Logger; import org.slf4j.LoggerFactory; /** - * APCO25 Phase 1 Linear Simulcast Modulation (LSM) decoder. Decimates incoming sample buffers to as close as possible - * to 25 kHz for ~5 samples per symbol. Employs baseband and pulse shaping filters. Incorporates an demodulator to - * process complex baseband (I/Q) sample stream into soft symbols with soft sync detection and message framing. A - * registered message listener receives the detected and framed messages. + * APCO25 Phase 1 Linear Simulcast Modulation (LSM) decoder. Decimates incoming sample buffers to 4-8 samples per + * symbol. Employs a flat baseband (channel) filter without a matched (root raised cosine) filter since the differential + * demodulator is a limiter-discriminator receiver where narrowing the bandwidth ahead of the discriminator distorts the + * phase trajectory and degrades performance (see BASEBAND_FILTER_PASS_BAND). The differential demodulator + * (mDemodulator) removes the coarse carrier offset and differentially demodulates the filtered complex baseband (I/Q) + * sample stream into the phase change across each symbol period, in the same radians per sample domain as the C4FM + * decoder. This allows the LSM decoder to share the C4FM symbol processor (mSymbolProcessor) that performs sync + * detection, timing correction, equalization and balance (PLL) residual mistune correction. Message framer performs + * message framing and message creation, using the demodulator's soft symbol reliabilities to weight soft decision + * decoding. A registered message listener receives the detected and framed messages. * * As a child of the FeedbackDecoder, this decoder provides periodic PLL measurements to the tuner for automatic PPM * correction. It also provides a stream of demodulated soft symbols (in radians) for display to the user. @@ -64,18 +66,23 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi { private static final Logger LOGGER = LoggerFactory.getLogger(P25P1DecoderLSM.class); private static final DecimalFormat DECIMAL_FORMAT = new DecimalFormat("#.##"); - private static final Map BASEBAND_FILTERS = new HashMap<>(); + //Baseband filter settings. Testing against LSM recordings: a 7250 Hz pass band passes adjacent channel + //interference and a matched (4800 baud root raised cosine) filter increases symbol errors. Without injected noise, + //flat pass bands from 4800 to 5600 Hz perform best, while narrower pass bands perform better with injected noise. + //A 4400 Hz pass band balances both: relative to 5200 Hz, 0.3% fewer valid messages without injected noise and + //25-50% more valid messages at 3-0 dB carrier to noise. + private static final double BASEBAND_FILTER_PASS_BAND = 4400.0; //Hertz + private static final double BASEBAND_FILTER_STOP_BAND = 5600.0; //Hertz private static final int SYMBOL_RATE = 4800; - private final P25P1DemodulatorLSM mDemodulator; + private final P25P1SymbolProcessorLSM mSymbolProcessor; + private final P25P1DifferentialDemodulatorLSM mDemodulator = new P25P1DifferentialDemodulatorLSM(); private final P25P1MessageFramer mMessageFramer = new P25P1MessageFramer(); private final P25P1MessageProcessor mMessageProcessor = new P25P1MessageProcessor(); private IRealDecimationFilter mDecimationFilterI; private IRealDecimationFilter mDecimationFilterQ; private IRealFilter mBasebandFilterI; private IRealFilter mBasebandFilterQ; - private IRealFilter mPulseShapingFilterI; - private IRealFilter mPulseShapingFilterQ; @Override public DecoderType getDecoderType() @@ -83,10 +90,14 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi return DecoderType.P25_PHASE1; } + /** + * Constructs an instance + */ public P25P1DecoderLSM() { mMessageProcessor.setMessageListener(getMessageListener()); - mDemodulator = new P25P1DemodulatorLSM(mMessageFramer, this); + mMessageFramer.setPhaseWrappedSoftSymbols(true); + mSymbolProcessor = new P25P1SymbolProcessorLSM(this, mMessageFramer); } @Override @@ -111,35 +122,27 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi //Identify decimation that gets as close to 4.0 Samples Per Symbol as possible (4800 x 4.0 = 19.2 kHz) while((sampleRate / decimation) >= 38400) -// while((sampleRate / decimation) > 50000) { decimation *= 2; } - try - { - mDecimationFilterI = DecimationFilterFactory.getRealDecimationFilter(decimation); - mDecimationFilterQ = DecimationFilterFactory.getRealDecimationFilter(decimation); - } - catch(Exception e) - { - LOGGER.error("Error getting decimation filter for sample rate [" + sampleRate + "] decimation [" + decimation + "]"); - } + mDecimationFilterI = DecimationFilterFactory.getRealDecimationFilter(decimation); + mDecimationFilterQ = DecimationFilterFactory.getRealDecimationFilter(decimation); float decimatedSampleRate = (float)sampleRate / decimation; //Set the decimated sample rate to use for PLL error reporting. setDecimatedSampleRate(decimatedSampleRate); - int symbolLength = 16; - float rolloff = 0.2f; + float[] basebandTaps = FilterFactory.getRemezLowPass(decimatedSampleRate, BASEBAND_FILTER_PASS_BAND, + BASEBAND_FILTER_STOP_BAND); + mBasebandFilterI = FilterFactory.getRealFilter(basebandTaps); + mBasebandFilterQ = FilterFactory.getRealFilter(basebandTaps); + float samplesPerSymbol = decimatedSampleRate / (float)SYMBOL_RATE; + + mDemodulator.setSamplesPerSymbol(samplesPerSymbol); + mSymbolProcessor.setSamplesPerSymbol(samplesPerSymbol); - float[] taps = FilterFactory.getRootRaisedCosine(decimatedSampleRate / 4800.0, symbolLength, rolloff); - mPulseShapingFilterI = new RealFIRFilter(taps); - mPulseShapingFilterQ = new RealFIRFilter(taps); - mBasebandFilterI = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); - mBasebandFilterQ = FilterFactory.getRealFilter(getBasebandFilter(decimatedSampleRate)); - mDemodulator.setSamplesPerSymbol(decimatedSampleRate / (float)SYMBOL_RATE); mMessageFramer.setListener(mMessageProcessor); } @@ -162,49 +165,25 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi i = mBasebandFilterI.filter(i); q = mBasebandFilterQ.filter(q); - i = mPulseShapingFilterI.filter(i); - q = mPulseShapingFilterQ.filter(q); + //Differential demodulation: phase change across each symbol period (radians per sample) + float[] demodulated = mDemodulator.demodulate(i, q, mSymbolProcessor.isSynchronized()); - //Demodulate samples into symbols with timing, sync detection, and message framing. - mDemodulator.process(i, q); + //Process demodulated samples into symbols and apply message sync detection and framing. The soft symbol + //reliabilities weight the soft decision decoding. + mSymbolProcessor.process(demodulated, mDemodulator.getReliabilities()); } /** - * Constructs a baseband filter for this decoder using the current sample rate + * Includes the coarse carrier offset that the differential demodulator removes ahead of the symbol processor, + * since the symbol processor's PLL error only represents the residual offset. The symbol processor invokes this + * method when the message framer decodes a valid NID, which also verifies the demodulator's carrier offset range. + * @param pllError residual carrier offset (radians per sample) from the symbol processor */ - private float[] getBasebandFilter(double sampleRate) + @Override + public void processPLLError(float pllError) { - if(BASEBAND_FILTERS.containsKey(sampleRate)) - { - return BASEBAND_FILTERS.get(sampleRate); - } - - FIRFilterSpecification specification = FIRFilterSpecification - .lowPassBuilder() - .sampleRate(sampleRate) - .passBandCutoff(7250) - .passBandAmplitude(1.0).passBandRipple(0.01) //.01 - .stopBandAmplitude(0.0).stopBandStart(8000) //6500 - .stopBandRipple(0.01).build(); - - float[] coefficients = null; - - try - { - coefficients = FilterFactory.getTaps(specification); - BASEBAND_FILTERS.put(sampleRate, coefficients); - } - catch(Exception fde) //FilterDesignException - { - System.out.println("Error"); - } - - if(coefficients == null) - { - throw new IllegalStateException("Unable to design low pass filter for sample rate [" + sampleRate + "]"); - } - - return coefficients; + mDemodulator.offsetVerified(); + super.processPLLError(pllError - mDemodulator.getFrequencyOffset()); } /** @@ -213,7 +192,7 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi @Override public void setBufferListener(Listener listener) { - mDemodulator.setBufferListener(listener); + mSymbolProcessor.setBufferListener(listener); } /** @@ -222,7 +201,7 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi @Override public void removeBufferListener(Listener listener) { - mDemodulator.setBufferListener(null); + setBufferListener(null); } /** @@ -231,7 +210,7 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi @Override public boolean hasBufferListeners() { - return mDemodulator.hasBufferListener(); + return mSymbolProcessor.hasBufferListener(); } @Override @@ -252,6 +231,7 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi { super.stop(); mMessageFramer.stop(); + mSymbolProcessor.dispose(); } /** @@ -262,8 +242,13 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi switch(sourceEvent.getEvent()) { case NOTIFICATION_FREQUENCY_CHANGE: + mSymbolProcessor.resetBalance(); + mDemodulator.reset(); + break; case NOTIFICATION_FREQUENCY_CORRECTION_CHANGE: - mDemodulator.resetPLL(); + //Tuner PPM changes are small (auto-PPM steps are <= 10 Hz) and the carrier offset estimator and symbol + //processor balance track the change. Resetting would discard the carrier offset estimate and force an + //offset ambiguity search when the offset exceeds the +/- 600 Hz fourth power estimator range. break; case NOTIFICATION_SAMPLE_RATE_CHANGE: setSampleRate(sourceEvent.getValue().doubleValue()); @@ -282,17 +267,17 @@ public class P25P1DecoderLSM extends FeedbackDecoder implements IByteBufferProvi LOGGER.info("Starting ..."); // String directory = "D:\\DQPSK Equalizer Research - P25\\"; //Windows - String directory = "/media/denny/T9/DQPSK Equalizer Research - P25/"; //Linux + String directory = "/run/media/denny/T9/DQPSK Equalizer Research - P25/"; //Linux // String file = directory + "P25-S2-LSM-20241225_040119_460500000_CNYICC_Onondaga_Onondaga_CC_0_baseband.wav"; // String file = directory + "P25-S4-LSM-TCH-Data-20250105_141051_453587500_CNYICC_Onondaga_T-Onondaga_CC_38_baseband.wav"; // String file = directory + "P25-S5-LSM-CCH-20250923_034304_460500000_CNYICC_Onondaga_Onondaga-CC_0_baseband.wav"; - String file = directory + "P25-S6-LSM-CCH-20251025_071546_460500000_CNYICC_Onondaga_Onondaga-CC_0_baseband.wav"; - - -// String file = directory + "P25-S3-C4FM-20241225_040459_152517500_NYSEG_Onondaga_Control_30_baseband.wav"; - - +// String file = directory + "P25-S6-LSM-CCH-20251025_071546_460500000_CNYICC_Onondaga_Onondaga-CC_0_baseband.wav"; +// String file = directory + "P25-S7-LSM-CCH-20261027_030041_460475000_Central-New-York-Interoperable-Communications-Consortium-(CNYICC)_Madison-Simulcast_Control_51_baseband.wav"; +// String file = directory + "P25-S8-LSM-TCH-20261027_030138_460250000_Central-New-York-Interoperable-Communications-Consortium-(CNYICC)_Madison-Simulcast_T-Control_52_baseband.wav"; +// String file = directory + "P25-S9-LSM-TCH-20261027_030248_460175000_Central-New-York-Interoperable-Communications-Consortium-(CNYICC)_Madison-Simulcast_T-Control_53_baseband.wav"; +// String file = directory + "P25-S10-LSM-TCH-20261027_031432_460175000_Central-New-York-Interoperable-Communications-Consortium-(CNYICC)_Madison-Simulcast_T-Control_11_baseband.wav"; + String file = directory + "P25-S11-LSM-TCH-20261027_031603_460250000_Central-New-York-Interoperable-Communications-Consortium-(CNYICC)_Madison-Simulcast_T-Control_12_baseband.wav"; boolean autoReplay = false; P25P1DecoderLSM decoder = new P25P1DecoderLSM(); diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorC4FM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorC4FM.java deleted file mode 100644 index 0ed6ebba..00000000 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorC4FM.java +++ /dev/null @@ -1,1278 +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 - * **************************************************************************** - */ - -package io.github.dsheirer.module.decode.p25.phase1; - -import io.github.dsheirer.bits.CorrectedBinaryMessage; -import io.github.dsheirer.bits.IntField; -import io.github.dsheirer.dsp.filter.interpolator.LinearInterpolator; -import io.github.dsheirer.dsp.symbol.Dibit; -import io.github.dsheirer.dsp.symbol.DibitDelayLine; -import io.github.dsheirer.dsp.symbol.DibitToByteBufferAssembler; -import io.github.dsheirer.edac.bch.BCH_63_16_23_P25; -import io.github.dsheirer.gui.viewer.sync.FMSyncResultsViewer; -import io.github.dsheirer.gui.viewer.sync.ISyncResultsListener; -import io.github.dsheirer.module.decode.FeedbackDecoder; -import io.github.dsheirer.module.decode.p25.phase1.sync.P25P1SoftSyncDetector; -import io.github.dsheirer.module.decode.p25.phase1.sync.P25P1SoftSyncDetectorFactory; -import io.github.dsheirer.sample.Listener; -import java.nio.ByteBuffer; -import java.util.Arrays; -import java.util.concurrent.CountDownLatch; -import org.apache.commons.math3.stat.descriptive.moment.StandardDeviation; - -/** - * P25 Phase 1 decoder for C4FM/4FSK signals. Processes PI/4 differential QPSK demodulated samples with soft sync - * detection. Optimizes symbol timing and PLL at each sync detection. - * - * Note: this demodulator applies a gain of ~1.2x to the decoded phase of each symbol to correct for an imbalance in the - * decoded constellation and adjusts this gain at each sync detection based on how the decoded sync symbols align with - * the ideal symbols in the sync pattern. This compressed symbol constellation/gain imbalance is observed in both P25 - * and DMR 4FSK modulated signals when using PI/4 DQPSK demodulation and is possibly introduced by the pulse shaping - * transmit filter. This imbalance is not observed in LSM signals which likely use direct phase manipulation. - */ -public class P25P1DemodulatorC4FM -{ - private static final float EQUALIZER_LOOP_GAIN = 0.15f; - private static final float EQUALIZER_MAXIMUM_PLL = (float)(Math.PI / 3.0); //+/- 800 Hz - private static final float EQUALIZER_MAXIMUM_GAIN = 1.25f; - private static final float EQUALIZER_RECALIBRATE_THRESHOLD = (float)(Math.PI / 8.0); - private static final float SOFT_SYMBOL_QUADRANT_BOUNDARY = (float)(Math.PI / 2.0); - private static final float SYNC_THRESHOLD_DETECTION = 80; - private static final float SYNC_THRESHOLD_OPTIMIZED = 80; - private static final float SYNC_THRESHOLD_EQUALIZED = 110; - private static final float TWO_PI = (float)(Math.PI * 2.0); - private static final int BUFFER_WORKSPACE_LENGTH = 1024; - private static final int DIBIT_LENGTH_NID = 33; //32 dibits (64 bits) +1 status - private static final int DIBIT_LENGTH_SYNC = 24; - private static final Correction INVALID_SYNC_DETECTION = new Correction(0d, Double.MAX_VALUE, 0f, 0f, 0f, 0f, 0); - private static final IntField NAC_FIELD = IntField.length12(0); - private static final IntField DUID_FIELD = IntField.length4(12); - - private final BCH_63_16_23_P25 mBCHDecoder = new BCH_63_16_23_P25(); - private final DibitDelayLine mSymbolDelayLine = new DibitDelayLine(DIBIT_LENGTH_SYNC); - private final DibitToByteBufferAssembler mDibitAssembler = new DibitToByteBufferAssembler(300); - private final Equalizer mEqualizer = new Equalizer(); - private final FeedbackDecoder mFeedbackDecoder; - private final NACTracker mNACTracker = new NACTracker(); - private final P25P1MessageFramer mMessageFramer; - private final P25P1SoftSyncDetector mSyncDetector = P25P1SoftSyncDetectorFactory.getDetector(); - private final P25P1SoftSyncDetector mSyncDetectorLagging = P25P1SoftSyncDetectorFactory.getDetector(); - private boolean mFineSync = false; - private double mMaxFineSyncTimingAdjustment; - private double mNoiseStandardDeviationThreshold; - private double mSamplePoint; - private double mSamplePointAdjustment; - private double mSamplePointAdjustmentIncrement; - private double mSamplePointAdjustmentMax; - private double mSamplesPerSymbol; - private float[] mBuffer; - private float mLaggingSyncOffset; - private float mOptimizeFineIncrement; - private int mBufferPointer; - private int mBufferReloadThreshold; - private int mSymbolsSinceLastSync = 0; - - //Debug for waveform and symbol analysis - private ISyncResultsListener mSyncResultsViewer; - private final boolean mDebugVisualize = false; - private final int mDebugSymbolStart = 0; - private int mDebugSymbolCounter = 0; - private int mDebugSampleCounter = 0; - - /** - * Constructs an instance - * @param messageFramer to receive symbol decisions (dibits) and sync notifications. - */ - public P25P1DemodulatorC4FM(P25P1MessageFramer messageFramer, FeedbackDecoder feedbackDecoder) - { - mMessageFramer = messageFramer; - mFeedbackDecoder = feedbackDecoder; - } - - /** - * Resets the PLL - */ - public void resetPLL() - { - mEqualizer.reset(); - } - - /** - * Primary input method for receiving a stream of demodulated samples to process into symbols. - * @param samples to process - */ - public void process(float[] samples) - { - //Create shadow copy of heap variables on stack - boolean fineSync = mFineSync; - double samplePoint = mSamplePoint; - double samplesPerSymbol = mSamplesPerSymbol; - int bufferPointer = mBufferPointer; - int bufferReloadThreshold = mBufferReloadThreshold; - int symbolsSinceLastSync = mSymbolsSinceLastSync; - - int samplesPointer = 0; - float softSymbol, scorePrimary; - Dibit delayedSymbol; - Correction correctionCandidate; - - while(samplesPointer < samples.length) - { - //Note: buffer pointer can become greater than reload threshold during timing optimization at sync detect - if(bufferPointer >= bufferReloadThreshold) - { - //Do reload - int copyLength = Math.min(BUFFER_WORKSPACE_LENGTH, samples.length - samplesPointer); - System.arraycopy(mBuffer, copyLength, mBuffer, 0, mBuffer.length - copyLength); - System.arraycopy(samples, samplesPointer, mBuffer, mBuffer.length - copyLength, copyLength); - samplesPointer += copyLength; - bufferPointer -= copyLength; - - //Unwrap phases - for(int x = mBuffer.length - copyLength; x < mBuffer.length; x++) - { - if(mBuffer[x - 1] > 1.5f && mBuffer[x] < -1.5f) - { - mBuffer[x] += TWO_PI; - } - else if(mBuffer[x - 1] < -1.5f && mBuffer[x] > 1.5f) - { - mBuffer[x] -= TWO_PI; - } - } - } - - while(bufferPointer < bufferReloadThreshold) - { - bufferPointer++; - samplePoint--; - - if(samplePoint < 1) - { - symbolsSinceLastSync++; - softSymbol = mEqualizer.getEqualizedSymbol(mBuffer[bufferPointer], mBuffer[bufferPointer + 1], samplePoint); - - //Debug utility for viewing demodulation process -// getSyncResultsViewer().symbol(softSymbol); - - //Broadcast the decoded soft symbol to an optionally registered listener (ie Symbol view in channel panel). - mFeedbackDecoder.broadcast(softSymbol); - Dibit symbol = toSymbol(softSymbol); - samplePoint += mEqualizer.getAdjustment(softSymbol, symbol, bufferPointer); - - mMessageFramer.process(symbol); - - //We delay sending the symbol to the dibit assembler so that we have a chance to fully correct - //detected sync patterns in the delay buffer before they are sent downstream for recording. This - //ensures that replay of the recorded demodulated bit stream has fully corrected sync patterns - //when possible. - delayedSymbol = mSymbolDelayLine.insert(symbol); - mDibitAssembler.receive(delayedSymbol); - - scorePrimary = mSyncDetector.process(softSymbol); - correctionCandidate = INVALID_SYNC_DETECTION; - - if(fineSync) - { - if(symbolsSinceLastSync > 1 && scorePrimary > SYNC_THRESHOLD_DETECTION) - { - correctionCandidate = mEqualizer.optimize(0, scorePrimary, - bufferPointer + samplePoint); - } - } - else - { - //Check for lagging sync pattern - float lag = (float)(bufferPointer + samplePoint - mLaggingSyncOffset); - int lagIntegral = (int)Math.floor(lag); - float softSymbolLag = mEqualizer.getEqualizedSymbol(mBuffer[lagIntegral], mBuffer[lagIntegral + 1], lag - lagIntegral); - float scoreLag = mSyncDetectorLagging.process(softSymbolLag); - - if(symbolsSinceLastSync > 1) - { - if(scorePrimary > SYNC_THRESHOLD_DETECTION && scorePrimary > scoreLag) - { - correctionCandidate = mEqualizer.optimize(0.0f, scorePrimary, - bufferPointer + samplePoint); - } - - if(correctionCandidate == INVALID_SYNC_DETECTION && scoreLag > SYNC_THRESHOLD_DETECTION) - { - correctionCandidate = mEqualizer.optimize(-mLaggingSyncOffset, scoreLag, - bufferPointer + samplePoint); - } - } - } - - if(correctionCandidate.hasValidTiming()) - { - validateNID(correctionCandidate, bufferPointer + samplePoint); - - //Detect when the correction candidate fails to validate the NID and we have a high-quality - //original sync detection score. Attempt to validate the NID with no timing correction, using - //the original sync detection timing offset. - if(!correctionCandidate.isValid() && correctionCandidate.hasHighQualityDetectionScore()) - { - //Remove the candidate timing adjustment from the correction and use just the detection - //score as both detection and optimization. This limits the correction to just the equalizer - //update/adjustment. - correctionCandidate = mEqualizer.getCorrection(correctionCandidate.getAdditionalOffset(), - 0.0, correctionCandidate.getDetectionScore(), - correctionCandidate.getDetectionScore(), bufferPointer + samplePoint); - validateNID(correctionCandidate, bufferPointer + samplePoint); - } - - //Apply candidate correction values to timing and equalizer when we have a valid SYNC/NID - if(correctionCandidate.isValid()) - { - if(fineSync) - { - double adjustment = correctionCandidate.getTimingCorrection(); - - if(adjustment < 1) - { - adjustment = Math.max(adjustment, -mMaxFineSyncTimingAdjustment); - adjustment = Math.min(adjustment, mMaxFineSyncTimingAdjustment); - } - - samplePoint += adjustment; - } - else - { - samplePoint += correctionCandidate.getTimingCorrection(); - } - - //Reset the sample point intra-sync adjustment value. - mSamplePointAdjustment = 0; - - mEqualizer.apply(correctionCandidate); - - //broadcast equalizer balance as the current frequency offset - mFeedbackDecoder.processPLLError(mEqualizer.mPll); - } - - if(correctionCandidate.isValid() || correctionCandidate.hasHighQualityDetectionScore() || - (correctionCandidate.hasHighQualityOptimizationScore() && mMessageFramer.isComplete())) - { - mMessageFramer.syncDetected(); - fineSync = true; - symbolsSinceLastSync = 0; - - //Overwrite the sync symbols in the symbol delay line to eliminate any decoding errors - for(Dibit syncSymbol: mSyncDetector.getPatternDibits()) - { - mSymbolDelayLine.insert(syncSymbol); - } - } - } - - if(fineSync ^ mMessageFramer.isAssembling()) - { - fineSync = mMessageFramer.isAssembling(); - } - - //Add another symbol's worth of samples to the counter for processing - samplePoint += samplesPerSymbol; - } - } - } - - //Copy shadow variables back to the heap - mBufferPointer = bufferPointer; - mFineSync = fineSync; - mSamplePoint = samplePoint; - mSymbolsSinceLastSync = symbolsSinceLastSync; - } - - /** - * Registers the listener to receive demodulated bit stream buffers. - * @param listener to register - */ - public void setBufferListener(Listener listener) - { - mDibitAssembler.setBufferListener(listener); - } - - /** - * Indicates if there is a registered buffer listener - */ - public boolean hasBufferListener() - { - return mDibitAssembler.hasBufferListeners(); - } - - /** - * Sets or updates the samples per symbol - * @param samplesPerSymbol to apply. - */ - public void setSamplesPerSymbol(float samplesPerSymbol) - { - mSamplesPerSymbol = samplesPerSymbol; - mMaxFineSyncTimingAdjustment = samplesPerSymbol * .2; // 1/5th of a symbol period - mNoiseStandardDeviationThreshold = Dibit.D01_PLUS_3.getIdealPhase() * 2 / mSamplesPerSymbol * 1.2; //120% of optimal - mSamplePoint = samplesPerSymbol; - mSamplePointAdjustmentMax = samplesPerSymbol / 2; - mSamplePointAdjustmentIncrement = samplesPerSymbol / 100.0; - mOptimizeFineIncrement = samplesPerSymbol / 200f; - mLaggingSyncOffset = samplesPerSymbol / 2; - int bufferLength = BUFFER_WORKSPACE_LENGTH + (int)(Math.ceil((DIBIT_LENGTH_SYNC + DIBIT_LENGTH_NID + 2) * samplesPerSymbol)); - mBuffer = new float[bufferLength]; - mBufferReloadThreshold = mBuffer.length - (int)Math.ceil(samplesPerSymbol * (DIBIT_LENGTH_NID + 1)); - mBufferPointer = mBufferReloadThreshold; - } - - /** - * Resamples the NID dibits using the supplied adjustment. If the NID dibits pass error correction, notifies the - * message framer with the extracted NAC and DUID values and overwrites the NID dibit delay buffer with the - * static sync symbols and the resampled NID dibits so that the message framer and dibit assembler receive the - * optimally sampled sync and NID dibit sequences. - * @param correction (timing) to apply when resampling the NID. - */ - private void validateNID(Correction correction, double bufferOffset) - { - //Current sample point is pointing to the not yet corrected final sync symbol. Add in the candidate timing - // correction and resample the NID starting at the next symbol. - double pointer = bufferOffset + correction.getTimingCorrection() + mSamplesPerSymbol; - double fractional; - int integral; - float softSymbol = 0; - Dibit symbol; - - //Capture just the 63-bit BCH protected NID codeword including the 64th parity bit which we ignore. - CorrectedBinaryMessage candidateNID = new CorrectedBinaryMessage(64); - Dibit[] resampledNIDSymbols = new Dibit[33]; - - for(int x = 0; x < DIBIT_LENGTH_NID; x++) - { - integral = (int) Math.floor(pointer); - fractional = pointer - integral; - softSymbol = mEqualizer.getEqualizedSymbol(mBuffer[integral], mBuffer[integral + 1], fractional, correction); - symbol = toSymbol(softSymbol); - resampledNIDSymbols[x] = symbol; - - //Skip the status symbol at dibit 11 - if(x != 11) - { - candidateNID.add(symbol.getBit1(), symbol.getBit2()); - } - - pointer += mSamplesPerSymbol; - } - - int trackedNAC = mNACTracker.getTrackedNAC(); - mBCHDecoder.decode(candidateNID, trackedNAC); - - //If error correction fails, return the original correction candidate - if(candidateNID.getCorrectedBitCount() < 0) - { - correction.addCorrectedBitCount(candidateNID.getCorrectedBitCount()); - return; - } - - int nac = candidateNID.getInt(NAC_FIELD); - P25P1DataUnitID duid = P25P1DataUnitID.fromValue(candidateNID.getInt(DUID_FIELD)); - - //The BCH decoder can over-correct the NID and produce an invalid NAC. Compare it against the tracked NID to - //flag it as invalid NID when this happens. The NAC tracker will give us a value of 0 until it has enough - //observations of a valid NID value. - mNACTracker.track(nac); - - if(trackedNAC > 0 && trackedNAC != nac) - { - return; - } - - correction.setNID(nac, duid); - correction.addCorrectedBitCount(candidateNID.getCorrectedBitCount()); - } - - /** - * Decodes the sample value to determine the correct QPSK quadrant and maps the value to a Dibit symbol. - * @param sample in radians. - * @return symbol decision. - */ - public static Dibit toSymbol(float sample) - { - if(sample > 0) - { - return sample > SOFT_SYMBOL_QUADRANT_BOUNDARY ? Dibit.D01_PLUS_3 : Dibit.D00_PLUS_1; - } - else - { - return sample < -SOFT_SYMBOL_QUADRANT_BOUNDARY ? Dibit.D11_MINUS_3 : Dibit.D10_MINUS_1; - } - } - - /** - * Debug method to visualize the contents of the sync detection for the samples and symbols and symbol timing. - * @param score from the sync detector - * @param primary sync detector (true) or secondary (false) - */ -// public void visualizeSyncDetect(float score, boolean primary, String tag, int bufferPointer, double samplePoint) -// { -// int offset = 3; -// int length = (int)Math.ceil(mSamplesPerSymbol * 23) + (2 * offset); -// int end = bufferPointer + offset; -// int start = end - length; -// float[] symbols = new float[24]; -// float[] samples = Arrays.copyOfRange(mBuffer, start, end); -// for(int i = 0; i < samples.length; i++) -// { -// samples[i] = mEqualizer.equalize(samples[i]); -// } -// -// double symbolPointer = bufferPointer + samplePoint - (mSamplesPerSymbol * 23); -// int symbolIntegral = (int)Math.floor(symbolPointer); -// double mu = symbolPointer - symbolIntegral; -// -// for(int x = 0; x < 24; x++) -// { -// symbols[x] = mEqualizer.getEqualizedSymbol(mBuffer[symbolIntegral], mBuffer[symbolIntegral + 1], mu); -// symbolPointer += mSamplesPerSymbol; -// symbolIntegral = (int)Math.floor(symbolPointer); -// mu = symbolPointer - symbolIntegral; -// } -// -// float[] syncIntervals = new float[24]; -// int adjust = bufferPointer - length + offset; -// double pointer = bufferPointer + samplePoint - adjust; -// -// for(int x = 23; x >= 0; x--) -// { -// syncIntervals[x] = (float)pointer; -// pointer -= mSamplesPerSymbol; -// } -// -// CountDownLatch countDownLatch = new CountDownLatch(1); -// getSyncResultsViewer().receive(symbols, mSyncDetector.getSymbols(), samples, syncIntervals, mEqualizer.mPll, mEqualizer.mGain, -// "Score: " + score + (primary ? " PRIMARY " : " SECONDARY ") + (mFineSync ? "FINE " : "COARSE ") + -// " EQ-B:" + mEqualizer.mPll + " EQ-G:" + mEqualizer.mGain + " " + tag, countDownLatch); -// -// try -// { -// countDownLatch.await(); -// } -// catch(InterruptedException e) -// { -// throw new RuntimeException(e); -// } -// } - - /** - * Debug method to visualize the symbol decisions - */ -// public void visualizeSymbols(int bufferPointer, double samplePoint) -// { -// int offset = 3; -// int length = (int)Math.ceil(mSamplesPerSymbol * 23) + (2 * offset); -// int end = bufferPointer + offset; -// int start = end - length; -// float[] symbols = new float[24]; -// float[] decisions = new float[24]; -// float[] samples = Arrays.copyOfRange(mBuffer, start, end); -// for(int i = 0; i < samples.length; i++) -// { -// samples[i] = mEqualizer.equalize(samples[i]); -// } -// -// double symbolPointer = bufferPointer + samplePoint - (mSamplesPerSymbol * 23); -// int symbolIntegral = (int)Math.floor(symbolPointer); -// double mu = symbolPointer - symbolIntegral; -// -// for(int x = 0; x < 24; x++) -// { -// symbols[x] = mEqualizer.getEqualizedSymbol(mBuffer[symbolIntegral], mBuffer[symbolIntegral + 1], mu); -// decisions[x] = toSymbol(symbols[x]).getIdealPhase(); -// symbolPointer += mSamplesPerSymbol; -// symbolIntegral = (int)Math.floor(symbolPointer); -// mu = symbolPointer - symbolIntegral; -// } -// -// float[] syncIntervals = new float[24]; -// int adjust = bufferPointer - length + offset; -// double pointer = bufferPointer + samplePoint - adjust; -// -// for(int x = 23; x >= 0; x--) -// { -// syncIntervals[x] = (float)pointer; -// pointer -= mSamplesPerSymbol; -// } -// -// CountDownLatch countDownLatch = new CountDownLatch(1); -// getSyncResultsViewer().receive(symbols, decisions, samples, syncIntervals, mEqualizer.mPll, mEqualizer.mGain, -// "PLL:" + mEqualizer.mPll + " Gain:" + mEqualizer.mGain, countDownLatch); -// -// try -// { -// countDownLatch.await(); -// } -// catch(InterruptedException e) -// { -// throw new RuntimeException(e); -// } -// } - - /** - * Indicates if the samples in the sample buffer that contain a detected sync pattern have a standard deviation - * that is lower than the expected sample to sample deviation for a modulated signal. False detects against - * noise tend to have a standard deviation that is 2x the expected value. - * @param offset to the sample representing the final symbol in the detected sync pattern. - * @return true if the standard deviation is more than expected. - */ - public boolean isNoisy(double offset) - { - StandardDeviation standardDeviation = new StandardDeviation(); - int start = (int)Math.floor(offset - (23 * mSamplesPerSymbol)); - int end = (int)Math.ceil(offset); - end = Math.min(end, mBuffer.length - 1); - - for(int i = start; i < end; i++) - { - standardDeviation.increment(mBuffer[i] - mBuffer[i + 1]); - } - - return standardDeviation.getResult() > mNoiseStandardDeviationThreshold; - } - - /** - * Symbol timing and equalizer correction candidate settings - */ - public static class Correction - { - private final double mAdditionalOffset; - private final float mPllAdjustment; - private final float mGainAdjustment; - private final float mDetectionScore; - private final float mOptimizationScore; - private final double mTimingAdjustment; - private int mNAC; - private P25P1DataUnitID mDataUnitID = P25P1DataUnitID.PLACE_HOLDER; - private int mCorrectedBitCount = 0; - - /** - * Constructs an instance - * @param timingAdjustment correction for symbol sampling - * @param pllAdjustment correction for equalizer - * @param gainAdjustment correction for equalizer - * @param bitErrorCount in the detected sync pattern - */ - public Correction(double additionalOffset, double timingAdjustment, float pllAdjustment, float gainAdjustment, - float detectionScore, float optimizationScore, int bitErrorCount) - { - mAdditionalOffset = additionalOffset; - mTimingAdjustment = timingAdjustment; - mPllAdjustment = pllAdjustment; - mGainAdjustment = gainAdjustment; - mDetectionScore = detectionScore; - mOptimizationScore = optimizationScore; - mCorrectedBitCount = bitErrorCount; - } - - @Override - public String toString() - { - StringBuilder sb = new StringBuilder(); - - if(hasValidTiming()) - { - if(isValid()) - { - sb.append("\tCORRECTION - NAC [").append(mNAC); - sb.append("] DUID [").append(mDataUnitID.toString()); - sb.append("] Bit Errors [").append(mCorrectedBitCount).append("]"); - } - else - { - sb.append("CORRECTION - INVALID NID"); - } - - sb.append(" At Offset [").append(mAdditionalOffset); - sb.append("] Adjustments-Timing [").append(mTimingAdjustment); - sb.append("] PLL [").append(mPllAdjustment); - sb.append("] Gain [").append(mGainAdjustment); - sb.append("] Detect Score [").append(mDetectionScore); - sb.append("] Optimize Score [").append(mOptimizationScore); - sb.append("]"); - } - else - { - sb.append("Correction - INVALID TIMING"); - } - return sb.toString(); - } - - /** - * Candidate timing correction - * @return correction - */ - public double getTimingCorrection() - { - return mAdditionalOffset + mTimingAdjustment; - } - - /** - * Additional offset adjustment based on the location of the triggering sync detector. The primary sync - * detector has no additional offset, but a lagging sync detector has a lagging offset. - * @return addtional offset from lagging sync detector. - */ - public double getAdditionalOffset() - { - return mAdditionalOffset; - } - - /** - * Original sync detection score. - * @return score - */ - public float getDetectionScore() - { - return mDetectionScore; - } - - /** - * Timing optimization score. - * @return score - */ - public float getOptimizationScore() - { - return mOptimizationScore; - } - - /** - * Equalizer PLL correction value. - * @return correction - */ - public float getPllAdjustment() - { - return mPllAdjustment; - } - - /** - * Equalizer gain correction value - * @return correction. - */ - public float getGainAdjustment() - { - return mGainAdjustment; - } - - /** - * Sets the detected data unit ID. - * @param dataUnitID that was detected. - */ - public void setNID(int nac, P25P1DataUnitID dataUnitID) - { - mNAC = nac; - mDataUnitID = dataUnitID; - } - - /** - * Adds to the bit error count with the number of bit errors corrected in the NID by the BCH decoder - * @param correctedBitCount for the corrected NID - */ - public void addCorrectedBitCount(int correctedBitCount) - { - mCorrectedBitCount += correctedBitCount; - } - - /** - * Number of bits that were corrected across the SYNC and NID - * @return corrected bit count. - */ - public int getCorrectedBitCount() - { - return mCorrectedBitCount; - } - - /** - * Detected data unit ID - * @return data unit ID - */ - public P25P1DataUnitID getDataUnitID() - { - return mDataUnitID; - } - - public int getNAC() - { - return mNAC; - } - - /** - * Indicates if this correction has valid timing adjustment information. - * @return true if the timing correction is not the default (ie no timing) double max value. - */ - public boolean hasValidTiming() - { - return getTimingCorrection() != Double.MAX_VALUE; - } - - /** - * Indicates if this correction has a valid NID or the optimization (or detection) score is high quality. - * @return true if valid. - */ - public boolean isValid() - { - return mDataUnitID != P25P1DataUnitID.PLACE_HOLDER && getOptimizationScore() > SYNC_THRESHOLD_EQUALIZED; - } - - /** - * Indicates if the original sync detection score was high quality and is a higher score thant the - * optimization score. This indicates a failed optimization that can be reevaluated for NID detection using - * no timing adjustment. - * @return true if the detection score is high quality. - */ - public boolean hasHighQualityDetectionScore() - { - return getDetectionScore() > SYNC_THRESHOLD_EQUALIZED; - } - - /** - * Indicates if the optimization score is high quality, higher than the equalization threshold. - * @return true if high quality - */ - public boolean hasHighQualityOptimizationScore() - { - return getOptimizationScore() >= getDetectionScore() && getOptimizationScore() > SYNC_THRESHOLD_EQUALIZED; - } - } - - /** - * Equalizer and Phase-Locked Loop (PLL) - * - * Note: equalizer is initialized with a gain value of 1.219, as determined through significant testing against - * both P25 and DMR C4FM modulated samples that employ pulse shaping filters. Gain value is adjusted and - * optimized with each sync detection. - */ - class Equalizer - { - private boolean mInitialized = false; - private float mPll = 0.0f; - private float mGain = 1.219f; - - /** - * Reset after center frequency change or PPM change. - */ - public void reset() - { - mInitialized = false; - mPll = 0.0f; - } - - /** - * On sync detection, calculates the optimal timing adjustment that achieves the highest sync correlation score - * or returns an INVALID_SYNC_DETECTION sentinel value if the correlation score doesn't exceed the threshold. - * - * @param additionalOffset from current bufferPointer and samplePoint to adjust for which sync detector - * triggered the event. This can be zero offset for the primary sync detector or a lagging (ie half symbol) - * offset for the lagging sync detector. - * @return optimized timing adjustment or NO_OPTIMIZATION sentinel value. - */ - public Correction optimize(double additionalOffset, float detectionScore, double bufferOffset) - { - //Offset points to the final symbol/sample in the buffer for the detected sync pattern. - double offset = bufferOffset + additionalOffset; - - //Reject any sync detections where the sample:sample standard deviation exceeds the noise threshold. - if(isNoisy(offset)) - { - return INVALID_SYNC_DETECTION; - } - - //Find the optimal symbol timing - boolean fineSync = mFineSync; - double samplesPerSymbol = mSamplesPerSymbol; - double stepSize = samplesPerSymbol / (fineSync ? 16.0 : 8.0); //Start at 1/8th for coarse & 1/16th for fine - double stepSizeMin = mOptimizeFineIncrement; - double adjustment = 0.0; - - //In coarse sync mode, constrain max adjustment to half a symbol period so that a lagging sync detect doesn't - //preempt a primary sync detect prematurely and cause a significant timing correction. So, we constrain the - //max correction to a half symbol period for coarse sync mode. - double adjustmentMax = fineSync ? samplesPerSymbol : (samplesPerSymbol / 2.0); - double candidate = offset; - - float scoreCenter = score(candidate, samplesPerSymbol); - - candidate = offset - stepSize; - float scoreLeft = score(candidate, samplesPerSymbol); - - candidate = offset + stepSize; - float scoreRight = score(candidate, samplesPerSymbol); - - while(stepSize > stepSizeMin && Math.abs(adjustment) <= adjustmentMax) - { - if(scoreLeft > scoreRight && scoreLeft > scoreCenter) - { - adjustment -= stepSize; - scoreRight = scoreCenter; - scoreCenter = scoreLeft; - - candidate = offset + adjustment - stepSize; - scoreLeft = score(candidate, samplesPerSymbol); - } - else if(scoreRight > scoreLeft && scoreRight > scoreCenter) - { - adjustment += stepSize; - scoreLeft = scoreCenter; - scoreCenter = scoreRight; - - candidate = offset + adjustment + stepSize; - scoreRight = score(candidate, samplesPerSymbol); - } - else - { - stepSize *= 0.5f; - - if(stepSize > stepSizeMin) - { - candidate = offset + adjustment - stepSize; - scoreLeft = score(candidate, samplesPerSymbol); - - candidate = offset + adjustment + stepSize; - scoreRight = score(candidate, samplesPerSymbol); - } - } - } - - //If we didn't find an optimal correlation score above the threshold, return a false sync. - if(scoreCenter < SYNC_THRESHOLD_OPTIMIZED) - { - return INVALID_SYNC_DETECTION; - } - - return getCorrection(additionalOffset, adjustment, detectionScore, scoreCenter, bufferOffset); - } - - /** - * Equalize the symbol by applying PLL correction and multiplying with gain value. - * @param symbol to equalize. - * @return equalized symbol - */ - public float equalize(float symbol) - { - return (symbol + mPll) * mGain; - } - - /** - * Adjusts the symbol balance from the symbol error to address multi-path and fading between sync periods. - */ - public double getAdjustment(float softSymbol, Dibit symbol, int bufferPointer) - { - double adjustment; - - if(softSymbol < symbol.getIdealPhase()) - { - if(mBuffer[bufferPointer] > mBuffer[bufferPointer + 1]) - { - adjustment = -mSamplePointAdjustmentIncrement; - } - else - { - adjustment = mSamplePointAdjustmentIncrement; - } - } - else - { - if(mBuffer[bufferPointer] > mBuffer[bufferPointer + 1]) - { - adjustment = mSamplePointAdjustmentIncrement; - } - else - { - adjustment = -mSamplePointAdjustmentIncrement; - } - } - - if(Math.abs(mSamplePointAdjustment + adjustment) <= mSamplePointAdjustmentMax) - { - mSamplePointAdjustment += adjustment; - return adjustment; - } - - return 0; - } - - /** - * Equalizes both samples with optional phase wrap detection and returns the interpolated value at the mu offset. - * - * @param sample1 of the sample stream - * @param sample2 of the sample stream - * @param mu interpolation point between sample 1 and sample 2 - * @return equalized symbol - */ - public float getEqualizedSymbol(float sample1, float sample2, double mu) - { - sample1 = equalize(sample1); - sample2 = equalize(sample2); - return LinearInterpolator.calculate(sample1, sample2, mu); - } - - /** - * Equalizes both samples using the candidate correction values and returns the interpolated value at the - * mu offset. - * @param sample1 to equalize - * @param sample2 to equalize - * @param mu interpolation between samples 1 and 2 - * @param correction to apply when equalizing samples 1 and 2 - * @return interpolated value. - */ - public float getEqualizedSymbol(float sample1, float sample2, double mu, Correction correction) - { - sample1 = (sample1 + mPll + correction.getPllAdjustment()) * (mGain + correction.getGainAdjustment()); - sample2 = (sample2 + mPll + correction.getPllAdjustment()) * (mGain + correction.getGainAdjustment()); - - if(sample1 > Math.PI) - { - sample1 -= TWO_PI; - } - else if(sample1 < -Math.PI) - { - sample1 += TWO_PI; - } - - if(sample2 > Math.PI) - { - sample2 -= TWO_PI; - } - else if(sample2 < -Math.PI) - { - sample2 += TWO_PI; - } - - return LinearInterpolator.calculate(sample1, sample2, mu); - } - - /** - * Calculates the sync correlation score at the specified offset and samples per symbol interval using the - * current equalizer balance and gain values. - * - * @param offset to the final symbol in the soft symbol buffer - * @param samplesPerSymbol spacing to test for. - * @return correlation score. - */ - public float score(double offset, double samplesPerSymbol) - { - return score(offset, samplesPerSymbol, mPll, mGain); - } - - /** - * Calculates the sync correlation score at the specified offset and samples per symbol interval. - * @param offset to the final symbol in the soft symbol buffer - * @param samplesPerSymbol spacing to test for. - * @return correlation score. - */ - public float score(double offset, double samplesPerSymbol, float balance, float gain) - { - int maxPointer = mBuffer.length - 1; - float softSymbol; - - double pointer = offset - (samplesPerSymbol * 23.0); - int bufferPointer = (int)Math.floor(pointer); - double fractional = pointer - bufferPointer; - - float score = 0; - - for(int x = 0; x < 24; x++) - { - if(bufferPointer < maxPointer) - { - softSymbol = LinearInterpolator.calculate(mBuffer[bufferPointer], mBuffer[bufferPointer + 1], fractional); - softSymbol = (softSymbol + balance) * gain; - } - else - { - softSymbol = 0.0f; - } - - score += softSymbol * mSyncDetector.getPatternSymbols()[x]; - pointer += samplesPerSymbol; - bufferPointer = (int)Math.floor(pointer); - fractional = pointer - bufferPointer; - } - - return score; - } - - /** - * Indicates if this equalizer is initialized, meaning it has processed at least 1x sync detection to establish - * initial gain and balance values. - * - * @return true if initialized. - */ - public boolean isInitialized() - { - return mInitialized; - } - - /** - * Applies equalizer correction settings. - * - * @param correction settings to apply - */ - public void apply(Correction correction) - { - //Re-initialize the equalizer any time the balance correction value exceeds the threshold. - if(mInitialized && Math.abs(correction.getPllAdjustment()) > EQUALIZER_RECALIBRATE_THRESHOLD) - { - mInitialized = false; - } - - //Limit equalizer adjustments using a control loop at each sync detect after the initial equalizer setup or - //apply as new settings on initial sync detect or any time an excess balance correction is supplied. - if(mInitialized) - { - mPll += (correction.getPllAdjustment() * EQUALIZER_LOOP_GAIN); - mGain += (correction.getGainAdjustment() * EQUALIZER_LOOP_GAIN); - } - else - { - mPll += correction.getPllAdjustment(); - mGain += correction.getGainAdjustment(); - } - - //Constrain to +/- PI/2 (+/- 1200 Hertz of offset) - mPll = Math.min(mPll, EQUALIZER_MAXIMUM_PLL); - mPll = Math.max(mPll, -EQUALIZER_MAXIMUM_PLL); - - //Constrain gain between 1.0f and 1.25f - mGain = Math.min(mGain, EQUALIZER_MAXIMUM_GAIN); - mGain = Math.max(mGain, 1.0f); - mInitialized = true; - } - - /** - * Calculates a correction update for the equalizer balance and gain when the sync pattern is detected in the - * sample buffer, using the supplied timing correction argument. Resamples the sync symbols and compares each - * soft symbol to the ideal symbol phase to develop average error measurements for balance and gain. On initial - * sync detection, the equalizer settings are applied to the samples in the buffer allowing the symbols to be - * resampled during coarse sync acquisition. - */ - public Correction getCorrection(double additionalOffset, double timingCorrection, float detectionScore, - float optimizationScore, double offset) - { - double resampleStart = offset + additionalOffset + timingCorrection; - int resampleStartIntegral = (int)Math.floor(resampleStart); - float symbol = mSyncDetector.getPatternSymbols()[23]; - float resampledSoftSymbol = LinearInterpolator.calculate(mBuffer[resampleStartIntegral], - mBuffer[resampleStartIntegral + 1], resampleStart - resampleStartIntegral); - resampledSoftSymbol = (resampledSoftSymbol + mPll) * mGain; - - float balancePlus3Symbols = 0; - float balanceMinus3Symbols = resampledSoftSymbol - symbol; - float gainAccumulator = Math.abs(symbol) - Math.abs(resampledSoftSymbol); - Dibit resampledDibit = toSymbol(resampledSoftSymbol); - int bitErrorCount = mSyncDetector.getPatternDibits()[23].getBitErrorFrom(resampledDibit); - - resampleStart -= (23 * mSamplesPerSymbol); - resampleStartIntegral = (int)Math.floor(resampleStart); - - for(int x = 0; x < 23; x++) - { - if(resampleStartIntegral >= 0) - { - symbol = mSyncDetector.getPatternSymbols()[x]; - resampledSoftSymbol = LinearInterpolator.calculate(mBuffer[resampleStartIntegral], - mBuffer[resampleStartIntegral + 1], resampleStart - resampleStartIntegral); - resampledSoftSymbol = (resampledSoftSymbol + mPll) * mGain; - resampledDibit = toSymbol(resampledSoftSymbol); - bitErrorCount += mSyncDetector.getPatternDibits()[x].getBitErrorFrom(resampledDibit); - - if(symbol > 0) - { - balancePlus3Symbols += (resampledSoftSymbol - symbol); - } - else - { - balanceMinus3Symbols += (resampledSoftSymbol - symbol); - } - - gainAccumulator += Math.abs(symbol) - Math.abs(resampledSoftSymbol); - } - - resampleStart += mSamplesPerSymbol; - resampleStartIntegral = (int)Math.floor(resampleStart); - } - - balancePlus3Symbols /= -11.0f; //There are 11x Plus 3 and 13x Minus 3 symbols in the sync pattern. - balanceMinus3Symbols /= -13.0f; - float balanceCorrection = (balancePlus3Symbols + balanceMinus3Symbols) / 2f; - balanceCorrection = Math.min(balanceCorrection, SOFT_SYMBOL_QUADRANT_BOUNDARY); - balanceCorrection = Math.max(balanceCorrection, -SOFT_SYMBOL_QUADRANT_BOUNDARY); - - // System.out.println("Balance [" + balanceAverage + "] Plus3 [" + balancePlus3Symbols + "] Minus3 [" + balanceMinus3Symbols + "]"); - gainAccumulator /= (24.0f * Dibit.D01_PLUS_3.getIdealPhase()); - - //Recalculate the optimization score using the adjusted equalizer settings - if(!isInitialized()) - { - optimizationScore = score(offset + timingCorrection, mSamplesPerSymbol, - mPll + balanceCorrection, mGain + gainAccumulator); - } - - return new Correction(additionalOffset, timingCorrection, balanceCorrection, gainAccumulator, detectionScore, - optimizationScore, bitErrorCount); - } - } - - /** - * Debug viewer UI for visualizing sync detections in the sample stream/waveform. - */ - private ISyncResultsListener getSyncResultsViewer() - { - if(mSyncResultsViewer == null) - { - mSyncResultsViewer = new FMSyncResultsViewer(); - } - - return mSyncResultsViewer; - } - - /** - * Debug method to visualize the contents of the sync detection for the samples and symbols and symbol timing. - * @param score from the sync detector - * @param primary sync detector (true) or secondary (false) - */ - public void visualizeSyncDetect(float score, boolean primary, String tag, int bufferPointer, double samplePoint) - { - if(!mDebugVisualize) - { - return; - } - - int integral = (int)Math.floor(samplePoint); - samplePoint -= integral; - bufferPointer += integral; - - int offset = 3; - int length = (int)Math.ceil(mSamplesPerSymbol * 9) + (2 * offset); - int end = bufferPointer + offset; - int start = end - length; - float[] symbols = new float[10]; - float[] samples = Arrays.copyOfRange(mBuffer, start, end); - for(int i = 0; i < samples.length; i++) - { -// samples[i] = mSampleEqualizer.equalize(samples[i]); - } - - double symbolPointer = bufferPointer + samplePoint - (mSamplesPerSymbol * 9); - int symbolIntegral = (int)Math.floor(symbolPointer); - double mu = symbolPointer - symbolIntegral; - - for(int x = 0; x < 10; x++) - { -// symbols[x] = mSampleEqualizer.getEqualizedSymbol(mBuffer[symbolIntegral], mBuffer[symbolIntegral + 1], mu); - symbolPointer += mSamplesPerSymbol; - symbolIntegral = (int)Math.floor(symbolPointer); - mu = symbolPointer - symbolIntegral; - } - - float[] syncIntervals = new float[10]; - int adjust = bufferPointer - length + offset; - double pointer = bufferPointer + samplePoint - adjust; - - for(int x = 9; x >= 0; x--) - { - syncIntervals[x] = (float)pointer; - pointer -= mSamplesPerSymbol; - } - - CountDownLatch countDownLatch = new CountDownLatch(1); -// getSyncResultsViewer().receive(symbols, mSyncDetector.getSyncSymbols(), samples, syncIntervals, mSampleEqualizer.mBalance, mSampleEqualizer.mGain, -// tag + " Score: " + score + (primary ? " PRIMARY " : " SECONDARY ") + (mSynchronized ? "FINE " : "COARSE ") + -// " EQ-B:" + mSampleEqualizer.mBalance + " EQ-G:" + mSampleEqualizer.mGain + " " + -// " SAM:" + mDebugSampleCounter + " SYM:" + mDebugSymbolCounter + " DUR:" + DF.format(mDebugSymbolCounter / 2400.0), -// countDownLatch); - - try - { - countDownLatch.await(); - } - catch(InterruptedException e) - { - throw new RuntimeException(e); - } - } - - /** - * Debug method to visualize the symbol decisions - */ - public void visualizeSymbols(int bufferPointer, double samplePoint) - { - if(!mDebugVisualize) - { - return; - } - - if(mDebugSymbolCounter < mDebugSymbolStart) - { - return; - } - - int offset = 3; - int length = (int)Math.ceil(mSamplesPerSymbol * 9) + (2 * offset); - int end = bufferPointer + offset; - int start = end - length; - float[] symbols = new float[10]; - float[] decisions = new float[10]; - float[] samples = Arrays.copyOfRange(mBuffer, start, end); - for(int i = 0; i < samples.length; i++) - { -// samples[i] = mSampleEqualizer.equalize(samples[i]); - } - - double symbolPointer = bufferPointer + samplePoint - (mSamplesPerSymbol * 9); - int symbolIntegral = (int)Math.floor(symbolPointer); - double mu = symbolPointer - symbolIntegral; - - for(int x = 0; x < 10; x++) - { -// symbols[x] = mSampleEqualizer.getEqualizedSymbol(mBuffer[symbolIntegral], mBuffer[symbolIntegral + 1], mu); - decisions[x] = Dibit.fromSample(symbols[x]).getIdealPhase(); - symbolPointer += mSamplesPerSymbol; - symbolIntegral = (int)Math.floor(symbolPointer); - mu = symbolPointer - symbolIntegral; - } - - float[] syncIntervals = new float[10]; - int adjust = bufferPointer - length + offset; - double pointer = bufferPointer + samplePoint - adjust; - - for(int x = 9; x >= 0; x--) - { - syncIntervals[x] = (float)pointer; - pointer -= mSamplesPerSymbol; - } - - CountDownLatch countDownLatch = new CountDownLatch(1); -// getSyncResultsViewer().receive(symbols, decisions, samples, syncIntervals, mSampleEqualizer.mBalance, mSampleEqualizer.mGain, -// "Balance:" + mSampleEqualizer.mBalance + " Gain:" + mSampleEqualizer.mGain + -// " SAM:" + mDebugSampleCounter + " SYM:" + mDebugSymbolCounter + " DUR:" + DF.format(mDebugSymbolCounter / 2400.0), countDownLatch); - - try - { - countDownLatch.await(); - } - catch(InterruptedException e) - { - throw new RuntimeException(e); - } - } -} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorLSM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorLSM.java deleted file mode 100644 index 3b266ab6..00000000 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DemodulatorLSM.java +++ /dev/null @@ -1,391 +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 - * **************************************************************************** - */ - -package io.github.dsheirer.module.decode.p25.phase1; - -import io.github.dsheirer.dsp.filter.interpolator.LinearInterpolator; -import io.github.dsheirer.dsp.symbol.Dibit; -import io.github.dsheirer.dsp.symbol.DibitToByteBufferAssembler; -import io.github.dsheirer.gui.viewer.symbol.SymbolViewerFX; -import io.github.dsheirer.module.decode.FeedbackDecoder; -import io.github.dsheirer.sample.Listener; -import java.nio.ByteBuffer; -import java.util.Arrays; - -/** - * Demodulates filtered LSM I/Q samples to feed message framer for sync detection and framing. - */ -public class P25P1DemodulatorLSM -{ - private static final float HALF_PI = (float)(Math.PI / 2.0); - private static final float TWO_PI = (float)(Math.PI * 2.0); - private static final float MAX_PLL = (float)(Math.PI / 3.0); //+/- 800 Hz - private static final float OBJECTIVE_MAGNITUDE = 1.0f; - - private final DibitToByteBufferAssembler mDibitAssembler = new DibitToByteBufferAssembler(300); - private final FeedbackDecoder mFeedbackDecoder; - private final P25P1MessageFramer mMessageFramer; - private SymbolViewerFX mDebugSymbolViewer; - private double mSamplePoint; - private double mSamplesPerSymbol; - private double mSamplesPerHalfSymbol; - private float[] mBufferI; - private float[] mBufferQ; - private float mPLL = 0f; - private float mSampleGain = 1f; - private float mPreviousMiddleI, mPreviousMiddleQ, mPreviousCurrentI, mPreviousCurrentQ, mPreviousSymbolI = 0.7f, mPreviousSymbolQ = 0.7f; - private int mBufferPointer; - private int mBufferReserve; - - /** - * Constructs an instance - * @param messageFramer for receiving demodulated symbol stream and providing sync detection events. - * @param feedbackDecoder (parent) for receiving PLL carrier offset error reports. - */ - public P25P1DemodulatorLSM(P25P1MessageFramer messageFramer, FeedbackDecoder feedbackDecoder) - { - mMessageFramer = messageFramer; - mFeedbackDecoder = feedbackDecoder; - } - - /** - * Reset the PLL when the tuner source either changes center frequency or adjust the PPM value. - */ - public void resetPLL() - { - mPLL = 0f; - } - - /** - * Debug utility for viewing I/Q waveform, demodulated waveform and symbol decision details, sequentially. - */ - protected SymbolViewerFX getDebugSymbolViewer() - { - if(mDebugSymbolViewer == null) - { - mDebugSymbolViewer = new SymbolViewerFX(); - } - - return mDebugSymbolViewer; - } - - /** - * Primary input method for receiving a stream of filtered, pulse-shaped samples to process into symbols. - * @param i inphase samples to process - * @param q quadrature samples to process - */ - public void process(float[] i, float[] q) - { - //Shadow copy heap member variables onto the stack - double samplePoint = mSamplePoint; - double samplesPerSymbol = mSamplesPerSymbol; - double samplesPerHalfSymbol = mSamplesPerHalfSymbol; - float pll = mPLL; - float previousSymbolI = mPreviousSymbolI; - float previousSymbolQ = mPreviousSymbolQ; - float previousMiddleI = mPreviousMiddleI; - float previousMiddleQ = mPreviousMiddleQ; - float previousCurrentI = mPreviousCurrentI; - float previousCurrentQ = mPreviousCurrentQ; - float sampleGain = mSampleGain; - int bufferPointer = mBufferPointer; - - float pllGain = 0.1f; - double tedGain = samplesPerSymbol / 4; - double maxTimingAdjustment = samplesPerSymbol / 25; - double pointer, residual, timingAdjustment; - float magnitude, phaseError = 0, requiredGain, softSymbol, pllI, pllQ, pllTemp; - float iMiddle, qMiddle, iCurrent, qCurrent, iMiddleDemodulated, qMiddleDemodulated, iSymbol, qSymbol; - int offset; - Dibit hardSymbol; - - //I/Q buffers are the same length as incoming samples padded by an extra symbol length for processing space. - int bufferReserve = mBufferReserve; - int requiredLength = i.length + bufferReserve; - - if(mBufferI.length != requiredLength) - { - mBufferI = Arrays.copyOf(mBufferI, requiredLength); - mBufferQ = Arrays.copyOf(mBufferQ, requiredLength); - } - - //Transfer extra reserve samples from last processing iteration to the front of the I/Q buffers - System.arraycopy(mBufferI, bufferPointer, mBufferI, 0, bufferReserve); - System.arraycopy(mBufferQ, bufferPointer, mBufferQ, 0, bufferReserve); - - //Copy incoming I/Q samples to the processing buffers - System.arraycopy(i, 0, mBufferI, bufferReserve, i.length); - System.arraycopy(q, 0, mBufferQ, bufferReserve, q.length); - bufferPointer = 0; - int bufferReload = mBufferI.length - bufferReserve; - - while(bufferPointer < bufferReload) - { - bufferPointer++; - samplePoint--; - - if(samplePoint < 1) - { - //Sample point is the middle sample, between the previous and current symbols. - iMiddle = LinearInterpolator.calculate(mBufferI[bufferPointer], mBufferI[bufferPointer + 1], samplePoint); - qMiddle = LinearInterpolator.calculate(mBufferQ[bufferPointer], mBufferQ[bufferPointer + 1], samplePoint); - - //Calculate offset to next symbol. - pointer = bufferPointer + samplePoint + samplesPerHalfSymbol; - offset = (int)Math.floor(pointer); - residual = pointer - offset; - iCurrent = LinearInterpolator.calculate(mBufferI[offset], mBufferI[offset + 1], residual); - qCurrent = LinearInterpolator.calculate(mBufferQ[offset], mBufferQ[offset + 1], residual); - - //Adjust sample gain based on highest magnitude and apply to both middle and current samples. - magnitude = (float)(Math.sqrt(Math.pow(iCurrent, 2.0) + Math.pow(qCurrent, 2.0))); - - if(magnitude > 0 && !Float.isInfinite(magnitude)) - { - requiredGain = constrain(OBJECTIVE_MAGNITUDE / magnitude, 500); - sampleGain += (requiredGain - sampleGain) * .05f; - sampleGain = Math.min(sampleGain, requiredGain); - sampleGain = Math.min(sampleGain, 500); - } - - //Apply gain to the samples - iMiddle *= sampleGain; - qMiddle *= sampleGain; - iCurrent *= sampleGain; - qCurrent *= sampleGain; - - //Create I/Q representation of current PLL state - pllI = (float)Math.cos(pll); - pllQ = (float)Math.sin(pll); - - //Differential demodulation of middle sample - iMiddleDemodulated = (previousMiddleI * iMiddle) - (-previousMiddleQ * qMiddle); - qMiddleDemodulated = (previousMiddleI * qMiddle) + (-previousMiddleQ * iMiddle); - - //Rotate middle sample by the PLL offset - pllTemp = (iMiddleDemodulated * pllI) - (qMiddleDemodulated * pllQ); - qMiddleDemodulated = (qMiddleDemodulated * pllI) + (iMiddleDemodulated * pllQ); - iMiddleDemodulated = pllTemp; - - //Differential demodulation of symbol - iSymbol = (previousCurrentI * iCurrent) - (-previousCurrentQ * qCurrent); - qSymbol = (previousCurrentI * qCurrent) + (-previousCurrentQ * iCurrent); - - //Rotate symbol by the PLL offset - pllTemp = (iSymbol * pllI) - (qSymbol * pllQ); - qSymbol = (qSymbol * pllI) + (iSymbol * pllQ); - iSymbol = pllTemp; - - //Calculate the symbol (radians) from the I/Q values - softSymbol = (float)Math.atan2(qSymbol, iSymbol); - - //Apply Gardner timing error detection and correction - timingAdjustment = ((previousSymbolI - iSymbol) * iMiddleDemodulated) + ((previousSymbolQ - qSymbol) * qMiddleDemodulated); - timingAdjustment = constrain(timingAdjustment, maxTimingAdjustment); - timingAdjustment *= tedGain; - samplePoint += timingAdjustment; - - //Phase Locked Loop adjustment. Note: don't adjust PLL when soft symbol value is zero. - if(softSymbol != 0) - { - hardSymbol = toDibit(softSymbol); - phaseError = constrain(softSymbol - hardSymbol.getIdealPhase(), .3f); - pll -= (phaseError * pllGain); - pll = constrain(pll, MAX_PLL); - } - else - { - hardSymbol = Dibit.D00_PLUS_1; - } - - //Message framer returns boolean if valid sync and NID were detected/decoded. Send a PLL measured - //error report up to the tuner source so that it can self correct for mistuned channel and also to - //display the measured carrier offset value in the channel display - if(mMessageFramer.processWithSoftSyncDetect(softSymbol, hardSymbol)) - { - mFeedbackDecoder.processPLLError(pll); - } - - mDibitAssembler.receive(hardSymbol); - mFeedbackDecoder.broadcast(softSymbol); - -// if(bufferPointer >= 15) -// { -// SymbolViewerFX viewer = getDebugSymbolViewer(); -// CountDownLatch latch = new CountDownLatch(1); -// float[] rawI = Arrays.copyOfRange(mBufferI, bufferPointer - 15, bufferPointer + 15); -// float[] rawQ = Arrays.copyOfRange(mBufferQ, bufferPointer - 15, bufferPointer + 15); -// double middle = 15 + samplePoint;// - timingAdjustment; -// double previous = middle - samplesPerHalfSymbol; -// double current = middle + samplesPerHalfSymbol; -// double[] points = {previous, middle, current, //Timing -// previousCurrentI, iMiddle, iCurrent, //RawI -// previousCurrentQ, qMiddle, qCurrent, //RawQ -// previousSymbolI, iMiddleDemodulated, iSymbol, //DemodI -// previousSymbolQ, qMiddleDemodulated, qSymbol}; //DemodQ -// float middleSymbol = normalize((float)Math.atan2(qMiddleDemodulated, iMiddleDemodulated) + pll); -// float[] symbols = {previousSymbol, middleSymbol, softSymbol}; -// float adjustedPLL = pll + (phaseError * pllGain); -// viewer.receive(samplesPerSymbol, rawI, rawQ, sampleGain, adjustedPLL, points, symbols, latch); -// try -// { -// latch.await(); -// } -// catch(InterruptedException e) -// { -// throw new RuntimeException(e); -// } -// previousSymbol = softSymbol; -// } - - //Shuffle the values - previousSymbolI = iSymbol; - previousSymbolQ = qSymbol; - previousMiddleI = iMiddle; - previousMiddleQ = qMiddle; - previousCurrentI = iCurrent; - previousCurrentQ = qCurrent; - - //Add another symbol's worth of samples to the counter - samplePoint += samplesPerSymbol; - } - } - - //Copy shadow variables back to member variables. - mPLL = pll; - mBufferPointer = bufferPointer; - mPreviousMiddleI = previousMiddleI; - mPreviousMiddleQ = previousMiddleQ; - mPreviousCurrentI = previousCurrentI; - mPreviousCurrentQ = previousCurrentQ; - mPreviousSymbolI = previousSymbolI; - mPreviousSymbolQ = previousSymbolQ; - mSampleGain = sampleGain; - mSamplePoint = samplePoint; - } - - /** - * Normalizes the value in radians to the range of +/-(2 * PI) - * @param value to normalize - * @return normalized value. - */ - private static float normalize(float value) - { - if(value > Math.PI) - { - return value - TWO_PI; - } - else if(value < -Math.PI) - { - return value + TWO_PI; - } - - return value; - } - - /** - * Constrains value to range: -constraint to constrain - * @param value to constrain - * @param constraint as max positive value - * @return constrained value. - */ - private static double constrain(double value, double constraint) - { - if(Double.isNaN(value) || Double.isInfinite(value)) - { - return 0.0; - } - else if(value > constraint) - { - return constraint; - } - - return Math.max(value, -constraint); - } - - /** - * Constrains value to range: -constraint to constrain - * @param value to constrain - * @param constraint as max positive value - * @return constrained value. - */ - private static float constrain(float value, float constraint) - { - if(Float.isNaN(value) || Float.isInfinite(value)) - { - return 0.0f; - } - - if(value > constraint) - { - return constraint; - } - - return Math.max(value, -constraint); - } - - /** - * Sets or updates the samples per symbol - * @param samplesPerSymbol to apply. - */ - public void setSamplesPerSymbol(float samplesPerSymbol) - { - mSamplesPerSymbol = samplesPerSymbol; - mSamplesPerHalfSymbol = samplesPerSymbol / 2.0f; - mSamplePoint = samplesPerSymbol; - mBufferReserve = (int)Math.ceil(samplesPerSymbol); - mBufferI = new float[mBufferReserve]; - mBufferQ = new float[mBufferReserve]; - mBufferPointer = 0; - } - - /** - * Registers the listener to receive demodulated bit stream buffers. - * @param listener to register - */ - public void setBufferListener(Listener listener) - { - mDibitAssembler.setBufferListener(listener); - } - - /** - * Indicates if there is a registered buffer listener - */ - public boolean hasBufferListener() - { - return mDibitAssembler.hasBufferListeners(); - } - - /** - * Decodes the sample value to determine the correct QPSK quadrant and maps the value to a Dibit symbol. - * @param softSymbol in radians. - * @return symbol decision. - */ - public static Dibit toDibit(float softSymbol) - { - if(softSymbol > 0) - { - return softSymbol > HALF_PI ? Dibit.D01_PLUS_3 : Dibit.D00_PLUS_1; - } - else - { - return softSymbol < -HALF_PI ? Dibit.D11_MINUS_3 : Dibit.D10_MINUS_1; - } - } -} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSM.java new file mode 100644 index 00000000..5ac16507 --- /dev/null +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSM.java @@ -0,0 +1,371 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import io.github.dsheirer.dsp.filter.FilterFactory; +import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; +import io.github.dsheirer.dsp.filter.fir.real.RealFIRFilter; + +/** + * P25 Phase 1 LSM (CQPSK) differential demodulator that produces a real valued sample stream for the C4FM symbol + * processor. Each output sample is the phase change across the preceding symbol period, in radians per sample (ie + * the symbol phase divided by the samples per symbol), which is the same domain as the C4FM decoder's FM demodulated + * and receive filtered sample stream. + * + * Processing steps: + * + * 1. FM discriminator (atan2) produces the phase change between consecutive samples. A full range (atan2) + * discriminator is required because the LSM signal trajectory passes near the origin between symbols where the phase + * change between samples can exceed PI/2, which would fold with an atan(Q/I) discriminator. + * 2. Coarse frequency offset removal: subtracts the estimated carrier offset from each discriminator sample. + * 3. One symbol period integrator: the moving sum of the discriminator output is the phase change across the symbol + * period: arg(x[n] * conj(x[n - samplesPerSymbol])), modulo 2 * PI. + * 4. Phase wrap to +/- PI. The moving sum is only correct modulo 2 * PI since a discriminator sample near the origin + * can be wrapped by 2 * PI, and that error persists across the integrator span. + * 5. Carrier offset estimation: the fourth power of the differential QPSK symbols removes the modulation so that the + * residual phase is four times the carrier offset (radians per symbol). This is a constellation domain estimator, + * which isn't biased by the spectral tilt of a frequency selective (ie simulcast multipath) channel the way that a + * spectral centroid (ie mean discriminator output) estimator is. + * + * The fourth power estimator is ambiguous modulo PI/2 radians per symbol (1200 Hz at 4800 baud), so the estimator only + * resolves offsets within +/- 600 Hz. When a signal is present (ie the fourth power coherence exceeds the noise + * threshold) but the symbol processor is unable to synchronize, the estimate is stepped by +/- PI/2 radians per symbol + * to search the adjacent offset ranges. An offset just beyond the range aliases to an estimate near the opposite edge + * of the range (ie -650 Hz estimates as +550 Hz), so the search first steps toward the adjacent range nearest to the + * estimate (ie -1200 Hz for a positive estimate) and then the other adjacent range. The fourth power coherence is the + * same at each ambiguous offset, so high coherence without synchronization indicates an incorrect offset range, while + * low coherence indicates that there is no signal. The search holds the current offset range while there is no signal + * (ie between transmissions on a traffic channel) so that the next transmission starts with the offset range that last + * synchronized. Once an offset range is verified by a valid NID decode (see offsetVerified()), the search uses a + * longer dwell period before stepping away from it. Symbol processor synchronization alone doesn't verify the offset + * range since a false sync detection in an incorrect offset range would otherwise hold the search in that range for the + * longer dwell period. + * + * Soft symbol reliability: the demodulated phase change across the symbol period is the difference between the phases + * of the samples at either end of the period, and the phase noise variance of each sample is inversely proportional to + * its power, so the noise variance of the phase change is proportional to 1/|x[n]|^2 + 1/|x[n - samplesPerSymbol]|^2. + * The reliability (inverse noise variance) of each demodulated sample, (|x[n]|^2 * |x[n - samplesPerSymbol]|^2) / + * (|x[n]|^2 + |x[n - samplesPerSymbol]|^2), weights the soft decision decoding (see getReliabilities()). The LSM + * signal envelope passes near the origin between symbols and simulcast multipath produces deep envelope fades, so the + * symbol error rate varies by orders of magnitude with the envelope. Testing against a simulcast control channel + * recording found that 60% of the symbol errors occur in the 15% of the symbols with the lowest envelope. + * + * Since the carrier offset is removed ahead of the symbol processor, the symbol processor's balance (PLL) value only + * represents the residual offset. Use getFrequencyOffset() to include the coarse offset when reporting the carrier + * offset to the tuner. + */ +public class P25P1DifferentialDemodulatorLSM +{ + private static final double TWO_PI = 2.0 * Math.PI; + //Carrier offset ambiguity step of the fourth power estimator, in radians per symbol + private static final double AMBIGUITY_STEP = Math.PI / 2.0; + //The estimator, signal presence detection and ambiguity search update once per block of symbols, independent of + //the sample buffer size. The coherence (magnitude of the averaged unit vectors) of noise depends on the number of + //samples averaged, so a fixed block length also keeps the signal presence threshold independent of buffer size. + private static final int ESTIMATOR_BLOCK_SYMBOLS = 200; + //Coarse offset estimator loop gain applied to each block's fourth power phase error while the symbol processor is + //synchronized (~6,600 symbol time constant) + private static final double ESTIMATOR_GAIN = 0.03; + //Estimator loop gain while the symbol processor is unsynchronized, so that the estimate converges within the + //ambiguity search dwell period. With the tracking gain, the estimate converged over several dwell periods and an + //offset near the edge of the estimator range (ie -550 Hz) took longer to acquire than offsets beyond the range. + private static final double ESTIMATOR_ACQUISITION_GAIN = 0.2; + //Block coherence at and above which the estimator applies the full loop gain + private static final double ESTIMATOR_FULL_GAIN_COHERENCE = 0.3; + //Minimum average fourth power coherence (magnitude of the averaged unit vectors) to detect signal presence, to update + //the estimate and to advance the ambiguity search. LSM recordings produce values above 0.3 without injected noise + //and average 0.11-0.14 at 0 dB carrier to noise (50 kHz). Single block coherence values for noise reach 0.10, so + //the coherence is averaged across blocks, where the maximum noise value is 0.041 across 10 minutes of noise. + private static final double SIGNAL_MINIMUM_COHERENCE = 0.07; + //Coherence averaging gain, applied per block (~10 blocks or ~2,000 symbols) + private static final double COHERENCE_AVERAGING_GAIN = 0.1; + //Symbol count with a signal present but without synchronization before stepping to the next ambiguity search offset + private static final int SEARCH_DWELL_SYMBOLS = 2000; + //Symbol count with a signal present but without synchronization before stepping away from an offset range that + //was verified (~5 seconds). At low SNR, sync is frequently lost and reacquired, and the shorter search dwell caused + //the search to spend half of the time in an incorrect offset range at 0 dB carrier to noise. + private static final int VERIFIED_DWELL_SYMBOLS = 24000; + //Ambiguity search offsets, as multiples of the ambiguity step, relative to the estimator's offset. The search + //direction sets the sign of the adjacent range steps. + private static final int[] SEARCH_STEPS = {0, 1, -1}; + + private IRealFilter mIntegrator; + //Sample power history spanning the reliability symbol period delay, carried across sample buffers + private float[] mPowerHistory = new float[1]; + private float[] mReliabilities = new float[0]; + private double mSamplesPerSymbol = 1.0; + private float mPreviousI; + private float mPreviousQ; + //Fourth power estimator offset, radians per symbol + private double mEstimatorOffset = 0.0; + private int mSearchIndex = 0; + private int mSearchDirection = 1; + private double mUnsynchronizedSignalSymbols = 0; + //Indicates the current search offset range was verified by a valid NID decode since it was selected + private boolean mSearchIndexVerified = false; + private double mAverageCoherence = 0.0; + //Block length in samples and the fourth power accumulation for the block under accumulation + private int mBlockLength = ESTIMATOR_BLOCK_SYMBOLS; + private int mBlockPointer = 0; + private double mBlockI = 0.0; + private double mBlockQ = 0.0; + + /** + * Constructs an instance + */ + public P25P1DifferentialDemodulatorLSM() + { + } + + /** + * Sets or updates the samples per symbol. + * @param samplesPerSymbol for the incoming sample stream + */ + public void setSamplesPerSymbol(double samplesPerSymbol) + { + mSamplesPerSymbol = samplesPerSymbol; + mIntegrator = new RealFIRFilter(FilterFactory.getDifferentialIntegrator(samplesPerSymbol)); + mPowerHistory = new float[Math.max(1, (int)Math.round(samplesPerSymbol))]; + mBlockLength = (int)Math.round(ESTIMATOR_BLOCK_SYMBOLS * samplesPerSymbol); + resetBlock(); + } + + /** + * Resets the carrier offset estimate after a tuner center frequency change. + */ + public void reset() + { + mEstimatorOffset = 0.0; + mSearchIndex = 0; + mSearchDirection = 1; + mUnsynchronizedSignalSymbols = 0; + mSearchIndexVerified = false; + mAverageCoherence = 0.0; + resetBlock(); + } + + /** + * Resets the fourth power accumulation for the block under accumulation. + */ + private void resetBlock() + { + mBlockPointer = 0; + mBlockI = 0.0; + mBlockQ = 0.0; + } + + /** + * Notification that a valid NID was decoded, verifying the current carrier offset range so that the ambiguity search + * uses the longer dwell period before stepping away from it. + */ + public void offsetVerified() + { + mSearchIndexVerified = true; + } + + /** + * Carrier offset (coarse) removed by this demodulator. + * @return offset in radians per sample. + */ + public float getFrequencyOffset() + { + return (float)(getOffset() / mSamplesPerSymbol); + } + + /** + * Soft symbol reliabilities for each of the samples produced by the most recent demodulate() call. + * @return reliabilities (inverse of the relative phase noise variance), parallel to the demodulated samples. + */ + public float[] getReliabilities() + { + return mReliabilities; + } + + /** + * Current total offset including the ambiguity search step. + * @return offset in radians per symbol. + */ + private double getOffset() + { + return mEstimatorOffset + SEARCH_STEPS[mSearchIndex] * mSearchDirection * AMBIGUITY_STEP; + } + + /** + * Demodulates the filtered complex samples. + * @param i inphase samples + * @param q quadrature samples + * @param processorSynchronized indicates if the symbol processor is currently synchronized to the signal, which + * controls the ambiguity search. + * @return differentially demodulated samples in radians per sample. + */ + public float[] demodulate(float[] i, float[] q, boolean processorSynchronized) + { + float[] discriminated = new float[i.length]; + float previousI = mPreviousI; + float previousQ = mPreviousQ; + float demodulatedI, demodulatedQ; + + for(int x = 0; x < i.length; x++) + { + //Multiply the current sample by the complex conjugate of the previous sample + demodulatedI = (i[x] * previousI) + (q[x] * previousQ); + demodulatedQ = (q[x] * previousI) - (i[x] * previousQ); + discriminated[x] = (float)Math.atan2(demodulatedQ, demodulatedI); + previousI = i[x]; + previousQ = q[x]; + } + + mPreviousI = previousI; + mPreviousQ = previousQ; + mReliabilities = calculateReliabilities(i, q); + + float[] demodulated = mIntegrator.filter(discriminated); + + //Remove the carrier offset (the integrator has unity gain at DC, so the offset is removed after integration to + //apply each block's updated estimate), wrap each sample to +/- PI radians per symbol, and accumulate the fourth + //power of each symbol phase. + double samplesPerSymbol = mSamplesPerSymbol; + float wrap = (float)(TWO_PI / samplesPerSymbol); + float offset = getFrequencyOffset(); + double theta; + + for(int x = 0; x < demodulated.length; x++) + { + demodulated[x] -= offset; + demodulated[x] -= wrap * Math.round(demodulated[x] / wrap); + theta = 4.0 * demodulated[x] * samplesPerSymbol; + mBlockI += Math.cos(theta); + mBlockQ += Math.sin(theta); + + if(++mBlockPointer >= mBlockLength) + { + updateEstimate(processorSynchronized); + offset = getFrequencyOffset(); + } + } + + return demodulated; + } + + /** + * Updates the carrier offset estimate, the signal presence detection and the ambiguity search from the fourth power + * accumulation for the completed block. + * @param processorSynchronized state of the symbol processor + */ + private void updateEstimate(boolean processorSynchronized) + { + //Ideal symbols (+/- PI/4 and +/- 3PI/4) all map to PI in the fourth power, so the residual angle is four times + //the residual carrier offset (radians per symbol). + double coherence = Math.sqrt(mBlockI * mBlockI + mBlockQ * mBlockQ) / mBlockPointer; + mAverageCoherence += (coherence - mAverageCoherence) * COHERENCE_AVERAGING_GAIN; + boolean signalPresent = mAverageCoherence > SIGNAL_MINIMUM_COHERENCE; + + if(signalPresent) + { + //Weight the update by this block's coherence so that noise blocks (ie after a transmission ends, while the + //average coherence decays) don't move the estimate. + double error = Math.atan2(-mBlockQ, -mBlockI) / 4.0; + double weight = Math.min(1.0, coherence / ESTIMATOR_FULL_GAIN_COHERENCE); + double gain = processorSynchronized ? ESTIMATOR_GAIN : ESTIMATOR_ACQUISITION_GAIN; + mEstimatorOffset += error * gain * weight; + //Constrain the estimator to the fundamental range. The ambiguity search covers the adjacent ranges. + mEstimatorOffset = Math.clamp(mEstimatorOffset, -AMBIGUITY_STEP / 2, AMBIGUITY_STEP / 2); + } + + updateSearch(mBlockPointer, processorSynchronized, signalPresent); + resetBlock(); + } + + /** + * Calculates the reliability of each demodulated sample from the power of the samples at either end of the symbol + * period: (|x[n]|^2 * |x[n - samplesPerSymbol]|^2) / (|x[n]|^2 + |x[n - samplesPerSymbol]|^2) + * @param i inphase samples + * @param q quadrature samples + * @return reliabilities + */ + private float[] calculateReliabilities(float[] i, float[] q) + { + float[] history = mPowerHistory; + int delay = history.length; + float[] power = new float[i.length]; + float[] reliabilities = new float[i.length]; + float current, delayed, sum; + + for(int x = 0; x < i.length; x++) + { + power[x] = i[x] * i[x] + q[x] * q[x]; + } + + for(int x = 0; x < i.length; x++) + { + current = power[x]; + delayed = x >= delay ? power[x - delay] : history[x]; + sum = current + delayed; + reliabilities[x] = sum > 0.0f ? current * delayed / sum : 0.0f; + } + + //Retain the final (delay) sample powers for the next buffer, shifting forward when the buffer is short + if(i.length >= delay) + { + System.arraycopy(power, i.length - delay, history, 0, delay); + } + else + { + System.arraycopy(history, i.length, history, 0, delay - i.length); + System.arraycopy(power, 0, history, delay - i.length, i.length); + } + + return reliabilities; + } + + /** + * Steps the ambiguity search offset when a signal is present and the symbol processor hasn't synchronized within + * the dwell period. Holds the current search offset while there is no signal. + * @param sampleCount in the completed block + * @param processorSynchronized state of the symbol processor + * @param signalPresent indicates if a signal is present + */ + private void updateSearch(int sampleCount, boolean processorSynchronized, boolean signalPresent) + { + if(processorSynchronized) + { + mUnsynchronizedSignalSymbols = 0; + return; + } + + if(signalPresent) + { + mUnsynchronizedSignalSymbols += sampleCount / mSamplesPerSymbol; + + if(mUnsynchronizedSignalSymbols > (mSearchIndexVerified ? VERIFIED_DWELL_SYMBOLS : SEARCH_DWELL_SYMBOLS)) + { + mUnsynchronizedSignalSymbols = 0; + mSearchIndexVerified = false; + mSearchIndex = (mSearchIndex + 1) % SEARCH_STEPS.length; + + //Leaving the estimator's range: step first toward the adjacent range nearest to the estimate + if(mSearchIndex == 1) + { + mSearchDirection = mEstimatorOffset > 0.0 ? -1 : 1; + } + } + } + } +} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1FilterDesigner.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1FilterDesigner.java deleted file mode 100644 index 91ede962..00000000 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1FilterDesigner.java +++ /dev/null @@ -1,244 +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 - * **************************************************************************** - */ - -package io.github.dsheirer.module.decode.p25.phase1; - -/** - * Designs the P25 Phase 1 C4FM receive filter applied to the FM discriminator output. - * - * 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 (2880 Hz), 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 designer 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. - */ -public class P25P1FilterDesigner -{ - /** P25 C4FM raised cosine roll-off factor. */ - public static final double P25_ALPHA = 0.2; - - /** Least-squares weight applied in the stopband, relative to a weight of 1.0 in the passband. */ - private static final double STOPBAND_WEIGHT = 1.0; - - /** Frequency grid points per filter coefficient used for the least-squares fit. */ - private static final int GRID_DENSITY = 16; - - /** - * Designs the band limited integrate and dump receive filter. - * - * @param sampleRate in Hertz - * @param symbolRate in symbols per second - * @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[] designReceiveFilter(double sampleRate, double symbolRate, 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 + P25_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 + "]"); - } - - //Odd tap count so the filter has an integer group delay of (length - 1) / 2 samples - int length = (int)Math.ceil(symbolSpan * sampleRate / symbolRate); - - if(length % 2 == 0) - { - length++; - } - - int middle = (length - 1) / 2; - int coefficientCount = middle + 1; - double symbolPeriod = 1.0 / symbolRate; - - //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 from DC to Nyquist, skipping the transition band. - double[][] normal = new double[coefficientCount][coefficientCount]; - double[] rhs = new double[coefficientCount]; - double[] basis = new double[coefficientCount]; - int gridSize = GRID_DENSITY * coefficientCount; - - for(int g = 0; g <= gridSize; g++) - { - double frequency = nyquist * g / gridSize; - double desired; - double weight; - - if(frequency <= passBandEdge) - { - desired = sinc(frequency * symbolPeriod); - weight = 1.0; - } - else if(frequency >= stopBandStart) - { - desired = 0.0; - weight = STOPBAND_WEIGHT; - } - else - { - continue; //Transition band - } - - 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 and normalize to unity gain at DC - double[] taps = new double[length]; - double sum = 0.0; - - for(int m = 0; m < coefficientCount; m++) - { - taps[middle - m] = coefficients[m]; - taps[middle + m] = coefficients[m]; - } - - for(double tap : taps) - { - sum += tap; - } - - float[] output = new float[length]; - - for(int n = 0; n < 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; - } - - /** - * 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; - } -} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1Interleave.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1Interleave.java index 812ea518..b0c36dbf 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1Interleave.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1Interleave.java @@ -126,6 +126,7 @@ public class P25P1Interleave */ /** * Deinterleaves the soft symbols for a 196-bit (98 dibit) data block (TSBK, PDU header, unconfirmed data block). + * Also deinterleaves any other per-dibit values, such as the soft symbol reliabilities. * @param interleaved 98 soft symbols in transmitted order * @return 98 soft symbols in deinterleaved order. */ diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssembler.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssembler.java index 30acfb4e..b6871d31 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssembler.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssembler.java @@ -43,11 +43,19 @@ public class P25P1MessageAssembler implements Listener private static final Logger LOGGER = LoggerFactory.getLogger(P25P1MessageAssembler.class); private static final LoggingSuppressor LOGGING_SUPPRESSOR = new LoggingSuppressor(LOGGER); private static final int MINIMUM_LDU_BIT_LENGTH = 1500; //actual length is 1568 + //A place holder (failed NID) assembly also captures the sync (48 bits) and NID (64 bits) of the next frame, up to + //the point that the next frame's NID is decoded. + private static final int NEXT_FRAME_SYNC_NID_BIT_LENGTH = 112; + //Place holder length tolerance (bits) for identifying a missed TDU or TDULC frame. Testing against P25 recordings + //found missed TDU and TDULC place holders 1-12 bits shorter than the nominal length. + private static final int MISSED_FRAME_BIT_LENGTH_TOLERANCE = 12; private CorrectedBinaryMessage mMessage; private P25P1DataUnitID mDataUnitID; private int mNac; //Soft symbols (radians) for each assembled dibit, in received order, or NaN when not provided. private float[] mSoftSymbols = new float[256]; + //Soft symbol reliabilities for each assembled dibit, in received order, or NaN when not provided. + private float[] mReliabilities = new float[256]; private int mSoftSymbolCount = 0; /** @@ -96,11 +104,15 @@ public class P25P1MessageAssembler implements Listener } /** - * Forces this assembler to complete assembly of the current message. Attempts to identify + * Forces this assembler to complete assembly of the current message. Attempts to identify the message by either + * retaining the existing data unit ID when there was a valid NID decoded, or updating the place holder data unit + * ID based on the quantity of bits assembled thus far and the previous and next data unit IDs. * - * - * ion by either retaining the existing data unit ID when there was a valid NID - * decoded, or updating the data unit ID based on the quantity of dibits assembled thus far. + * A place holder is identified as a TDU only when its length matches a single missed TDU frame. A TDU has no + * error detection, so a TDU created from a place holder of any other length, ie inter-frame padding or a missed + * frame of another type, is always valid and falsely indicates the end of a call. A place holder that matches a + * single missed TDULC frame is identified as a TDULC, which is validated by its own error correction. A place + * holder that can't be identified remains a place holder, which is discarded. * @param previous data unit ID for the message prior to the currently assembled message * @param next data unit ID for the message that will be assembled next * @return number of bits that were dropped. @@ -113,11 +125,9 @@ public class P25P1MessageAssembler implements Listener { int length = getMessage().currentSize(); - if(length <= 28) - { - mDataUnitID = P25P1DataUnitID.TERMINATOR_DATA_UNIT; - } - else + //A place holder that's shorter than the shortest (TDU) missed frame is inter-frame padding or timing slack + //ahead of the next sync, not a missed frame, so it's left as a place holder and discarded. + if(length >= getMissedFrameMinimumLength(P25P1DataUnitID.TERMINATOR_DATA_UNIT)) { switch(next) { @@ -186,17 +196,55 @@ public class P25P1MessageAssembler implements Listener break; } } + + //If we're still PLACE_HOLDER, identify a single missed TDU or TDULC frame from the length + if(mDataUnitID == P25P1DataUnitID.PLACE_HOLDER) + { + if(isMissedFrame(length, P25P1DataUnitID.TERMINATOR_DATA_UNIT)) + { + mDataUnitID = P25P1DataUnitID.TERMINATOR_DATA_UNIT; + } + else if(isMissedFrame(length, P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL)) + { + mDataUnitID = P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL; + } + } } + //An unidentified place holder is discarded, so there are no dropped bits to report. The place holder bits + //are reported as sync loss. if(mDataUnitID == P25P1DataUnitID.PLACE_HOLDER) { - mDataUnitID = P25P1DataUnitID.TERMINATOR_DATA_UNIT; + return 0; } getMessage().setSize(mDataUnitID.getMessageLength()); return getMessage().size() - getMessage().currentSize(); } + /** + * Indicates if the place holder length matches a single missed frame of the data unit type: the message and null + * bits plus the sync and NID of the next frame. + * @param length of the place holder in bits + * @param dataUnitID of the missed frame + * @return true if the length is within tolerance of the missed frame length. + */ + private static boolean isMissedFrame(int length, P25P1DataUnitID dataUnitID) + { + return length >= getMissedFrameMinimumLength(dataUnitID) && + length <= dataUnitID.getMessageLength() + NEXT_FRAME_SYNC_NID_BIT_LENGTH; + } + + /** + * Minimum place holder length for a single missed frame of the data unit type. + * @param dataUnitID of the missed frame + * @return minimum length in bits. + */ + private static int getMissedFrameMinimumLength(P25P1DataUnitID dataUnitID) + { + return dataUnitID.getMessageLength() + NEXT_FRAME_SYNC_NID_BIT_LENGTH - MISSED_FRAME_BIT_LENGTH_TOLERANCE; + } + /** * The decoded NAC or zero if the NAC hasn't yet been fully decoded. * @return NAC @@ -231,6 +279,18 @@ public class P25P1MessageAssembler implements Listener * @param softSymbol the dibit was decided from (radians) or NaN if unavailable. */ public void receive(Dibit dibit, float softSymbol) + { + receive(dibit, softSymbol, Float.NaN); + } + + /** + * Primary input method for demodulated dibits/symbols, their soft symbol values and the soft symbol reliabilities + * to be appended to the message under assembly. Reliabilities weight the soft decision decoding. + * @param dibit to add to the message. + * @param softSymbol the dibit was decided from (radians) or NaN if unavailable. + * @param reliability of the soft symbol or NaN if unavailable. + */ + public void receive(Dibit dibit, float softSymbol, float reliability) { if(mMessage.isFull()) { @@ -244,8 +304,10 @@ public class P25P1MessageAssembler implements Listener if(mSoftSymbolCount >= mSoftSymbols.length) { mSoftSymbols = java.util.Arrays.copyOf(mSoftSymbols, mSoftSymbols.length * 2); + mReliabilities = java.util.Arrays.copyOf(mReliabilities, mReliabilities.length * 2); } + mReliabilities[mSoftSymbolCount] = reliability; mSoftSymbols[mSoftSymbolCount++] = softSymbol; } } @@ -279,6 +341,35 @@ public class P25P1MessageAssembler implements Listener return softSymbols; } + /** + * Soft symbol reliabilities for the assembled message bit range. + * @param startBit of the range (must be even/dibit aligned) + * @param endBit of the range, exclusive + * @return reliabilities for the range or null if the range isn't fully assembled or any reliability is unavailable. + */ + public float[] getReliabilities(int startBit, int endBit) + { + int startDibit = startBit / 2; + int endDibit = endBit / 2; + + if(startBit % 2 != 0 || endBit % 2 != 0 || startDibit < 0 || endDibit > mSoftSymbolCount || endDibit <= startDibit) + { + return null; + } + + float[] reliabilities = java.util.Arrays.copyOfRange(mReliabilities, startDibit, endDibit); + + for(float reliability : reliabilities) + { + if(Float.isNaN(reliability)) + { + return null; + } + } + + return reliabilities; + } + /** * Indicates if the message is fully assembled meaning that the quantity of assembled dibits matches the data unit * ID message length. diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageFramer.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageFramer.java index d101968b..35358081 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageFramer.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageFramer.java @@ -59,6 +59,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer //Soft symbols (radians) for the NID dibits, or NaN when the soft symbol wasn't provided private float[] mNIDSoftBuffer = new float[DIBIT_LENGTH_NID]; private float mCurrentSoftSymbol = Float.NaN; + private float mCurrentReliability = Float.NaN; //NID dibit index that carries the embedded status symbol private static final int NID_STATUS_DIBIT_INDEX = 11; //Soft decision NID decoding: maximum hard bit disagreements (of 63) allowed for the best candidate codeword. @@ -76,6 +77,8 @@ public class P25P1MessageFramer extends C4FMMessageFramer //Indicates the NID under capture follows a sync detected while a message was still assembling private boolean mPreemptiveNIDCheck = false; + //Indicates soft symbols are phase wrapped to +/- PI (ie P25 LSM) for soft decision decoding + private boolean mPhaseWrappedSoftSymbols = false; private Listener mMessageListener; private boolean mMessageAssemblyRequired = false; private boolean mRunning = false; @@ -99,6 +102,16 @@ public class P25P1MessageFramer extends C4FMMessageFramer { } + /** + * Sets the soft symbol domain for soft decision decoding. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie differentially demodulated P25 LSM) + * or false (default) for C4FM soft symbols that are not phase wrapped. + */ + public void setPhaseWrappedSoftSymbols(boolean phaseWrapped) + { + mPhaseWrappedSoftSymbols = phaseWrapped; + } + /** * Process soft symbol and apply soft symbol sync pattern detection. * @param softSymbol demodulated soft symbol @@ -162,6 +175,20 @@ public class P25P1MessageFramer extends C4FMMessageFramer return validNIDDetected; } + /** + * Process symbol decision, the soft symbol it was decided from and the soft symbol reliability, without sync + * detection. The reliability weights the soft decision decoding. + * @param symbol that was demodulated. + * @param softSymbol equalized soft symbol in radians. + * @param reliability of the soft symbol. + */ + @Override + public boolean process(Dibit symbol, float softSymbol, float reliability) + { + mCurrentReliability = reliability; + return process(symbol, softSymbol); + } + /** * Process hard symbol decision without sync detection. * @param symbol that was demodulated. @@ -171,6 +198,8 @@ public class P25P1MessageFramer extends C4FMMessageFramer boolean validNIDDetected = false; float softSymbol = mCurrentSoftSymbol; mCurrentSoftSymbol = Float.NaN; + float reliability = mCurrentReliability; + mCurrentReliability = Float.NaN; mDebugSymbolCount++; mDibitSinceTimestampCounter++; @@ -215,12 +244,12 @@ public class P25P1MessageFramer extends C4FMMessageFramer //If we still have an assembler, feed it the current dibit (e.g. TSBK and PDU continuation block assembly) if(mMessageAssembler != null) { - mMessageAssembler.receive(symbol, softSymbol); + mMessageAssembler.receive(symbol, softSymbol, reliability); } } else { - mMessageAssembler.receive(symbol, softSymbol); + mMessageAssembler.receive(symbol, softSymbol, reliability); } } //Start a message assembler after ignoring 24x Sync, 32x NID, and 1x status dibits. Trigger assembler @@ -406,7 +435,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer { case TRUNKING_SIGNALING_BLOCK_1: CorrectedBinaryMessage message1 = mMessageAssembler.getMessage().getSubMessage(0, 196); - TSBKMessage tsbk1 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message1, mMessageAssembler.getSoftSymbols(0, 196), mMessageAssembler.getNAC(), getTimestamp()); + TSBKMessage tsbk1 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message1, mMessageAssembler.getSoftSymbols(0, 196), mMessageAssembler.getReliabilities(0, 196), mPhaseWrappedSoftSymbols, mMessageAssembler.getNAC(), getTimestamp()); if(tsbk1 != null) { @@ -437,7 +466,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer break; case TRUNKING_SIGNALING_BLOCK_2: CorrectedBinaryMessage message2 = mMessageAssembler.getMessage().getSubMessage(196, 392); - TSBKMessage tsbk2 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message2, mMessageAssembler.getSoftSymbols(196, 392), mMessageAssembler.getNAC(), getTimestamp()); + TSBKMessage tsbk2 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message2, mMessageAssembler.getSoftSymbols(196, 392), mMessageAssembler.getReliabilities(196, 392), mPhaseWrappedSoftSymbols, mMessageAssembler.getNAC(), getTimestamp()); if(tsbk2 != null) { @@ -466,7 +495,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer break; case TRUNKING_SIGNALING_BLOCK_3: CorrectedBinaryMessage message3 = mMessageAssembler.getMessage().getSubMessage(392, 588); - TSBKMessage tsbk3 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message3, mMessageAssembler.getSoftSymbols(392, 588), mMessageAssembler.getNAC(), getTimestamp()); + TSBKMessage tsbk3 = TSBKMessageFactory.create(mChannelStatusProcessor.getDirection(), mMessageAssembler.getDataUnitID(), message3, mMessageAssembler.getSoftSymbols(392, 588), mMessageAssembler.getReliabilities(392, 588), mPhaseWrappedSoftSymbols, mMessageAssembler.getNAC(), getTimestamp()); broadcast(tsbk3); adjustDibitCounterFromMessageAssembler(); @@ -511,7 +540,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer { case PACKET_DATA_UNIT: CorrectedBinaryMessage message = mMessageAssembler.getMessage().getSubMessage(0, 196); - PDUHeader header = PDUMessageFactory.createHeader(message, mMessageAssembler.getSoftSymbols(0, 196)); + PDUHeader header = PDUMessageFactory.createHeader(message, mMessageAssembler.getSoftSymbols(0, 196), mMessageAssembler.getReliabilities(0, 196), mPhaseWrappedSoftSymbols); if(header != null) { @@ -546,7 +575,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer } else { - mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB1, mMessageAssembler.getSoftSymbols(196, 392))); + mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB1, mMessageAssembler.getSoftSymbols(196, 392), mMessageAssembler.getReliabilities(196, 392), mPhaseWrappedSoftSymbols)); } if(mPDUSequence.isComplete()) @@ -577,7 +606,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer } else { - mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB2, mMessageAssembler.getSoftSymbols(392, 588))); + mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB2, mMessageAssembler.getSoftSymbols(392, 588), mMessageAssembler.getReliabilities(392, 588), mPhaseWrappedSoftSymbols)); } if(mPDUSequence.isComplete()) @@ -608,7 +637,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer } else { - mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB3, mMessageAssembler.getSoftSymbols(588, 784))); + mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB3, mMessageAssembler.getSoftSymbols(588, 784), mMessageAssembler.getReliabilities(588, 784), mPhaseWrappedSoftSymbols)); } if(mPDUSequence.isComplete()) @@ -639,7 +668,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer } else { - mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB4, mMessageAssembler.getSoftSymbols(784, 980))); + mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB4, mMessageAssembler.getSoftSymbols(784, 980), mMessageAssembler.getReliabilities(784, 980), mPhaseWrappedSoftSymbols)); } if(mPDUSequence.isComplete()) @@ -670,7 +699,7 @@ public class P25P1MessageFramer extends C4FMMessageFramer } else { - mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB5, mMessageAssembler.getSoftSymbols(980, 1176))); + mPDUSequence.addDataBlock(PDUMessageFactory.createUnconfirmedDataBlock(messageB5, mMessageAssembler.getSoftSymbols(980, 1176), mMessageAssembler.getReliabilities(980, 1176), mPhaseWrappedSoftSymbols)); } adjustDibitCounterFromMessageAssembler(); diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageProcessor.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageProcessor.java index aab05ffd..9b4851d2 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageProcessor.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageProcessor.java @@ -119,9 +119,12 @@ public class P25P1MessageProcessor implements Listener //Reassemble extended source link control messages. if(message instanceof LDU1Message ldu1) { - LinkControlWord lcw = ldu1.getLinkControlWord();; + //The LDU1 message is valid for its voice frames, but the link control word that it carries is only used + //when it passes error correction. + LinkControlWord lcw = ldu1.getLinkControlWord(); + boolean lcwValid = lcw.isValid(); - if(lcw instanceof IExtendedSourceMessage esm) + if(lcwValid && lcw instanceof IExtendedSourceMessage esm) { if(lcw instanceof LCSourceIDExtension extension) { @@ -134,7 +137,7 @@ public class P25P1MessageProcessor implements Listener return; } } - else if(lcw instanceof LCHarrisTalkerAliasBase harrisTalkerAlias) + else if(lcwValid && lcw instanceof LCHarrisTalkerAliasBase harrisTalkerAlias) { //Send the LCW to the harris talker alias assembler dispatch(mHarrisTalkerAliasAssembler.process(harrisTalkerAlias, ldu1.getTimestamp())); @@ -146,12 +149,12 @@ public class P25P1MessageProcessor implements Listener dispatch(ldu1); //Process Harris Talker GPS messages - if(lcw instanceof LCHarrisTalkerGPSBlock1 block1) + if(lcwValid && lcw instanceof LCHarrisTalkerGPSBlock1 block1) { mHeldHarrisGPSBlock1 = block1; mHeldHarrisGPSBlock1Timestamp = ldu1.getTimestamp(); } - else if(lcw instanceof LCHarrisTalkerGPSBlock2 block2 && mHeldHarrisGPSBlock1 != null) + else if(lcwValid && lcw instanceof LCHarrisTalkerGPSBlock2 block2 && mHeldHarrisGPSBlock1 != null) { //Ensure block 1 and block 2 are within 2 seconds of each other if(Math.abs(ldu1.getTimestamp() - mHeldHarrisGPSBlock1Timestamp) < 2000) diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorC4FM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorC4FM.java index 88793125..1ccf0386 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorC4FM.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorC4FM.java @@ -2,7 +2,6 @@ package io.github.dsheirer.module.decode.p25.phase1; import io.github.dsheirer.dsp.c4fm.C4FMMessageFramer; import io.github.dsheirer.dsp.c4fm.C4FMSymbolProcessor; -import io.github.dsheirer.dsp.filter.interpolator.LinearInterpolator; import io.github.dsheirer.dsp.symbol.Dibit; import io.github.dsheirer.module.decode.FeedbackDecoder; import io.github.dsheirer.module.decode.p25.phase1.sync.P25P1SoftSyncDetectorFactory; @@ -205,7 +204,7 @@ public class P25P1SymbolProcessorC4FM extends C4FMSymbolProcessor } //If we didn't find an optimal correlation score above the threshold, return a false sync. - if(scoreCenter < SYNC_THRESHOLD_OPTIMIZED) + if(scoreCenter < getSyncThresholdOptimized()) { return INVALID_SYNC_DETECTION; } @@ -265,8 +264,8 @@ public class P25P1SymbolProcessorC4FM extends 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.gain() + mSampleEqualizer.balance(); resampledDibit = Dibit.fromSample(resampledSoftSymbol); bitErrorCount += patternDibits[x].getBitErrorFrom(resampledDibit); @@ -299,10 +298,11 @@ public class P25P1SymbolProcessorC4FM extends C4FMSymbolProcessor gainAccumulator = getGainAdjustment(mSyncDetector.getDibitLength() * Dibit.D01_PLUS_3.getIdealPhase(), gainAccumulator); - //Recalculate the optimization score using the adjusted equalizer settings + //Recalculate the optimization score using the adjusted equalizer settings, at the same (detector and timing + //adjusted) offset where the sync symbols were resampled. if(!mSynchronized) { - optimizationScore = score(offset + timingCorrection, mSamplesPerSymbol, + optimizationScore = score(offset + additionalOffset + timingCorrection, mSamplesPerSymbol, mSampleEqualizer.balance() + balanceCorrection, mSampleEqualizer.gain() + gainAccumulator); } diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSM.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSM.java new file mode 100644 index 00000000..8488195a --- /dev/null +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSM.java @@ -0,0 +1,131 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import io.github.dsheirer.dsp.c4fm.C4FMMessageFramer; +import io.github.dsheirer.module.decode.FeedbackDecoder; + +/** + * P25 Phase 1 LSM implementation of the symbol processor. Processes the differentially demodulated sample stream + * produced by the P25P1DifferentialDemodulatorLSM, which is in the same radians per sample domain as the C4FM decoder's + * sample stream, so this processor shares the P25 C4FM sync detection, optimization and tuning parameters. + * + * Note: testing against LSM recordings without injected noise showed no meaningful difference across equalizer spans + * of 1-11 symbols, with dynamic balance tracking disabled, or with the noise gate disabled. With injected noise, a + * single symbol equalizer span performed slightly better and dynamic balance tracking performed better. Disabling the + * continuous timing tracking significantly reduced valid messages on a simulcast recording. + * + * The demodulated sample stream is phase wrapped, so this processor interpolates across the phase wrap (see + * interpolate()). Measuring the sync correction balance error and the noise detection sample differences across the + * phase wrap made no measurable difference, so those use the C4FM implementations. + */ +public class P25P1SymbolProcessorLSM extends P25P1SymbolProcessorC4FM +{ + //Sync detection thresholds and the equalizer span were tuned with noise injection testing against LSM recordings at + //0-12 dB carrier to noise (50 kHz bandwidth). Relative to the C4FM values (80/100/100, 85% validation, 7 symbol + //equalizer), these increased valid messages 52% at 0 dB with 0.2% fewer valid messages without injected noise, and + //no false message detections against noise. The reference sync waveform validation threshold is the largest + //contributor since sync is frequently lost and reacquired at low SNR (see below). + private static final float SYNC_THRESHOLD_DETECTED = 70; + private static final float SYNC_THRESHOLD_OPTIMIZED = 85; + private static final float SYNC_THRESHOLD_EQUALIZED = 85; + //Reference sync waveform validation threshold for the first two sync detections after a sync loss. At low SNR, + //sync is frequently lost and reacquired, and true sync detections score 0.6 (median) at 6-9 dB carrier to noise + //(in-band), so a higher threshold rejects most reacquisitions: 0.75 rejected 17% of true syncs at 9 dB and 44% at + //6 dB. False sync detections are rejected by the optimization score threshold before validation. 0.45 increased + //valid messages 9% at 9 dB and 134% at 6 dB without false message detections against noise, while 0.3 produced a + //false NID decode against noise. + private static final double SYNC_WAVEFORM_VALIDATION_THRESHOLD = 0.45; + private static final int EQUALIZER_SYMBOL_LENGTH = 1; + //Phase wrap period of the demodulated sample stream: 2 * PI radians per symbol, in radians per sample + private float mWrapPeriod = (float)(2.0 * Math.PI); + + /** + * Constructs an instance + * + * @param feedbackDecoder parent decoder + * @param messageFramer to receive symbol decisions (dibits) and sync notifications. + */ + public P25P1SymbolProcessorLSM(FeedbackDecoder feedbackDecoder, C4FMMessageFramer messageFramer) + { + super(feedbackDecoder, messageFramer); + } + + @Override + public void setSamplesPerSymbol(float samplesPerSymbol) + { + mWrapPeriod = (float)(2.0 * Math.PI / samplesPerSymbol); + super.setSamplesPerSymbol(samplesPerSymbol); + } + + /** + * Interpolates across the phase wrap. The differentially demodulated sample stream is wrapped to +/- PI radians + * per symbol and the outer symbols (+/- 3PI/4) are only PI/4 from the wrap, so consecutive samples near an outer + * symbol can straddle the wrap (ie +3.0 and -3.0) where linear interpolation produces a value near zero. This + * unwraps the second sample relative to the first, interpolates, and wraps the result to +/- PI radians per symbol. + */ + @Override + protected float interpolate(float sample1, float sample2, double mu) + { + float wrap = mWrapPeriod; + float delta = sample2 - sample1; + delta -= wrap * Math.round(delta / wrap); + float interpolated = (float)(sample1 + delta * mu); + return interpolated - wrap * Math.round(interpolated / wrap); + } + + @Override + public int getEqualizerSymbolLength() + { + return EQUALIZER_SYMBOL_LENGTH; + } + + @Override + public float getSyncThresholdDetected() + { + return SYNC_THRESHOLD_DETECTED; + } + + @Override + public float getSyncThresholdOptimized() + { + return SYNC_THRESHOLD_OPTIMIZED; + } + + @Override + public float getSyncThresholdEqualized() + { + return SYNC_THRESHOLD_EQUALIZED; + } + + @Override + public double getSyncWaveformValidationThreshold() + { + return SYNC_WAVEFORM_VALIDATION_THRESHOLD; + } + + /** + * Indicates if the symbol processor is synchronized to the signal (ie sync detected within the max message length) + */ + public boolean isSynchronized() + { + return mSynchronized; + } +} diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/hdu/HDUMessage.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/hdu/HDUMessage.java index f0f477a9..cd4af4ab 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/hdu/HDUMessage.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/hdu/HDUMessage.java @@ -88,6 +88,17 @@ public class HDUMessage extends P25P1Message return P25P1DataUnitID.HEADER_DATA_UNIT; } + /** + * Indicates if the header data passed Reed-Solomon error detection and correction. Decodes the header data when + * it's not already decoded, since the decoding sets the message validity. + */ + @Override + public boolean isValid() + { + getHeaderData(); + return super.isValid(); + } + public HeaderData getHeaderData() { if(mHeaderData == null) diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/PDUMessageFactory.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/PDUMessageFactory.java index 0b7b0931..8e803b8e 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/PDUMessageFactory.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/PDUMessageFactory.java @@ -108,6 +108,20 @@ public class PDUMessageFactory * @return header or null. */ public static PDUHeader createHeader(CorrectedBinaryMessage message, float[] softSymbols) + { + return createHeader(message, softSymbols, null, false); + } + + /** + * Creates a PDU header message using soft decision Viterbi decoding when soft symbols are provided. + * @param message containing an interleaved, viterbi encoded message + * @param softSymbols for the 98 interleaved dibits of the header, or null. + * @param reliabilities for the 98 interleaved soft symbols to weight the soft decision decoding, or null. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie P25 LSM). + * @return header or null. + */ + public static PDUHeader createHeader(CorrectedBinaryMessage message, float[] softSymbols, float[] reliabilities, + boolean phaseWrapped) { if(softSymbols == null) { @@ -115,7 +129,7 @@ public class PDUMessageFactory } CorrectedBinaryMessage decoded = TSBKMessageFactory.deinterleaveViterbiAndCrc(message.getSubMessage(0, 196), - softSymbols); + softSymbols, reliabilities, phaseWrapped); if(decoded != null) { @@ -183,13 +197,29 @@ public class PDUMessageFactory * @param softSymbols for the 98 interleaved dibits of the block, or null. */ public static DataBlock createUnconfirmedDataBlock(CorrectedBinaryMessage interleaved, float[] softSymbols) + { + return createUnconfirmedDataBlock(interleaved, softSymbols, null, false); + } + + /** + * Creates an unconfirmed data block using soft decision Viterbi decoding when soft symbols are provided. + * @param interleaved data block + * @param softSymbols for the 98 interleaved dibits of the block, or null. + * @param reliabilities for the 98 interleaved soft symbols to weight the soft decision decoding, or null. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie P25 LSM). + */ + public static DataBlock createUnconfirmedDataBlock(CorrectedBinaryMessage interleaved, float[] softSymbols, + float[] reliabilities, boolean phaseWrapped) { if(softSymbols == null || softSymbols.length != P25P1Interleave.DATA_DEINTERLEAVE.length / 2) { return createUnconfirmedDataBlock(interleaved); } - return new UnconfirmedDataBlock(P25P1Interleave.deinterleaveDataSoftSymbols(softSymbols)); + float[] deinterleavedReliabilities = reliabilities != null && reliabilities.length == softSymbols.length ? + P25P1Interleave.deinterleaveDataSoftSymbols(reliabilities) : null; + return new UnconfirmedDataBlock(P25P1Interleave.deinterleaveDataSoftSymbols(softSymbols), + deinterleavedReliabilities, phaseWrapped); } public static DataBlock createUnconfirmedDataBlock(CorrectedBinaryMessage interleaved) diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/block/UnconfirmedDataBlock.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/block/UnconfirmedDataBlock.java index 25b45dfc..efaa081a 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/block/UnconfirmedDataBlock.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/pdu/block/UnconfirmedDataBlock.java @@ -47,7 +47,20 @@ public class UnconfirmedDataBlock extends DataBlock */ public UnconfirmedDataBlock(float[] deinterleavedSoftSymbols) { - mDecodedMessage = VITERBI_HALF_RATE_DECODER.decodeSoft(deinterleavedSoftSymbols); + this(deinterleavedSoftSymbols, null, false); + } + + /** + * Constructs an unconfirmed data block using soft decision Viterbi decoding. + * @param deinterleavedSoftSymbols 98 deinterleaved soft symbols for the 196-bit data block. + * @param deinterleavedReliabilities for each deinterleaved soft symbol to weight the soft decision decoding, or null. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie P25 LSM). + */ + public UnconfirmedDataBlock(float[] deinterleavedSoftSymbols, float[] deinterleavedReliabilities, + boolean phaseWrapped) + { + mDecodedMessage = VITERBI_HALF_RATE_DECODER.decodeSoft(deinterleavedSoftSymbols, deinterleavedReliabilities, + phaseWrapped); } /** diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tdu/TDULCMessage.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tdu/TDULCMessage.java index 90b0fbc4..024f682d 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tdu/TDULCMessage.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tdu/TDULCMessage.java @@ -74,6 +74,16 @@ public class TDULCMessage extends P25P1Message implements IFrequencyBandReceiver return P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL; } + /** + * Indicates if the link control word passed Reed-Solomon error detection and correction. Decodes the link control + * word when it's not already decoded. + */ + @Override + public boolean isValid() + { + return getLinkControlWord().isValid(); + } + public LinkControlWord getLinkControlWord() { if(mLinkControlWord == null) diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tsbk/TSBKMessageFactory.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tsbk/TSBKMessageFactory.java index 26941f72..82cd889d 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tsbk/TSBKMessageFactory.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase1/message/tsbk/TSBKMessageFactory.java @@ -157,13 +157,32 @@ public class TSBKMessageFactory * 2 bit errors is invalid. */ public static CorrectedBinaryMessage deinterleaveViterbiAndCrc(CorrectedBinaryMessage raw, float[] softSymbols) + { + return deinterleaveViterbiAndCrc(raw, softSymbols, null, false); + } + + /** + * Performs deinterleave, Viterbi decoding, and applies the CRC-CCITT error detection/correction. Uses soft + * decision Viterbi decoding when the soft symbols are provided, otherwise hard decision decoding. + * @param raw interleaved message + * @param softSymbols for the 98 interleaved dibits of the raw message, or null. + * @param reliabilities for the 98 interleaved soft symbols to weight the soft decision decoding, or null. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie P25 LSM). + * @return null or the decoded message with CRC results set in the message where 0 or 1 bit errors is valid and + * 2 bit errors is invalid. + */ + public static CorrectedBinaryMessage deinterleaveViterbiAndCrc(CorrectedBinaryMessage raw, float[] softSymbols, + float[] reliabilities, boolean phaseWrapped) { CorrectedBinaryMessage message; if(softSymbols != null && softSymbols.length == P25P1Interleave.DATA_DEINTERLEAVE.length / 2) { //Soft decision decode 1/2 rate trellis encoded block from the deinterleaved soft symbols - message = VITERBI_HALF_RATE_DECODER.decodeSoft(P25P1Interleave.deinterleaveDataSoftSymbols(softSymbols)); + float[] deinterleavedReliabilities = reliabilities != null && reliabilities.length == softSymbols.length ? + P25P1Interleave.deinterleaveDataSoftSymbols(reliabilities) : null; + message = VITERBI_HALF_RATE_DECODER.decodeSoft(P25P1Interleave.deinterleaveDataSoftSymbols(softSymbols), + deinterleavedReliabilities, phaseWrapped); } else { @@ -205,7 +224,21 @@ public class TSBKMessageFactory CorrectedBinaryMessage correctedBinaryMessage, float[] softSymbols, int nac, long timestamp) { - CorrectedBinaryMessage message = deinterleaveViterbiAndCrc(correctedBinaryMessage, softSymbols); + return create(direction, dataUnitID, correctedBinaryMessage, softSymbols, null, false, nac, timestamp); + } + + /** + * Creates a TSBK message, using soft decision Viterbi decoding when soft symbols are provided. + * @param softSymbols for the 98 interleaved dibits of the message, or null for hard decision decoding. + * @param reliabilities for the 98 interleaved soft symbols to weight the soft decision decoding, or null. + * @param phaseWrapped true when the soft symbols are phase wrapped to +/- PI (ie P25 LSM). + */ + public static TSBKMessage create(Direction direction, P25P1DataUnitID dataUnitID, + CorrectedBinaryMessage correctedBinaryMessage, float[] softSymbols, + float[] reliabilities, boolean phaseWrapped, int nac, long timestamp) + { + CorrectedBinaryMessage message = deinterleaveViterbiAndCrc(correctedBinaryMessage, softSymbols, reliabilities, + phaseWrapped); if(message == null) { diff --git a/src/main/java/io/github/dsheirer/module/decode/p25/phase2/P25P2DecoderHDQPSK.java b/src/main/java/io/github/dsheirer/module/decode/p25/phase2/P25P2DecoderHDQPSK.java index b1ca89a6..42b51983 100644 --- a/src/main/java/io/github/dsheirer/module/decode/p25/phase2/P25P2DecoderHDQPSK.java +++ b/src/main/java/io/github/dsheirer/module/decode/p25/phase2/P25P2DecoderHDQPSK.java @@ -23,8 +23,6 @@ import io.github.dsheirer.alias.AliasModel; import io.github.dsheirer.channel.state.MultiChannelState; import io.github.dsheirer.controller.channel.Channel; import io.github.dsheirer.dsp.filter.FilterFactory; -import io.github.dsheirer.dsp.filter.design.FilterDesignException; -import io.github.dsheirer.dsp.filter.fir.FIRFilterSpecification; import io.github.dsheirer.dsp.filter.fir.real.IRealFilter; import io.github.dsheirer.dsp.gain.complex.ComplexGainFactory; import io.github.dsheirer.dsp.gain.complex.IComplexGainControl; @@ -51,9 +49,7 @@ import io.github.dsheirer.source.wave.ComplexWaveSource; import java.io.IOException; import java.nio.file.Path; import java.nio.file.Paths; -import java.util.HashMap; import java.util.Iterator; -import java.util.Map; import org.slf4j.Logger; import org.slf4j.LoggerFactory; @@ -69,7 +65,6 @@ public class P25P2DecoderHDQPSK extends P25P2Decoder implements IdentifierUpdate protected CostasLoop mCostasLoop; protected P25P2MessageFramer mMessageFramer; protected IComplexGainControl mAGC = ComplexGainFactory.getComplexGainControl(); - private Map mBasebandFilters = new HashMap<>(); protected IRealFilter mIBasebandFilter; protected IRealFilter mQBasebandFilter; private DecodeConfigP25Phase2 mDecodeConfigP25Phase2; @@ -107,8 +102,10 @@ public class P25P2DecoderHDQPSK extends P25P2Decoder implements IdentifierUpdate { super.setSampleRate(sampleRate); - mIBasebandFilter = FilterFactory.getRealFilter(getBasebandFilter()); - mQBasebandFilter = FilterFactory.getRealFilter(getBasebandFilter()); + //Note: the baseband filter is designed for a 50 kHz sample rate, independent of the current sample rate. + float[] basebandTaps = FilterFactory.getRemezLowPass(50000.0, 6500, 0.005, 7200, 0.01); + mIBasebandFilter = FilterFactory.getRealFilter(basebandTaps); + mQBasebandFilter = FilterFactory.getRealFilter(basebandTaps); mCostasLoop = new CostasLoop(getSampleRate(), getSymbolRate()); mCostasLoop.setPLLBandwidth(PLLBandwidth.BW_300); @@ -152,48 +149,6 @@ public class P25P2DecoderHDQPSK extends P25P2Decoder implements IdentifierUpdate mQPSKDemodulator.receive(amplified); } - /** - * Constructs a baseband filter for this decoder using the current sample rate - */ - private float[] getBasebandFilter() - { - //Attempt to reuse a cached (ie already-designed) filter if available - float[] filter = mBasebandFilters.get(getSampleRate()); - - if(filter == null) - { - FIRFilterSpecification specification = FIRFilterSpecification.lowPassBuilder() - .sampleRate(50000.0) - .passBandCutoff(6500) - .passBandAmplitude(1.0) - .passBandRipple(0.005) - .stopBandAmplitude(0.0) - .stopBandStart(7200) - .stopBandRipple(0.01) - .build(); - - try - { - filter = FilterFactory.getTaps(specification); - } - catch(FilterDesignException fde) - { - mLog.error("Couldn't design low pass baseband filter for sample rate: " + getSampleRate()); - } - - if(filter != null) - { - mBasebandFilters.put(getSampleRate(), filter); - } - else - { - throw new IllegalStateException("Couldn't design a C4FM baseband filter for sample rate: " + getSampleRate()); - } - } - - return filter; - } - @Override protected void process(SourceEvent sourceEvent) { diff --git a/src/main/java/io/github/dsheirer/source/mixer/MixerReader.java b/src/main/java/io/github/dsheirer/source/mixer/MixerReader.java index ef298b0f..25e83865 100644 --- a/src/main/java/io/github/dsheirer/source/mixer/MixerReader.java +++ b/src/main/java/io/github/dsheirer/source/mixer/MixerReader.java @@ -260,7 +260,7 @@ public class MixerReader { while(bytesRead < buffer.length) { - bytesRead += mTargetDataLine.read(buffer, bytesRead, buffer.length); + bytesRead += mTargetDataLine.read(buffer, bytesRead, buffer.length - bytesRead); } } catch(ArrayIndexOutOfBoundsException aioobe) diff --git a/src/main/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManager.java b/src/main/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManager.java index aa780651..43eb38d6 100644 --- a/src/main/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManager.java +++ b/src/main/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManager.java @@ -42,6 +42,7 @@ public class ChannelFrequencyErrorManager implements ISourceEventListener, Liste private final TunerChannelSource mTunerChannelSource; private final TunerFrequencyErrorManager mParent; private ScheduledFuture mScheduledFuture; + //Guarded by this instance: updated by this manager's timer and by the tuner's PPM change notifications private long mAppliedFrequencyCorrection = 0; private long mAverageRequestedCorrection = 0; private boolean mCurrentForTunerProcessing = false; @@ -73,7 +74,7 @@ public class ChannelFrequencyErrorManager implements ISourceEventListener, Liste * Snapshot of the current applied frequency correction for the channel. * @return total frequency error in Hertz. */ - public long getFrequencyError() + public synchronized long getFrequencyError() { return mAppliedFrequencyCorrection; } @@ -81,12 +82,22 @@ public class ChannelFrequencyErrorManager implements ISourceEventListener, Liste /** * Receive a notification from the parent tuner is updating the PPM and frequency correction change value. * - * Note: on receiving this change notification, the tuner manager will expose the current PPM and correction values. + * The tuner PPM change shifts every channel by the change: a positive PPM change lowers the tuned frequency, which + * raises each channel's signal frequency by the change, the same as a positive channel frequency correction. The + * tuner absorbs the change from the average of the channel corrections, so remove the change from this channel's + * applied correction to keep the channel's net correction (and the decoder's view of the signal) unchanged. + * Otherwise, each PPM change appears to the decoder as a sudden carrier offset step that the channel correction + * removes over the following correction intervals, while the tuner continues to absorb the stale channel + * correction. + * + * Note: invoked by the tuner's frequency error manager while holding the tuner controller lock, so this method + * must not call back into the tuner. * @param change in Hertz resulting from the tuner's PPM update. */ - public void receive(long change) + public synchronized void receive(long change) { - //The TunerFrequencyErrorManager notifies this channel of the correction applied at the tuner + mAppliedFrequencyCorrection -= change; + mTunerChannelSource.setFrequencyCorrection(mAppliedFrequencyCorrection); } /** @@ -133,16 +144,25 @@ public class ChannelFrequencyErrorManager implements ISourceEventListener, Liste if(mCurrentForChannelProcessing) { long partialCorrection = Math.clamp(mAverageRequestedCorrection, -MAXIMUM_CHANNEL_ERROR_CORRECTION_PER_INTERVAL, MAXIMUM_CHANNEL_ERROR_CORRECTION_PER_INTERVAL); + long appliedFrequencyCorrection; - if(partialCorrection != 0) + //receive(long change) holds the tuner controller lock while it locks this instance, so the status + //notification below, which calls into the tuner, is created outside of this block to avoid a lock order + //inversion. + synchronized(this) { - mAppliedFrequencyCorrection += (long)(partialCorrection * ERROR_PROCESSING_GAIN); - mTunerChannelSource.setFrequencyCorrection(mAppliedFrequencyCorrection); - mRequestedCorrectAveragingBuffer.reset(); + if(partialCorrection != 0) + { + mAppliedFrequencyCorrection += (long)(partialCorrection * ERROR_PROCESSING_GAIN); + mTunerChannelSource.setFrequencyCorrection(mAppliedFrequencyCorrection); + mRequestedCorrectAveragingBuffer.reset(); + } + + appliedFrequencyCorrection = mAppliedFrequencyCorrection; } ChannelFrequencyCorrectionStatusNotification notification = ChannelFrequencyCorrectionStatusNotification - .create(mTunerChannelSource, mAverageRequestedCorrection, mAppliedFrequencyCorrection, mParent.getTunerPPM(), + .create(mTunerChannelSource, mAverageRequestedCorrection, appliedFrequencyCorrection, mParent.getTunerPPM(), mParent.getTunerFrequencyCorrection(), mParent.isEnabled()); mTunerChannelSource.broadcastConsumerSourceEvent(notification); mAverageRequestedCorrection = 0; @@ -152,7 +172,7 @@ public class ChannelFrequencyErrorManager implements ISourceEventListener, Liste else { ChannelFrequencyCorrectionStatusNotification notification = ChannelFrequencyCorrectionStatusNotification - .create(mTunerChannelSource, 0, mAppliedFrequencyCorrection, mParent.getTunerPPM(), + .create(mTunerChannelSource, 0, getFrequencyError(), mParent.getTunerPPM(), mParent.getTunerFrequencyCorrection(), mParent.isEnabled()); mTunerChannelSource.broadcastConsumerSourceEvent(notification); } diff --git a/src/main/java/io/github/dsheirer/source/tuner/frequency/TunerFrequencyErrorManager.java b/src/main/java/io/github/dsheirer/source/tuner/frequency/TunerFrequencyErrorManager.java index e1c01799..789c5694 100644 --- a/src/main/java/io/github/dsheirer/source/tuner/frequency/TunerFrequencyErrorManager.java +++ b/src/main/java/io/github/dsheirer/source/tuner/frequency/TunerFrequencyErrorManager.java @@ -132,6 +132,13 @@ public class TunerFrequencyErrorManager implements ISourceEventProcessor } } + //Skip the update when no channel has reported a measurement since the last update. A divide by zero + //exception would cancel all subsequent executions of this scheduled task. + if(count == 0) + { + return; + } + requestedChangeHz /= count; mTunerController.setMeasuredFrequencyError((int)requestedChangeHz); diff --git a/src/test/java/io/github/dsheirer/buffer/FloatNativeBufferTest.java b/src/test/java/io/github/dsheirer/buffer/FloatNativeBufferTest.java new file mode 100644 index 00000000..1d8759e0 --- /dev/null +++ b/src/test/java/io/github/dsheirer/buffer/FloatNativeBufferTest.java @@ -0,0 +1,109 @@ +/* + * ***************************************************************************** + * 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.buffer; + +import io.github.dsheirer.sample.complex.ComplexSamples; +import io.github.dsheirer.sample.complex.InterleavedComplexSamples; +import java.util.ArrayList; +import java.util.Iterator; +import java.util.List; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the float native buffer. + */ +public class FloatNativeBufferTest +{ + private static final int SAMPLE_COUNT = 5000; + + /** + * Creates interleaved complex samples where each value is unique and non-zero. + */ + private static float[] getSamples() + { + float[] samples = new float[SAMPLE_COUNT * 2]; + + for(int x = 0; x < samples.length; x++) + { + samples[x] = x + 1; + } + + return samples; + } + + /** + * Test: the complex samples iterator returns every sample exactly once, with a final chunk that holds the remaining + * samples rather than padding them to the full chunk length with zero valued samples. + */ + @Test + void testComplexSamplesIterator() + { + float[] samples = getSamples(); + List lengths = new ArrayList<>(); + int pointer = 0; + + Iterator iterator = new FloatNativeBuffer(samples, 0L, 0.0f).iterator(); + + while(iterator.hasNext()) + { + ComplexSamples chunk = iterator.next(); + lengths.add(chunk.i().length); + + for(int x = 0; x < chunk.i().length; x++) + { + assertEquals(samples[pointer++], chunk.i()[x], "Inphase sample " + (pointer / 2)); + assertEquals(samples[pointer++], chunk.q()[x], "Quadrature sample " + (pointer / 2)); + } + } + + assertEquals(samples.length, pointer, "Sample count"); + assertEquals(List.of(2048, 2048, 904), lengths, "Chunk lengths"); + } + + /** + * Test: the interleaved complex samples iterator returns every sample exactly once, without zero padding. + */ + @Test + void testInterleavedComplexSamplesIterator() + { + float[] samples = getSamples(); + List lengths = new ArrayList<>(); + int pointer = 0; + + Iterator iterator = new FloatNativeBuffer(samples, 0L, 0.0f).iteratorInterleaved(); + + while(iterator.hasNext()) + { + float[] chunk = iterator.next().samples(); + lengths.add(chunk.length / 2); + + for(float sample : chunk) + { + assertEquals(samples[pointer], sample, "Sample " + pointer); + pointer++; + } + } + + assertEquals(samples.length, pointer, "Sample count"); + assertEquals(List.of(2048, 2048, 904), lengths, "Chunk lengths"); + } +} diff --git a/src/test/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliabilityTest.java b/src/test/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliabilityTest.java new file mode 100644 index 00000000..969c23cd --- /dev/null +++ b/src/test/java/io/github/dsheirer/dsp/fm/FMDiscriminatorReliabilityTest.java @@ -0,0 +1,95 @@ +/* + * ***************************************************************************** + * 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.fm; + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertArrayEquals; +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the FM discriminator soft symbol reliability. + */ +public class FMDiscriminatorReliabilityTest +{ + private static final int BUFFER_LENGTH = 256; + + /** + * Test: inverse power filter taps are the squared differences of the post-discriminator filter taps. + */ + @Test + void testInversePowerFilterTaps() + { + float[] taps = FMDiscriminatorReliability.getInversePowerFilterTaps(new float[]{0.25f, 0.5f, 0.25f}); + assertArrayEquals(new float[]{0.0625f, 0.0625f, 0.0625f, 0.0625f}, taps, 1e-7f, "Taps"); + + //An integrate and dump (flat) filter depends only on the samples at either end of the filter span + float[] flat = FMDiscriminatorReliability.getInversePowerFilterTaps(new float[]{0.2f, 0.2f, 0.2f, 0.2f, 0.2f}); + assertArrayEquals(new float[]{0.04f, 0f, 0f, 0f, 0f, 0.04f}, flat, 1e-7f, "Flat filter taps"); + } + + /** + * Test: each reliability is 1 / sum(g[k]^2 / |x[n - k]|^2), including across the sample buffer boundary. + */ + @Test + void testReliabilities() + { + float[] filter = {0.1f, 0.2f, 0.4f, 0.2f, 0.1f}; + float[] g2 = FMDiscriminatorReliability.getInversePowerFilterTaps(filter); + FMDiscriminatorReliability reliability = new FMDiscriminatorReliability(filter); + float[] amplitudes = new float[2 * BUFFER_LENGTH]; + + for(int x = 0; x < amplitudes.length; x++) + { + amplitudes[x] = 0.1f + 0.9f * (float)Math.abs(Math.sin(x * 0.37)); + } + + float[] reliabilities = new float[amplitudes.length]; + + for(int buffer = 0; buffer < 2; buffer++) + { + float[] i = new float[BUFFER_LENGTH]; + float[] q = new float[BUFFER_LENGTH]; + + for(int x = 0; x < BUFFER_LENGTH; x++) + { + //Rotate the phase to show that the reliability depends only on the sample power + double phase = x * 0.5; + float amplitude = amplitudes[buffer * BUFFER_LENGTH + x]; + i[x] = (float)(amplitude * Math.cos(phase)); + q[x] = (float)(amplitude * Math.sin(phase)); + } + + System.arraycopy(reliability.calculate(i, q), 0, reliabilities, buffer * BUFFER_LENGTH, BUFFER_LENGTH); + } + + for(int x = g2.length; x < amplitudes.length; x++) + { + double sum = 0.0; + + for(int k = 0; k < g2.length; k++) + { + sum += g2[k] / (amplitudes[x - k] * amplitudes[x - k]); + } + + assertEquals(1.0 / sum, reliabilities[x], 1e-4 / sum, "Reliability at " + x); + } + } +} diff --git a/src/test/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25Test.java b/src/test/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25Test.java new file mode 100644 index 00000000..e071c960 --- /dev/null +++ b/src/test/java/io/github/dsheirer/edac/trellis/ViterbiDecoder_1_2_P25Test.java @@ -0,0 +1,238 @@ +/* + * ***************************************************************************** + * 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.edac.trellis; + +import io.github.dsheirer.bits.CorrectedBinaryMessage; +import io.github.dsheirer.dsp.symbol.Dibit; +import java.util.Random; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * JUnit tests for the P25 1/2 rate trellis soft decision Viterbi decoder. + */ +public class ViterbiDecoder_1_2_P25Test +{ + private static final int INPUT_DIBITS = 48; + //Faded symbol spacing for the reliability weighting test (ie one third of the symbols are faded) + private static final int SPACING = 3; + private static final Dibit[] DIBITS = {Dibit.D00_PLUS_1, Dibit.D01_PLUS_3, Dibit.D10_MINUS_1, Dibit.D11_MINUS_3}; + + /** + * Encodes the input dibit values (plus a flushing zero) into 98 ideal soft symbols. + */ + private static float[] encode(int[] inputs) + { + float[] symbols = new float[(inputs.length + 1) * 2]; + int state = 0; + + for(int step = 0; step <= inputs.length; step++) + { + int input = step < inputs.length ? inputs[step] : 0; + int output = P25_1_2_Node.TRANSITION_MATRIX[state][input]; + symbols[2 * step] = DIBITS[output >> 2].getIdealPhase(); + symbols[2 * step + 1] = DIBITS[output & 3].getIdealPhase(); + state = input; + } + + return symbols; + } + + /** + * Counts the decoded message bits that differ from the input dibit values. + */ + private static int countBitErrors(int[] inputs, CorrectedBinaryMessage decoded) + { + int errors = 0; + + for(int step = 0; step < inputs.length; step++) + { + if(decoded.get(2 * step) != ((inputs[step] & 2) == 2)) + { + errors++; + } + + if(decoded.get(2 * step + 1) != ((inputs[step] & 1) == 1)) + { + errors++; + } + } + + return errors; + } + + private static int[] randomInputs(Random random) + { + int[] inputs = new int[INPUT_DIBITS]; + + for(int x = 0; x < inputs.length; x++) + { + inputs[x] = random.nextInt(4); + } + + return inputs; + } + + /** + * Test: ideal soft symbols decode without error in both the linear (C4FM) and phase wrapped (LSM) modes. + */ + @Test + void testIdealSymbols() + { + ViterbiDecoder_1_2_P25 decoder = new ViterbiDecoder_1_2_P25(); + Random random = new Random(1); + + for(int trial = 0; trial < 50; trial++) + { + int[] inputs = randomInputs(random); + float[] symbols = encode(inputs); + assertEquals(0, countBitErrors(inputs, decoder.decodeSoft(symbols)), "Linear"); + assertEquals(0, countBitErrors(inputs, decoder.decodeSoft(symbols, true)), "Phase wrapped"); + } + } + + /** + * Test: phase wrapped decoding measures the soft symbol error around the circle, so outer symbols that noise + * pushed across the +/- PI phase wrap are decoded correctly. + */ + @Test + void testOuterSymbolsAcrossPhaseWrap() + { + ViterbiDecoder_1_2_P25 decoder = new ViterbiDecoder_1_2_P25(); + Random random = new Random(2); + int linearErrors = 0; + int wrappedErrors = 0; + + for(int trial = 0; trial < 200; trial++) + { + int[] inputs = randomInputs(random); + float[] symbols = encode(inputs); + int outerCount = 0; + + for(int x = 0; x < symbols.length; x++) + { + //Push every third outer (+/- 3PI/4) symbol 0.85 radians toward and across the phase wrap + if(Math.abs(symbols[x]) > Math.PI / 2 && (outerCount++ % 3) == 0) + { + double pushed = symbols[x] + Math.signum(symbols[x]) * 0.85; + symbols[x] = (float)(pushed - 2.0 * Math.PI * Math.round(pushed / (2.0 * Math.PI))); + } + } + + linearErrors += countBitErrors(inputs, decoder.decodeSoft(symbols)); + wrappedErrors += countBitErrors(inputs, decoder.decodeSoft(symbols, true)); + } + + assertEquals(0, wrappedErrors, "Phase wrapped bit errors"); + assertTrue(wrappedErrors <= linearErrors, "Phase wrapped decoding no worse than linear decoding"); + } + + /** + * Test: reliability weighting decodes corrupted symbols that have low reliability that unweighted soft decision + * decoding fails to decode. The corrupted symbols are spaced across the block, the way the deinterleaver spreads + * the symbols of a transmitted envelope fade. + */ + @Test + void testReliabilityWeightedFade() + { + ViterbiDecoder_1_2_P25 decoder = new ViterbiDecoder_1_2_P25(); + Random random = new Random(4); + int unweightedErrors = 0; + int weightedErrors = 0; + + for(int trial = 0; trial < 50; trial++) + { + int[] inputs = randomInputs(random); + float[] symbols = encode(inputs); + float[] reliabilities = new float[symbols.length]; + java.util.Arrays.fill(reliabilities, 1.0f); + //Replace every SPACING symbol with a random phase (ie a deep fade) and mark it with low reliability + for(int x = random.nextInt(SPACING); x < symbols.length; x += SPACING) + { + symbols[x] = (float)((random.nextDouble() * 2.0 - 1.0) * Math.PI); + reliabilities[x] = 0.01f; + } + + unweightedErrors += countBitErrors(inputs, decoder.decodeSoft(symbols, true)); + weightedErrors += countBitErrors(inputs, decoder.decodeSoft(symbols, reliabilities, true)); + } + + assertEquals(0, weightedErrors, "Reliability weighted bit errors"); + assertTrue(unweightedErrors > 0, "Unweighted decoding fails to decode the faded symbols"); + } + + /** + * Test: unavailable or unusable reliabilities (null, NaN or all zero) decode the same as unweighted decoding. + */ + @Test + void testUnusableReliabilities() + { + ViterbiDecoder_1_2_P25 decoder = new ViterbiDecoder_1_2_P25(); + Random random = new Random(5); + + for(int trial = 0; trial < 20; trial++) + { + int[] inputs = randomInputs(random); + float[] symbols = encode(inputs); + + for(int x = 0; x < symbols.length; x++) + { + symbols[x] += (float)(random.nextGaussian() * 0.5); + } + + CorrectedBinaryMessage unweighted = decoder.decodeSoft(symbols, true); + float[] nan = new float[symbols.length]; + java.util.Arrays.fill(nan, 1.0f); + nan[7] = Float.NaN; + assertEquals(unweighted, decoder.decodeSoft(symbols, null, true), "Null reliabilities"); + assertEquals(unweighted, decoder.decodeSoft(symbols, nan, true), "NaN reliability"); + assertEquals(unweighted, decoder.decodeSoft(symbols, new float[symbols.length], true), "Zero reliabilities"); + } + } + + /** + * Test: soft symbols that are not phase wrapped (C4FM) with outer symbols beyond +/- PI decode correctly with + * linear error measurement. + */ + @Test + void testUnwrappedOuterSymbolsBeyondPi() + { + ViterbiDecoder_1_2_P25 decoder = new ViterbiDecoder_1_2_P25(); + Random random = new Random(3); + + for(int trial = 0; trial < 50; trial++) + { + int[] inputs = randomInputs(random); + float[] symbols = encode(inputs); + + for(int x = 0; x < symbols.length; x++) + { + if(Math.abs(symbols[x]) > Math.PI / 2) + { + symbols[x] += (float)(Math.signum(symbols[x]) * 0.85); + } + } + + assertEquals(0, countBitErrors(inputs, decoder.decodeSoft(symbols)), "Linear"); + } + } +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/NACTrackerTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/NACTrackerTest.java new file mode 100644 index 00000000..21e1b5b7 --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/NACTrackerTest.java @@ -0,0 +1,75 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the NAC tracker. + */ +public class NACTrackerTest +{ + /** + * Test: the dominant NAC needs the minimum number of observations. + */ + @Test + void testDominantNAC() + { + NACTracker tracker = new NACTracker(); + assertEquals(0, tracker.getTrackedNAC(), "Empty"); + tracker.track(0x293); + tracker.track(0x293); + tracker.track(0x111); + assertEquals(0, tracker.getTrackedNAC(), "Below the minimum observations"); + tracker.track(0x293); + assertEquals(0x293, tracker.getTrackedNAC(), "Dominant"); + } + + /** + * Test: when a fourth NAC is observed, the least recently observed NAC is evicted, independent of the system clock. + * All of the observations occur within the same millisecond, as they do when a recording decodes faster than real + * time. + */ + @Test + void testLeastRecentlyObservedEvicted() + { + for(int repeat = 0; repeat < 100; repeat++) + { + NACTracker tracker = new NACTracker(); + + //0x293 is the most observed, but the least recently observed when 0x444 arrives + for(int x = 0; x < 5; x++) + { + tracker.track(0x293); + } + + for(int x = 0; x < 4; x++) + { + tracker.track(0x111); + } + + tracker.track(0x222); + tracker.track(0x444); + assertEquals(0x111, tracker.getTrackedNAC(), "Least recently observed NAC evicted, repeat " + repeat); + } + } +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSMTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSMTest.java new file mode 100644 index 00000000..b82bd51f --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DecoderLSMTest.java @@ -0,0 +1,409 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import io.github.dsheirer.dsp.symbol.Dibit; +import io.github.dsheirer.edac.CRCP25; +import io.github.dsheirer.edac.bch.BCH_63_16_23_P25; +import io.github.dsheirer.edac.trellis.P25_1_2_Node; +import io.github.dsheirer.message.IMessage; +import io.github.dsheirer.module.decode.p25.phase1.message.tsbk.TSBKMessage; +import io.github.dsheirer.sample.complex.ComplexSamples; +import java.util.ArrayList; +import java.util.BitSet; +import java.util.HashSet; +import java.util.List; +import java.util.Random; +import java.util.Set; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * End-to-end JUnit tests for the P25 Phase 1 LSM decoder: a synthetic LSM (CQPSK) control channel carrying single + * block TSBK frames is generated and decoded, and each valid decoded TSBK is checked against the transmitted TSBKs. + */ +public class P25P1DecoderLSMTest +{ + private static final double SAMPLE_RATE = 48000.0; + private static final double SYMBOL_RATE = 4800.0; + private static final int FRAME_COUNT = 150; + private static final int NAC = 0x293; + //Single block TSBK frame: 180 dibits (360 bits) including a status dibit every 36th dibit and null padding + private static final int FRAME_DIBITS = 180; + private static final int TSBK_DUID = 7; + private static final int OPCODE_RFSS_STATUS_BROADCAST = 58; + //Raised cosine transmit pulse + private static final double ROLLOFF = 0.2; + private static final int PULSE_SPAN_SYMBOLS = 8; + private static final int BUFFER_SIZE = 2048; + private static final Dibit[] DIBITS = {Dibit.D00_PLUS_1, Dibit.D01_PLUS_3, Dibit.D10_MINUS_1, Dibit.D11_MINUS_3}; + + /** + * Test: a clean signal decodes every TSBK frame, with no invalid messages. + */ + @Test + void testCleanSignal() + { + Result result = decode(0.0, 0.0, Double.NaN, 1, BUFFER_SIZE); + assertEquals(0, result.invalid(), "Invalid messages"); + assertEquals(FRAME_COUNT, result.valid(), "Valid messages"); + } + + /** + * Test: a carrier offset within the +/- 600 Hz fourth power offset estimator range and a 100 ppm symbol rate error + * decode every TSBK frame. + */ + @Test + void testCarrierOffsetAndSymbolRateError() + { + Result result = decode(300.0, 100.0, Double.NaN, 2, BUFFER_SIZE); + assertEquals(0, result.invalid(), "Invalid messages"); + assertEquals(FRAME_COUNT, result.valid(), "Valid messages"); + } + + /** + * Test: carrier offsets beyond the fourth power offset estimator range are resolved by the offset ambiguity search, + * after which every TSBK frame decodes. The search steps toward the adjacent offset range nearest to the estimate, + * so offsets on either side acquire after one search dwell period (~11 frames). + */ + @Test + void testOffsetAmbiguitySearch() + { + for(double carrierOffset : new double[]{700.0, -700.0, 900.0, -900.0}) + { + Result result = decode(carrierOffset, -100.0, Double.NaN, 3, BUFFER_SIZE); + assertEquals(0, result.invalid(), "Invalid messages at " + carrierOffset + " Hz"); + assertTrue(result.first() <= 16, "Acquisition at " + carrierOffset + " Hz: " + result.first() + " frames"); + assertEquals(0, result.missed(), "Missed frames after acquisition at " + carrierOffset + " Hz"); + } + } + + /** + * Test: acquisition of offsets near the edge of the fourth power offset estimator range is fast and independent of + * the sample buffer size. The estimator previously updated once per buffer, so a -550 Hz offset took 36 frames to + * acquire with 2048 sample buffers and 55 frames with 8192 sample buffers. + */ + @Test + void testAcquisitionIndependentOfBufferSize() + { + for(int bufferSize : new int[]{512, 2048, 8192}) + { + for(double carrierOffset : new double[]{550.0, -550.0}) + { + Result result = decode(carrierOffset, 50.0, Double.NaN, 5, bufferSize); + String condition = carrierOffset + " Hz with " + bufferSize + " sample buffers"; + assertEquals(0, result.invalid(), "Invalid messages at " + condition); + assertTrue(result.first() <= 5, "Acquisition at " + condition + ": " + result.first() + " frames"); + assertEquals(0, result.missed(), "Missed frames after acquisition at " + condition); + } + } + } + + /** + * Test: with additive white gaussian noise at 12 dB in-band (+/- 3 kHz) carrier to noise, nearly all TSBK frames + * decode and every valid TSBK matches a transmitted TSBK. + */ + @Test + void testNoise() + { + Result result = decode(200.0, 50.0, 12.0, 4, BUFFER_SIZE); + assertTrue(result.valid() >= FRAME_COUNT - 5, "Valid messages: " + result.valid()); + } + + /** + * Generates and decodes the synthetic signal. + * @param carrierOffset in Hertz + * @param symbolRateError in parts per million + * @param carrierToNoise in dB, in-band (+/- 3 kHz), or NaN for no noise + * @param seed for the random TSBK payloads and noise + * @param bufferSize of the sample buffers delivered to the decoder (even) + * @return decoded message counts, the index of the first decoded TSBK frame (the frame count when + * none decode) and the count of TSBK frames that + * didn't decode after the first decoded TSBK frame. + */ + private static Result decode(double carrierOffset, double symbolRateError, double carrierToNoise, long seed, + int bufferSize) + { + Random random = new Random(seed); + List tsbks = new ArrayList<>(); + Set transmitted = new HashSet<>(); + List dibits = new ArrayList<>(); + + for(int frame = 0; frame < FRAME_COUNT; frame++) + { + BitSet tsbk = getTSBK(random); + tsbks.add(tsbk); + transmitted.add(tsbk); + dibits.addAll(getFrame(tsbk)); + } + + //Trailing frame that isn't checked: the decoder dispatches a message after its trailing symbols are processed + dibits.addAll(getFrame(getTSBK(random))); + + float[][] samples = modulate(dibits, carrierOffset, symbolRateError, carrierToNoise, random); + + List received = new ArrayList<>(); + P25P1DecoderLSM decoder = new P25P1DecoderLSM(); + decoder.setMessageListener((IMessage message) -> { + if(message instanceof TSBKMessage tsbk) + { + received.add(tsbk); + } + }); + decoder.start(); + decoder.setSampleRate(SAMPLE_RATE); + + for(int offset = 0; offset < samples[0].length; offset += bufferSize) + { + int length = Math.min(bufferSize, samples[0].length - offset); + float[] i = new float[length]; + float[] q = new float[length]; + System.arraycopy(samples[0], offset, i, 0, length); + System.arraycopy(samples[1], offset, q, 0, length); + decoder.receive(new ComplexSamples(i, q, 0L)); + } + + decoder.stop(); + + int valid = 0; + int invalid = 0; + Set decoded = new HashSet<>(); + + for(TSBKMessage tsbk : received) + { + if(tsbk.isValid()) + { + BitSet bits = new BitSet(); + + for(int x = 0; x < 80; x++) + { + bits.set(x, tsbk.getMessage().get(x)); + } + + //A valid TSBK that doesn't match a transmitted TSBK is a CRC false positive + assertTrue(transmitted.contains(bits), "Decoded TSBK was not transmitted"); + assertTrue(decoded.add(bits), "Decoded TSBK was duplicated"); + valid++; + } + else + { + invalid++; + } + } + + int first = tsbks.size(); + int missed = 0; + + for(int x = 0; x < tsbks.size(); x++) + { + if(decoded.contains(tsbks.get(x))) + { + first = Math.min(first, x); + } + else if(first < x) + { + missed++; + } + } + + return new Result(valid, invalid, first, missed); + } + + /** + * Creates a random single block (last block) TSBK, 80 message bits. + */ + private static BitSet getTSBK(Random random) + { + BitSet tsbk = new BitSet(); + tsbk.set(0); //Last block + setValue(tsbk, 2, 6, OPCODE_RFSS_STATUS_BROADCAST); + + for(int x = 16; x < 80; x++) + { + tsbk.set(x, random.nextBoolean()); + } + + return tsbk; + } + + /** + * Creates the frame dibits for the TSBK: sync, NID, interleaved trellis encoded TSBK and null padding, with a + * status dibit at every 36th dibit. + */ + private static List getFrame(BitSet tsbk) + { + List data = new ArrayList<>(); + addBits(data, 0x5575F5FF77FFL, 48); + addBits(data, BCH_63_16_23_P25.encode(NAC, TSBK_DUID) << 1, 64); //64th (parity) bit is zero + data.addAll(encodeTSBK(tsbk)); + + List frame = new ArrayList<>(); + int pointer = 0; + + for(int x = 1; x <= FRAME_DIBITS; x++) + { + if(x % 36 == 0) + { + frame.add(Dibit.D01_PLUS_3); //Outbound (repeater) status + } + else + { + frame.add(pointer < data.size() ? data.get(pointer++) : Dibit.D00_PLUS_1); + } + } + + return frame; + } + + /** + * Adds the CRC-CCITT checksum, 1/2 rate trellis encodes and interleaves the TSBK into 98 dibits. + */ + private static List encodeTSBK(BitSet tsbk) + { + BitSet message = (BitSet)tsbk.clone(); + int checksum = 0xFFFF; + + for(int x = message.nextSetBit(0); x >= 0 && x < 80; x = message.nextSetBit(x + 1)) + { + checksum ^= CRCP25.CCITT_80_CHECKSUMS[x]; + } + + setValue(message, 80, 16, checksum); + + //48 input dibits plus a flushing zero input produce 98 encoded dibits + int[] encoded = new int[98]; + int state = 0; + + for(int step = 0; step < 49; step++) + { + int input = step < 48 ? (message.get(2 * step) ? 2 : 0) + (message.get(2 * step + 1) ? 1 : 0) : 0; + int output = P25_1_2_Node.TRANSITION_MATRIX[state][input]; + encoded[2 * step] = output >> 2; + encoded[2 * step + 1] = output & 3; + state = input; + } + + List interleaved = new ArrayList<>(); + + for(int dibit = 0; dibit < encoded.length; dibit++) + { + interleaved.add(DIBITS[encoded[P25P1Interleave.DATA_DEINTERLEAVE[2 * dibit] / 2]]); + } + + return interleaved; + } + + /** + * Modulates the dibits as LSM: differential QPSK with raised cosine pulse shaping. + * @return i and q sample arrays + */ + private static float[][] modulate(List dibits, double carrierOffset, double symbolRateError, + double carrierToNoise, Random random) + { + double[] symbolI = new double[dibits.size()]; + double[] symbolQ = new double[dibits.size()]; + double phase = 0.0; + + for(int x = 0; x < dibits.size(); x++) + { + phase += dibits.get(x).getIdealPhase(); + symbolI[x] = Math.cos(phase); + symbolQ[x] = Math.sin(phase); + } + + double samplesPerSymbol = SAMPLE_RATE / (SYMBOL_RATE * (1.0 + symbolRateError * 1e-6)); + //The decimation filter requires an even sample count + int length = (int)((dibits.size() + PULSE_SPAN_SYMBOLS) * samplesPerSymbol) / 2 * 2; + float[] i = new float[length]; + float[] q = new float[length]; + + //Unit power symbols: the in-band noise variance is the in-band fraction of the total noise variance + double noise = Double.isNaN(carrierToNoise) ? 0.0 : + Math.sqrt(Math.pow(10.0, -carrierToNoise / 10.0) * SAMPLE_RATE / 6000.0 / 2.0); + + for(int n = 0; n < length; n++) + { + double t = n / samplesPerSymbol - PULSE_SPAN_SYMBOLS / 2.0; //Symbol periods + double sumI = 0.0; + double sumQ = 0.0; + + for(int k = (int)Math.ceil(t - PULSE_SPAN_SYMBOLS / 2.0); k <= t + PULSE_SPAN_SYMBOLS / 2.0; k++) + { + if(k >= 0 && k < dibits.size()) + { + double pulse = getRaisedCosine(t - k); + sumI += symbolI[k] * pulse; + sumQ += symbolQ[k] * pulse; + } + } + + double carrier = 2.0 * Math.PI * carrierOffset * n / SAMPLE_RATE; + double cos = Math.cos(carrier); + double sin = Math.sin(carrier); + i[n] = (float)(sumI * cos - sumQ * sin + noise * random.nextGaussian()); + q[n] = (float)(sumI * sin + sumQ * cos + noise * random.nextGaussian()); + } + + return new float[][]{i, q}; + } + + /** + * Raised cosine pulse value + * @param t time in symbol periods + */ + private static double getRaisedCosine(double t) + { + double sinc = Math.abs(t) < 1e-9 ? 1.0 : Math.sin(Math.PI * t) / (Math.PI * t); + double denominator = 1.0 - Math.pow(2.0 * ROLLOFF * t, 2.0); + + if(Math.abs(denominator) < 1e-9) + { + return Math.PI / 4.0 * sinc; + } + + return sinc * Math.cos(Math.PI * ROLLOFF * t) / denominator; + } + + /** + * Adds the most significant bits first, as dibits. + */ + private static void addBits(List dibits, long value, int bitCount) + { + for(int bit = bitCount - 2; bit >= 0; bit -= 2) + { + dibits.add(DIBITS[(int)((value >> bit) & 3)]); + } + } + + /** + * Sets the value at the bit offset, most significant bit first. + */ + private static void setValue(BitSet bits, int offset, int bitCount, int value) + { + for(int x = 0; x < bitCount; x++) + { + bits.set(offset + x, ((value >> (bitCount - 1 - x)) & 1) == 1); + } + } + + private record Result(int valid, int invalid, int first, int missed) {} +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSMTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSMTest.java new file mode 100644 index 00000000..7c2a25a1 --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1DifferentialDemodulatorLSMTest.java @@ -0,0 +1,353 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import io.github.dsheirer.dsp.filter.FilterFactory; +import io.github.dsheirer.dsp.symbol.Dibit; +import java.util.Random; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * JUnit tests for the P25 Phase 1 LSM differential demodulator and differential integrator filter design. + */ +public class P25P1DifferentialDemodulatorLSMTest +{ + private static final int SAMPLE_RATE = 24000; + private static final int SAMPLES_PER_SYMBOL = 5; + private static final int BUFFER_LENGTH = 1024; + private static final Dibit[] DIBITS = {Dibit.D00_PLUS_1, Dibit.D01_PLUS_3, Dibit.D10_MINUS_1, Dibit.D11_MINUS_3}; + + /** + * Synthetic differential QPSK signal generator. Ramps the phase across each symbol period by the symbol's ideal + * phase so that the phase change across each symbol period is the ideal symbol phase, plus a carrier offset. + */ + private static class SignalGenerator + { + private final Random mRandom = new Random(1); + private final double mOffset; //Radians per sample + private final double mNoise; + private double mPhase = 0.0; + private Dibit mSymbol = Dibit.D00_PLUS_1; + private int mSampleCounter = 0; + private Dibit mPreviousSymbol; + + SignalGenerator(double offsetHertz, double noise) + { + mOffset = 2.0 * Math.PI * offsetHertz / SAMPLE_RATE; + mNoise = noise; + } + + /** + * Generates the next buffer of samples, recording the symbol completed at each sample (or null). + */ + float[][] next(Dibit[] completedSymbols) + { + float[] i = new float[BUFFER_LENGTH]; + float[] q = new float[BUFFER_LENGTH]; + + for(int x = 0; x < BUFFER_LENGTH; x++) + { + if(mSampleCounter % SAMPLES_PER_SYMBOL == 0) + { + mPreviousSymbol = mSymbol; + mSymbol = DIBITS[mRandom.nextInt(4)]; + completedSymbols[x] = mSampleCounter > 0 ? mPreviousSymbol : null; + } + else + { + completedSymbols[x] = null; + } + + i[x] = (float)(Math.cos(mPhase) + mRandom.nextGaussian() * mNoise); + q[x] = (float)(Math.sin(mPhase) + mRandom.nextGaussian() * mNoise); + mPhase += mSymbol.getIdealPhase() / SAMPLES_PER_SYMBOL + mOffset; + mSampleCounter++; + } + + return new float[][]{i, q}; + } + } + + private static double toHertz(float radiansPerSample) + { + return radiansPerSample * SAMPLE_RATE / (2.0 * Math.PI); + } + + /** + * Test: the integrator taps span one symbol period with the fractional remainder in the final tap and unity DC gain. + */ + @Test + void testDifferentialIntegratorDesign() + { + float[] taps = FilterFactory.getDifferentialIntegrator(5.2); + assertEquals(6, taps.length, "Tap count"); + float sum = 0; + + for(float tap : taps) + { + sum += tap; + } + + assertEquals(1.0, sum, 1e-5, "DC gain"); + assertEquals(1.0 / 5.2, taps[0], 1e-6, "Whole tap"); + assertEquals(0.2 / 5.2, taps[5], 1e-6, "Fractional tap"); + assertEquals(5, FilterFactory.getDifferentialIntegrator(5.0).length, "Integral samples per symbol"); + } + + /** + * Test: the carrier offset estimate converges and each symbol period's demodulated sample, scaled to radians per + * symbol, decides the transmitted symbol with the offset removed. + */ + @Test + void testOffsetRemovalAndSymbolRecovery() + { + SignalGenerator generator = new SignalGenerator(300, 0.0); + P25P1DifferentialDemodulatorLSM demodulator = new P25P1DifferentialDemodulatorLSM(); + demodulator.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + Dibit[] completed = new Dibit[BUFFER_LENGTH]; + int symbolCount = 0; + int errorCount = 0; + + for(int buffer = 0; buffer < 400; buffer++) + { + float[][] samples = generator.next(completed); + float[] demodulated = demodulator.demodulate(samples[0], samples[1], true); + + //Check symbol decisions after the estimator converges + if(buffer >= 300) + { + for(int x = 0; x < BUFFER_LENGTH; x++) + { + if(completed[x] != null) + { + symbolCount++; + + if(Dibit.fromSample(demodulated[x] * SAMPLES_PER_SYMBOL) != completed[x]) + { + errorCount++; + } + } + } + } + } + + assertEquals(300.0, toHertz(demodulator.getFrequencyOffset()), 5.0, "Offset estimate"); + assertTrue(symbolCount > 10000, "Symbol count"); + assertEquals(0, errorCount, "Symbol errors"); + } + + /** + * Test: an offset beyond the +/- 600 Hz fourth power estimator range is found by the ambiguity search while + * unsynchronized, and is held once synchronized. + */ + @Test + void testAmbiguitySearch() + { + SignalGenerator generator = new SignalGenerator(-750, 0.1); + P25P1DifferentialDemodulatorLSM demodulator = new P25P1DifferentialDemodulatorLSM(); + demodulator.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + Dibit[] completed = new Dibit[BUFFER_LENGTH]; + boolean resolved = false; + + for(int buffer = 0; buffer < 1000 && !resolved; buffer++) + { + float[][] samples = generator.next(completed); + demodulator.demodulate(samples[0], samples[1], false); + + //Emulate the symbol processor synchronizing once the offset is resolved + resolved = Math.abs(toHertz(demodulator.getFrequencyOffset()) + 750) < 20; + } + + assertTrue(resolved, "Search resolves the offset"); + + for(int buffer = 0; buffer < 200; buffer++) + { + float[][] samples = generator.next(completed); + demodulator.demodulate(samples[0], samples[1], true); + } + + assertEquals(-750.0, toHertz(demodulator.getFrequencyOffset()), 10.0, "Offset held while synchronized"); + } + + /** + * Test: the ambiguity search holds the resolved offset range across a gap with no signal while unsynchronized (ie + * between transmissions on a traffic channel) so that the next transmission starts in the correct offset range. + */ + @Test + void testSearchHoldsAcrossSignalGap() + { + SignalGenerator generator = new SignalGenerator(-750, 0.1); + P25P1DifferentialDemodulatorLSM demodulator = new P25P1DifferentialDemodulatorLSM(); + demodulator.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + Dibit[] completed = new Dibit[BUFFER_LENGTH]; + boolean resolved = false; + + for(int buffer = 0; buffer < 1000 && !resolved; buffer++) + { + float[][] samples = generator.next(completed); + demodulator.demodulate(samples[0], samples[1], false); + resolved = Math.abs(toHertz(demodulator.getFrequencyOffset()) + 750) < 20; + } + + assertTrue(resolved, "Search resolves the offset"); + + //Emulate the symbol processor synchronizing to the signal and the framer decoding a valid NID once the offset is + //resolved + for(int buffer = 0; buffer < 20; buffer++) + { + float[][] samples = generator.next(completed); + demodulator.demodulate(samples[0], samples[1], true); + demodulator.offsetVerified(); + } + + //Gap of ~20,000 symbols (10x the search dwell period) with no signal and no synchronization + Random random = new Random(3); + + for(int buffer = 0; buffer < 100; buffer++) + { + float[] i = new float[BUFFER_LENGTH]; + float[] q = new float[BUFFER_LENGTH]; + + for(int x = 0; x < BUFFER_LENGTH; x++) + { + i[x] = (float)random.nextGaussian(); + q[x] = (float)random.nextGaussian(); + } + + demodulator.demodulate(i, q, false); + } + + assertEquals(-750.0, toHertz(demodulator.getFrequencyOffset()), 20.0, "Offset held across the gap"); + } + + /** + * Test: symbol processor synchronization without a valid NID decode (ie a false sync detection in an incorrect + * offset range) doesn't hold the ambiguity search in the incorrect offset range for the longer dwell period, while + * a verified offset range is held. + */ + @Test + void testUnverifiedSynchronizationDoesNotHoldSearch() + { + SignalGenerator generatorUnverified = new SignalGenerator(-750, 0.1); + SignalGenerator generatorVerified = new SignalGenerator(-750, 0.1); + P25P1DifferentialDemodulatorLSM unverified = new P25P1DifferentialDemodulatorLSM(); + P25P1DifferentialDemodulatorLSM verified = new P25P1DifferentialDemodulatorLSM(); + unverified.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + verified.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + Dibit[] completed = new Dibit[BUFFER_LENGTH]; + + //Emulate a sync detection in the initial (incorrect) offset range, with and without a valid NID decode + for(int buffer = 0; buffer < 5; buffer++) + { + float[][] samples = generatorUnverified.next(completed); + unverified.demodulate(samples[0], samples[1], true); + samples = generatorVerified.next(completed); + verified.demodulate(samples[0], samples[1], true); + verified.offsetVerified(); + } + + //Unsynchronized for ~10,000 symbols: longer than the search dwell and shorter than the verified dwell + for(int buffer = 0; buffer < 50; buffer++) + { + float[][] samples = generatorUnverified.next(completed); + unverified.demodulate(samples[0], samples[1], false); + samples = generatorVerified.next(completed); + verified.demodulate(samples[0], samples[1], false); + } + + //The unverified search steps to an adjacent (1200 Hz) offset range while the verified search holds its range + double difference = Math.abs(toHertz(unverified.getFrequencyOffset()) - toHertz(verified.getFrequencyOffset())); + assertTrue(difference > 600, "Unverified offset range search stepped away [" + difference + " Hz]"); + } + + /** + * Test: each sample's reliability is (|x[n]|^2 * |x[n - sps]|^2) / (|x[n]|^2 + |x[n - sps]|^2), including across the + * sample buffer boundary. + */ + @Test + void testReliabilities() + { + P25P1DifferentialDemodulatorLSM demodulator = new P25P1DifferentialDemodulatorLSM(); + demodulator.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + float[] amplitudes = new float[2 * BUFFER_LENGTH]; + + for(int x = 0; x < amplitudes.length; x++) + { + amplitudes[x] = 0.1f + 0.9f * (float)Math.abs(Math.sin(x * 0.37)); + } + + float[] reliabilities = new float[amplitudes.length]; + + for(int buffer = 0; buffer < 2; buffer++) + { + float[] i = new float[BUFFER_LENGTH]; + float[] q = new float[BUFFER_LENGTH]; + + for(int x = 0; x < BUFFER_LENGTH; x++) + { + //Rotate the phase to show that the reliability depends only on the sample power + double phase = x * 0.5; + float amplitude = amplitudes[buffer * BUFFER_LENGTH + x]; + i[x] = (float)(amplitude * Math.cos(phase)); + q[x] = (float)(amplitude * Math.sin(phase)); + } + + demodulator.demodulate(i, q, true); + System.arraycopy(demodulator.getReliabilities(), 0, reliabilities, buffer * BUFFER_LENGTH, BUFFER_LENGTH); + } + + for(int x = SAMPLES_PER_SYMBOL; x < amplitudes.length; x++) + { + double current = amplitudes[x] * amplitudes[x]; + double delayed = amplitudes[x - SAMPLES_PER_SYMBOL] * amplitudes[x - SAMPLES_PER_SYMBOL]; + assertEquals(current * delayed / (current + delayed), reliabilities[x], 1e-5, "Reliability at " + x); + } + } + + /** + * Test: the offset estimate doesn't drift when the input is noise. + */ + @Test + void testNoiseDoesNotMoveEstimate() + { + Random random = new Random(2); + P25P1DifferentialDemodulatorLSM demodulator = new P25P1DifferentialDemodulatorLSM(); + demodulator.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + + for(int buffer = 0; buffer < 2000; buffer++) + { + float[] i = new float[BUFFER_LENGTH]; + float[] q = new float[BUFFER_LENGTH]; + + for(int x = 0; x < BUFFER_LENGTH; x++) + { + i[x] = (float)random.nextGaussian(); + q[x] = (float)random.nextGaussian(); + } + + demodulator.demodulate(i, q, true); + } + + assertEquals(0.0, toHertz(demodulator.getFrequencyOffset()), 1e-3, "Offset estimate"); + } +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssemblerTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssemblerTest.java new file mode 100644 index 00000000..0f0d5a61 --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1MessageAssemblerTest.java @@ -0,0 +1,138 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import io.github.dsheirer.dsp.symbol.Dibit; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the P25 Phase 1 message assembler forced completion of place holder (failed NID) assemblies. + */ +public class P25P1MessageAssemblerTest +{ + private static final P25P1DataUnitID TDU = P25P1DataUnitID.TERMINATOR_DATA_UNIT; + private static final P25P1DataUnitID TDULC = P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL; + private static final P25P1DataUnitID TSBK = P25P1DataUnitID.TRUNKING_SIGNALING_BLOCK_1; + + /** + * Creates a place holder assembly with the bit length. + */ + private static P25P1MessageAssembler getPlaceHolder(int bitLength) + { + P25P1MessageAssembler assembler = new P25P1MessageAssembler(0x293, P25P1DataUnitID.PLACE_HOLDER); + + for(int x = 0; x < bitLength / 2; x++) + { + assembler.receive(Dibit.D01_PLUS_3); + } + + return assembler; + } + + /** + * Forces completion of a place holder with the bit length. + * @return the identified data unit ID + */ + private static P25P1DataUnitID forceCompletion(int bitLength, P25P1DataUnitID previous, P25P1DataUnitID next) + { + P25P1MessageAssembler assembler = getPlaceHolder(bitLength); + assembler.forceCompletion(previous, next); + return assembler.getDataUnitID(); + } + + /** + * Test: a place holder that's shorter than a missed TDU frame (ie inter-frame padding or timing slack ahead of the + * next sync) isn't identified as a TDU, even when followed by an LDU1 (HDU resolution) or a TSBK. + */ + @Test + void testShortPlaceHolderDiscarded() + { + for(int length : new int[]{0, 2, 28, 100, 126}) + { + P25P1MessageAssembler assembler = getPlaceHolder(length); + assertEquals(0, assembler.forceCompletion(TDULC, TDULC), "Dropped bits at " + length); + assertEquals(P25P1DataUnitID.PLACE_HOLDER, assembler.getDataUnitID(), "Length " + length); + assertEquals(P25P1DataUnitID.PLACE_HOLDER, forceCompletion(length, P25P1DataUnitID.HEADER_DATA_UNIT, + P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_1), "Length " + length + " before LDU1"); + } + } + + /** + * Test: a place holder with the length of a single missed TDU frame (28 null bits plus the next frame's 112 sync and + * NID bits) is identified as a TDU. + */ + @Test + void testMissedTDU() + { + for(int length : new int[]{128, 136, 140}) + { + assertEquals(TDU, forceCompletion(length, TDU, TDU), "Length " + length); + assertEquals(TDU, forceCompletion(length, TDULC, TDU), "Length " + length); + } + + assertEquals(P25P1DataUnitID.PLACE_HOLDER, forceCompletion(144, TDU, TDU), "Length 144"); + } + + /** + * Test: a place holder with the length of a single missed TDULC frame (308 message and null bits plus the next + * frame's 112 sync and NID bits) is identified as a TDULC. + */ + @Test + void testMissedTDULC() + { + for(int length : new int[]{408, 416, 420}) + { + assertEquals(TDULC, forceCompletion(length, TDULC, TDULC), "Length " + length); + } + } + + /** + * Test: a place holder that doesn't match a missed frame isn't identified as a TDU. + */ + @Test + void testUnidentifiedPlaceHolderDiscarded() + { + for(int length : new int[]{208, 276, 556, 696, 1000}) + { + P25P1MessageAssembler assembler = getPlaceHolder(length); + assertEquals(0, assembler.forceCompletion(TDULC, TDU), "Dropped bits at " + length); + assertEquals(P25P1DataUnitID.PLACE_HOLDER, assembler.getDataUnitID(), "Length " + length); + } + + //A place holder between voice frames that is too short for an LDU + assertEquals(P25P1DataUnitID.PLACE_HOLDER, forceCompletion(1000, P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_2, + P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_1), "Length 1000 between LDUs"); + } + + /** + * Test: place holders are resolved from the next data unit ID before the missed frame lengths are applied, so a + * place holder ahead of a TSBK is a TSBK, even when its length matches a missed TDULC frame. + */ + @Test + void testResolvedFromNextDataUnit() + { + assertEquals(TSBK, forceCompletion(350, TSBK, TSBK), "Missed TSBK"); + assertEquals(TSBK, forceCompletion(416, TSBK, TSBK), "Missed TSBK with TDULC length"); + assertEquals(P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_2, forceCompletion(1560, + P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_1, P25P1DataUnitID.LOGICAL_LINK_DATA_UNIT_1), "Missed LDU2"); + } +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSMTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSMTest.java new file mode 100644 index 00000000..bc639856 --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/P25P1SymbolProcessorLSMTest.java @@ -0,0 +1,80 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1; + +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the P25 Phase 1 LSM symbol processor. + */ +public class P25P1SymbolProcessorLSMTest +{ + private static final float SAMPLES_PER_SYMBOL = 5.0f; + private static final double TOLERANCE = 1e-5; + + private static P25P1SymbolProcessorLSM getProcessor() + { + P25P1SymbolProcessorLSM processor = new P25P1SymbolProcessorLSM(new P25P1DecoderLSM(), new P25P1MessageFramer()); + processor.setSamplesPerSymbol(SAMPLES_PER_SYMBOL); + return processor; + } + + /** + * Converts radians per symbol to the radians per sample domain of the demodulated sample stream. + */ + private static float toSample(double radiansPerSymbol) + { + return (float)(radiansPerSymbol / SAMPLES_PER_SYMBOL); + } + + /** + * Test: samples that don't straddle the phase wrap are linearly interpolated. + */ + @Test + void testLinearInterpolation() + { + P25P1SymbolProcessorLSM processor = getProcessor(); + assertEquals(toSample(1.5), processor.interpolate(toSample(1.0), toSample(2.0), 0.5), TOLERANCE); + assertEquals(toSample(-1.25), processor.interpolate(toSample(-1.0), toSample(-2.0), 0.25), TOLERANCE); + } + + /** + * Test: samples that straddle the phase wrap near an outer symbol interpolate across the wrap instead of through + * zero, and the result is wrapped to +/- PI radians per symbol. + */ + @Test + void testInterpolationAcrossPhaseWrap() + { + P25P1SymbolProcessorLSM processor = getProcessor(); + + //+3.0 and -3.0 are 0.283 radians apart across the wrap at +/- PI. + assertEquals(toSample(Math.PI), Math.abs(processor.interpolate(toSample(3.0), toSample(-3.0), 0.5)), TOLERANCE); + + //A quarter of the way from +3.0 toward -3.0 (across the wrap) is still on the positive side of the wrap. + double expected = 3.0 + 0.25 * (2.0 * Math.PI - 6.0); + assertEquals(toSample(expected), processor.interpolate(toSample(3.0), toSample(-3.0), 0.25), TOLERANCE); + + //Three quarters of the way crosses the wrap to the negative side. + expected = 3.0 + 0.75 * (2.0 * Math.PI - 6.0) - 2.0 * Math.PI; + assertEquals(toSample(expected), processor.interpolate(toSample(3.0), toSample(-3.0), 0.75), TOLERANCE); + } +} diff --git a/src/test/java/io/github/dsheirer/module/decode/p25/phase1/message/P25P1MessageValidityTest.java b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/message/P25P1MessageValidityTest.java new file mode 100644 index 00000000..51e0b5bd --- /dev/null +++ b/src/test/java/io/github/dsheirer/module/decode/p25/phase1/message/P25P1MessageValidityTest.java @@ -0,0 +1,117 @@ +/* + * ***************************************************************************** + * 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.module.decode.p25.phase1.message; + +import io.github.dsheirer.bits.CorrectedBinaryMessage; +import io.github.dsheirer.module.decode.p25.phase1.P25P1DataUnitID; +import io.github.dsheirer.module.decode.p25.phase1.message.hdu.HDUMessage; +import io.github.dsheirer.module.decode.p25.phase1.message.tdu.TDULCMessage; +import java.util.Random; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertTrue; + +/** + * JUnit tests for the message validity of P25 Phase 1 messages that decode their error protected content on demand. + * The validity must reflect the error correction result on the first isValid() call, before any other accessor has + * decoded the content. + */ +public class P25P1MessageValidityTest +{ + /** + * Creates a message with all-zero bits (a valid codeword for the linear Golay and Reed-Solomon codes). + */ + private static CorrectedBinaryMessage getZeroMessage(P25P1DataUnitID dataUnitID) + { + CorrectedBinaryMessage message = new CorrectedBinaryMessage(dataUnitID.getMessageLength()); + message.setSize(dataUnitID.getMessageLength()); + return message; + } + + /** + * Creates a message with random bits. + */ + private static CorrectedBinaryMessage getRandomMessage(P25P1DataUnitID dataUnitID, long seed) + { + Random random = new Random(seed); + CorrectedBinaryMessage message = getZeroMessage(dataUnitID); + + for(int x = 0; x < message.size(); x++) + { + message.set(x, random.nextBoolean()); + } + + return message; + } + + /** + * Test: a TDULC is valid when its link control word passes error correction. + */ + @Test + void testValidTDULC() + { + TDULCMessage tdulc = new TDULCMessage(getZeroMessage(P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL), + 0x293, 0L); + assertTrue(tdulc.isValid(), "TDULC"); + assertTrue(tdulc.getLinkControlWord().isValid(), "Link control word"); + } + + /** + * Test: a TDULC is invalid when its link control word fails error correction. + */ + @Test + void testInvalidTDULC() + { + for(long seed = 1; seed <= 10; seed++) + { + TDULCMessage tdulc = new TDULCMessage(getRandomMessage( + P25P1DataUnitID.TERMINATOR_DATA_UNIT_LINK_CONTROL, seed), 0x293, 0L); + assertFalse(tdulc.isValid(), "TDULC with random bits, seed " + seed); + assertFalse(tdulc.getLinkControlWord().isValid(), "Link control word, seed " + seed); + } + } + + /** + * Test: an HDU is valid when its header data passes error correction. + */ + @Test + void testValidHDU() + { + HDUMessage hdu = new HDUMessage(getZeroMessage(P25P1DataUnitID.HEADER_DATA_UNIT), 0x293, 0L); + assertTrue(hdu.isValid(), "HDU"); + assertTrue(hdu.getHeaderData().isValid(), "Header data"); + } + + /** + * Test: an HDU is invalid when its header data fails error correction, including when isValid() is invoked before + * the header data is accessed. + */ + @Test + void testInvalidHDU() + { + for(long seed = 1; seed <= 10; seed++) + { + HDUMessage hdu = new HDUMessage(getRandomMessage(P25P1DataUnitID.HEADER_DATA_UNIT, seed), 0x293, 0L); + assertFalse(hdu.isValid(), "HDU with random bits, seed " + seed); + assertFalse(hdu.getHeaderData().isValid(), "Header data, seed " + seed); + } + } +} diff --git a/src/test/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManagerTest.java b/src/test/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManagerTest.java new file mode 100644 index 00000000..88622485 --- /dev/null +++ b/src/test/java/io/github/dsheirer/source/tuner/frequency/ChannelFrequencyErrorManagerTest.java @@ -0,0 +1,104 @@ +/* + * ***************************************************************************** + * 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.source.tuner.frequency; + +import io.github.dsheirer.sample.Listener; +import io.github.dsheirer.sample.complex.ComplexSamples; +import io.github.dsheirer.source.tuner.channel.TunerChannel; +import io.github.dsheirer.source.tuner.channel.TunerChannelSource; +import org.junit.jupiter.api.Test; + +import static org.junit.jupiter.api.Assertions.assertEquals; + +/** + * JUnit tests for the channel frequency error manager. + */ +public class ChannelFrequencyErrorManagerTest +{ + /** + * Test: a tuner PPM change is removed from the channel's applied frequency correction, and the updated correction + * is applied to the channel source. + */ + @Test + void testTunerPPMChange() + { + TestChannelSource source = new TestChannelSource(); + ChannelFrequencyErrorManager manager = new ChannelFrequencyErrorManager(source, null); + + manager.receive(-10); + assertEquals(10, manager.getFrequencyError(), "Applied correction"); + assertEquals(10, source.mFrequencyCorrection, "Channel source correction"); + + manager.receive(25); + assertEquals(-15, manager.getFrequencyError(), "Applied correction"); + assertEquals(-15, source.mFrequencyCorrection, "Channel source correction"); + } + + /** + * Channel source that records the applied frequency correction. + */ + private static class TestChannelSource extends TunerChannelSource + { + private long mFrequencyCorrection = Long.MIN_VALUE; + + public TestChannelSource() + { + super(null, new TunerChannel(460_000_000L, 12_500), "test", null); + } + + @Override + public void setFrequencyCorrection(long correction) + { + mFrequencyCorrection = correction; + } + + @Override + public void setFrequency(long frequency) + { + } + + @Override + protected void setSampleRate(double sampleRate) + { + } + + @Override + public void setListener(Listener complexSamplesListener) + { + } + + @Override + public double getSampleRate() + { + return 50_000.0; + } + + @Override + public long getFrequency() + { + return 460_000_000L; + } + + @Override + public void reset() + { + } + } +}