mirror of
https://github.com/f4exb/sdrangel.git
synced 2026-08-11 14:03:56 -04:00
Use MLSE demod to improve sensitivity.
Normalise threshold. Fix frequency deviation.
This commit is contained in:
Jon Beniston
@@ -583,8 +583,8 @@ void AISDemodGUI::on_fmDev_valueChanged(int value)
|
||||
|
||||
void AISDemodGUI::on_threshold_valueChanged(int value)
|
||||
{
|
||||
ui->thresholdText->setText(QString("%1").arg(value));
|
||||
m_settings.m_correlationThreshold = value;
|
||||
m_settings.m_correlationThreshold = value / 100.0f;
|
||||
ui->thresholdText->setText(QString("%1").arg(m_settings.m_correlationThreshold, 0, 'f', 2));
|
||||
applySettings(QStringList({"correlationThreshold"}));
|
||||
}
|
||||
|
||||
@@ -898,8 +898,8 @@ void AISDemodGUI::displaySettings()
|
||||
ui->fmDevText->setText(QString("%1%2k").arg(QChar(0xB1, 0x00)).arg(m_settings.m_fmDeviation / 1000.0, 0, 'f', 1));
|
||||
ui->fmDev->setValue(m_settings.m_fmDeviation / 100.0);
|
||||
|
||||
ui->thresholdText->setText(QString("%1").arg(m_settings.m_correlationThreshold));
|
||||
ui->threshold->setValue(m_settings.m_correlationThreshold);
|
||||
ui->thresholdText->setText(QString("%1").arg(m_settings.m_correlationThreshold, 0, 'f', 2));
|
||||
ui->threshold->setValue((int) std::round(m_settings.m_correlationThreshold * 100.0f));
|
||||
|
||||
updateIndexLabel();
|
||||
|
||||
|
||||
@@ -359,13 +359,13 @@
|
||||
</size>
|
||||
</property>
|
||||
<property name="toolTip">
|
||||
<string>Correlation threshold</string>
|
||||
<string>Normalised correlation threshold with the preamble, from 0 to 1</string>
|
||||
</property>
|
||||
<property name="minimum">
|
||||
<number>0</number>
|
||||
</property>
|
||||
<property name="maximum">
|
||||
<number>60</number>
|
||||
<number>100</number>
|
||||
</property>
|
||||
<property name="singleStep">
|
||||
<number>1</number>
|
||||
@@ -380,8 +380,14 @@
|
||||
</item>
|
||||
<item>
|
||||
<widget class="QLabel" name="thresholdText">
|
||||
<property name="minimumSize">
|
||||
<size>
|
||||
<width>28</width>
|
||||
<height>0</height>
|
||||
</size>
|
||||
</property>
|
||||
<property name="text">
|
||||
<string>60</string>
|
||||
<string>0.60</string>
|
||||
</property>
|
||||
</widget>
|
||||
</item>
|
||||
|
||||
@@ -37,8 +37,8 @@ void AISDemodSettings::resetToDefaults()
|
||||
m_baud = AISDEMOD_BAUD_RATE; // Fixed
|
||||
m_inputFrequencyOffset = 0;
|
||||
m_rfBandwidth = 16000.0f;
|
||||
m_fmDeviation = 4800.0f;
|
||||
m_correlationThreshold = 30;
|
||||
m_fmDeviation = 2400.0f;
|
||||
m_correlationThreshold = 0.6f;
|
||||
m_filterMMSI = "";
|
||||
m_udpEnabled = false;
|
||||
m_udpAddress = "127.0.0.1";
|
||||
@@ -73,7 +73,6 @@ QByteArray AISDemodSettings::serialize() const
|
||||
s.writeS32(1, m_inputFrequencyOffset);
|
||||
s.writeFloat(2, m_rfBandwidth);
|
||||
s.writeFloat(3, m_fmDeviation);
|
||||
s.writeFloat(4, m_correlationThreshold);
|
||||
s.writeString(5, m_filterMMSI);
|
||||
s.writeBool(6, m_udpEnabled);
|
||||
s.writeString(7, m_udpAddress);
|
||||
@@ -108,6 +107,7 @@ QByteArray AISDemodSettings::serialize() const
|
||||
s.writeBool(28, m_hidden);
|
||||
s.writeBool(29, m_showSlotMap);
|
||||
s.writeBool(30, m_useFileTime);
|
||||
s.writeFloat(32, m_correlationThreshold);
|
||||
|
||||
for (int i = 0; i < AISDEMOD_MESSAGE_COLUMNS; i++)
|
||||
s.writeS32(100 + i, m_messageColumnIndexes[i]);
|
||||
@@ -135,8 +135,7 @@ bool AISDemodSettings::deserialize(const QByteArray& data)
|
||||
|
||||
d.readS32(1, &m_inputFrequencyOffset, 0);
|
||||
d.readFloat(2, &m_rfBandwidth, 16000.0f);
|
||||
d.readFloat(3, &m_fmDeviation, 4800.0f);
|
||||
d.readFloat(4, &m_correlationThreshold, 30);
|
||||
d.readFloat(3, &m_fmDeviation, 2400.0f);
|
||||
d.readString(5, &m_filterMMSI, "");
|
||||
d.readBool(6, &m_udpEnabled);
|
||||
d.readString(7, &m_udpAddress, "127.0.0.1");
|
||||
@@ -196,6 +195,8 @@ bool AISDemodSettings::deserialize(const QByteArray& data)
|
||||
d.readBool(29, &m_showSlotMap, false);
|
||||
d.readBool(30, &m_useFileTime, false);
|
||||
|
||||
d.readFloat(32, &m_correlationThreshold, 0.6f);
|
||||
|
||||
for (int i = 0; i < AISDEMOD_MESSAGE_COLUMNS; i++) {
|
||||
d.readS32(100 + i, &m_messageColumnIndexes[i], i);
|
||||
}
|
||||
|
||||
@@ -36,8 +36,11 @@ struct AISDemodSettings
|
||||
qint32 m_baud;
|
||||
qint32 m_inputFrequencyOffset;
|
||||
Real m_rfBandwidth;
|
||||
Real m_fmDeviation; //!< Peak deviation to give modulation index of 0.5 for 9600 baud
|
||||
Real m_correlationThreshold;
|
||||
Real m_fmDeviation; //!< Peak deviation. M.1371-5 2.3.2 specifies a modulation index of
|
||||
//!< 0.5, and h = 2.dev/baud, so the peak deviation is 2400 Hz at
|
||||
//!< 9600 baud. 4800 Hz is the mark to space separation, not the
|
||||
//!< deviation
|
||||
Real m_correlationThreshold; //!< Normalised correlation with the preamble, 0 to 1
|
||||
QString m_filterMMSI;
|
||||
bool m_udpEnabled;
|
||||
QString m_udpAddress;
|
||||
@@ -74,6 +77,13 @@ struct AISDemodSettings
|
||||
static const int AISDEMOD_CHANNEL_SAMPLE_RATE = 57600; //!< 6x 9600 baud rate (use even multiple so Gaussian filter has odd number of taps)
|
||||
static const int m_scopeStreams = 9;
|
||||
|
||||
//! The trellis models the transmitted waveform, so this is the transmit BT-product of
|
||||
//! 0.4 from M.1371-5 2.3.1.2 - not the 0.5 receive BT-product of 2.3.1.3, which is what
|
||||
//! m_pulseShape uses for the preamble correlator. They are different numbers on purpose.
|
||||
//! A 3 symbol phase pulse gives a 16 state trellis; 4 measures no better and costs twice
|
||||
static constexpr float AISDEMOD_MLSE_BT = 0.4f;
|
||||
static const int AISDEMOD_MLSE_SPAN = 3;
|
||||
|
||||
AISDemodSettings();
|
||||
void resetToDefaults();
|
||||
void setChannelMarker(Serializable *channelMarker) { m_channelMarker = channelMarker; }
|
||||
|
||||
@@ -1,6 +1,7 @@
|
||||
///////////////////////////////////////////////////////////////////////////////////
|
||||
// Copyright (C) 2021, 2023 Jon Beniston, M7RCE <jon@beniston.com> //
|
||||
// Copyright (C) 2021-2022 Edouard Griffiths, F4EXB <f4exb06@gmail.com> //
|
||||
// Some code by AI //
|
||||
// //
|
||||
// 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 //
|
||||
@@ -40,6 +41,7 @@ AISDemodSink::AISDemodSink(AISDemod *aisDemod) :
|
||||
m_messageQueueToChannel(nullptr),
|
||||
m_rxBuf(nullptr),
|
||||
m_train(nullptr),
|
||||
m_trainEnergy(0.0f),
|
||||
m_sampleBufferIndex(0)
|
||||
{
|
||||
m_magsq = 0.0;
|
||||
@@ -116,6 +118,58 @@ void AISDemodSink::feed(const SampleVector::const_iterator& begin, const SampleV
|
||||
}
|
||||
}
|
||||
|
||||
// Per burst carrier frequency offset, radians per sample.
|
||||
//
|
||||
// The training sequence has equal numbers of ones and zeros, so its mean phase increment
|
||||
// is the carrier offset and nothing else. Summing the products before taking the argument
|
||||
// keeps the estimate out of the noisy per sample phase.
|
||||
//
|
||||
// This only has to get the MLSE's phase tracking loop into its pull in range - the loop
|
||||
// removes what is left, including any drift the estimate cannot see.
|
||||
double AISDemodSink::estimateFrequency() const
|
||||
{
|
||||
std::complex<double> sum(0.0, 0.0);
|
||||
|
||||
for (int i = 1; i < m_correlationLength; i++)
|
||||
{
|
||||
int j = (m_rxBufIdx + i) % m_rxBufLength;
|
||||
int k = (j + m_rxBufLength - 1) % m_rxBufLength;
|
||||
sum += m_iqBuf[j] * std::conj(m_iqBuf[k]);
|
||||
}
|
||||
|
||||
return std::arg(sum);
|
||||
}
|
||||
|
||||
// Run the sequence detector over n symbols starting at buffer position x, leaving the
|
||||
// symbol decisions in m_symSoft
|
||||
void AISDemodSink::computeSymbols(int x, int n)
|
||||
{
|
||||
const double dphi = estimateFrequency();
|
||||
const int sps = m_samplesPerSymbol;
|
||||
const int half = sps / 2;
|
||||
const int len = m_rxBufLength;
|
||||
const std::complex<double> *iq = m_iqBuf.data();
|
||||
|
||||
// The phase loop starts with no idea of the carrier phase, so run it over some of the
|
||||
// preamble first and throw those symbols away. x sits 18 symbols into the 24 bit
|
||||
// training sequence, so there is room for this - but only that much, so clamp rather
|
||||
// than reading backwards off the start of the buffer.
|
||||
const int xOffset = ((x - m_rxBufIdx) % len + len) % len;
|
||||
const int warmup = std::min(AISDEMOD_MLSE_WARMUP, (xOffset - half) / sps);
|
||||
|
||||
m_mlse.decode(n + warmup,
|
||||
[x, sps, half, len, iq, dphi, warmup](int k, int m) -> std::complex<double>
|
||||
{
|
||||
int off = (k - warmup)*sps - half + m;
|
||||
int idx = ((x + off) % len + len) % len;
|
||||
double a = -dphi * off;
|
||||
return iq[idx] * std::complex<double>(cos(a), sin(a));
|
||||
},
|
||||
m_mlseSoft);
|
||||
|
||||
m_symSoft.assign(m_mlseSoft.begin() + warmup, m_mlseSoft.end());
|
||||
}
|
||||
|
||||
void AISDemodSink::processOneSample(Complex &ci)
|
||||
{
|
||||
// FM demodulation
|
||||
@@ -150,10 +204,12 @@ void AISDemodSink::processOneSample(Complex &ci)
|
||||
// Buffer filtered samples. We buffer enough samples for a max length message
|
||||
// before trying to demod, so false triggering can't make us miss anything
|
||||
m_rxBuf[m_rxBufIdx] = filtClipped;
|
||||
m_iqBuf[m_rxBufIdx] = std::complex<double>(ci.real() / SDR_RX_SCALEF, ci.imag() / SDR_RX_SCALEF);
|
||||
m_rxBufIdx = (m_rxBufIdx + 1) % m_rxBufLength;
|
||||
m_rxBufCnt = std::min(m_rxBufCnt + 1, m_rxBufLength);
|
||||
|
||||
Real corr = 0.0f;
|
||||
Real metric = 0.0f;
|
||||
bool scopeCRCValid = false;
|
||||
bool scopeCRCInvalid = false;
|
||||
Real dcOffset = 0.0f;
|
||||
@@ -161,6 +217,7 @@ void AISDemodSink::processOneSample(Complex &ci)
|
||||
if (m_rxBufCnt >= m_rxBufLength)
|
||||
{
|
||||
Real trainingSum = 0.0f;
|
||||
Real energy = 0.0f;
|
||||
|
||||
// Correlate with training sequence
|
||||
// Note that DC offset doesn't matter for this
|
||||
@@ -170,182 +227,78 @@ void AISDemodSink::processOneSample(Complex &ci)
|
||||
int j = (m_rxBufIdx + i) % m_rxBufLength;
|
||||
corr += m_train[i] * m_rxBuf[j];
|
||||
trainingSum += m_rxBuf[j];
|
||||
energy += m_rxBuf[j] * m_rxBuf[j];
|
||||
}
|
||||
|
||||
// If we meet threshold, try to demod
|
||||
// Take abs value, to account for both initial phases
|
||||
thresholdMet = fabs(corr) >= m_settings.m_correlationThreshold;
|
||||
// Take abs value, to account for both initial phases.
|
||||
// Dividing by the geometric mean of the two energies gives a correlation
|
||||
// coefficient in 0..1, which unlike the raw correlation does not depend on the
|
||||
// signal level. That matters because the sequence detector is too expensive to run
|
||||
// on the false triggers an absolute threshold lets through - one fires on noise for
|
||||
// several percent of all samples.
|
||||
metric = (Real) (fabs(corr) / sqrt((double) energy * m_trainEnergy + 1e-12));
|
||||
thresholdMet = metric >= m_settings.m_correlationThreshold;
|
||||
|
||||
if (thresholdMet)
|
||||
{
|
||||
// Use mean of preamble as DC offset
|
||||
// Mean of the preamble, which is the carrier offset. Only reported to the
|
||||
// scope now - the sequence detector estimates and tracks it for itself.
|
||||
dcOffset = trainingSum/m_correlationLength;
|
||||
|
||||
// Start demod after (most of) preamble
|
||||
int x = (m_rxBufIdx + m_correlationLength*3/4 + 4) % m_rxBufLength;
|
||||
|
||||
// Attempt to demodulate
|
||||
bool gotSOP = false;
|
||||
int bits = 0;
|
||||
int bitCount = 0;
|
||||
int onesCount = 0;
|
||||
int byteCount = 0;
|
||||
int symbolPrev = 0;
|
||||
int totalBitCount = 0; // Count of bits after start flag, before bit stuffing removal, including stop flag
|
||||
for (int sampleIdx = 0; sampleIdx < m_rxBufLength; sampleIdx += m_samplesPerSymbol)
|
||||
int endSampleIdx = 0;
|
||||
|
||||
// Decode only far enough to tell whether there is a start flag. Nearly all
|
||||
// triggers are noise and abort here, which is what makes the sequence
|
||||
// detector affordable.
|
||||
computeSymbols(x, AISDEMOD_SOP_SYMBOLS + 8);
|
||||
|
||||
DemodResult result = deframe(endSampleIdx, scopeCRCValid, scopeCRCInvalid);
|
||||
|
||||
// Extend on demand. Most messages fit one slot, but type 5 and the binary
|
||||
// messages take two or more, and decoding every candidate to the full
|
||||
// buffer length would cost far more than decoding the few long ones twice.
|
||||
//
|
||||
// The look ahead must not run off the end of the circular buffer: symbol n-1
|
||||
// reads up to n*samplesPerSymbol - samplesPerSymbol/2 - 1 samples beyond x,
|
||||
// and anything past the end wraps onto pre-trigger samples, which the
|
||||
// Viterbi traceback would then start from. So cap the span at what is
|
||||
// genuinely buffered, and make sure that capped span is actually attempted
|
||||
// rather than being skipped by the doubling.
|
||||
int xOffset = ((x - m_rxBufIdx) % m_rxBufLength + m_rxBufLength) % m_rxBufLength;
|
||||
int maxSymbols = (m_rxBufLength - xOffset + m_samplesPerSymbol/2) / m_samplesPerSymbol;
|
||||
|
||||
// Never extend to fewer symbols than the first stage already decoded, or the
|
||||
// deframer would run out before it reaches the start flag it just found
|
||||
const int firstSpan = std::max(AISDEMOD_MLSE_SYMBOLS, AISDEMOD_SOP_SYMBOLS + 8);
|
||||
|
||||
for (int span = firstSpan; result == FrameTruncated; span *= 2)
|
||||
{
|
||||
// Sum and slice
|
||||
// Summing 3 samples seems to give a very small improvement vs just using 1
|
||||
int sampleCnt = 3;
|
||||
int sampleOffset = -1;
|
||||
Real sampleSum = 0.0f;
|
||||
for (int i = 0; i < sampleCnt; i++) {
|
||||
sampleSum += m_rxBuf[(x + sampleOffset + i + m_rxBufLength) % m_rxBufLength] - dcOffset;
|
||||
}
|
||||
int symbol = sampleSum >= 0.0f ? 1 : 0;
|
||||
int symbols = std::min(span, maxSymbols);
|
||||
|
||||
// Move to next symbol
|
||||
x = (x + m_samplesPerSymbol) % m_rxBufLength;
|
||||
computeSymbols(x, symbols);
|
||||
result = deframe(endSampleIdx, scopeCRCValid, scopeCRCInvalid);
|
||||
|
||||
// HDLC deframing
|
||||
|
||||
// NRZI decoding
|
||||
int bit;
|
||||
if (symbol != symbolPrev) {
|
||||
bit = 0;
|
||||
} else {
|
||||
bit = 1;
|
||||
}
|
||||
symbolPrev = symbol;
|
||||
|
||||
// Store in shift reg
|
||||
bits |= bit << bitCount;
|
||||
bitCount++;
|
||||
|
||||
if (bit == 1)
|
||||
{
|
||||
onesCount++;
|
||||
// Shouldn't ever get 7 1s in a row
|
||||
if ((onesCount == 7) && gotSOP)
|
||||
{
|
||||
gotSOP = false;
|
||||
byteCount = 0;
|
||||
break;
|
||||
}
|
||||
}
|
||||
else if (bit == 0)
|
||||
{
|
||||
if (onesCount == 5)
|
||||
{
|
||||
// Remove bit-stuffing (5 1s followed by a 0)
|
||||
bitCount--;
|
||||
}
|
||||
else if (onesCount == 6)
|
||||
{
|
||||
// Start/end of packet
|
||||
if (gotSOP && (bitCount == 8) && (bits == 0x7e) && (byteCount >= 3))
|
||||
{
|
||||
// End of packet
|
||||
// Check CRC is valid
|
||||
m_crc.init();
|
||||
m_crc.calculate(m_bytes, byteCount - 2);
|
||||
uint16_t calcCrc = m_crc.get();
|
||||
uint16_t rxCrc = m_bytes[byteCount-2] | (m_bytes[byteCount-1] << 8);
|
||||
if (calcCrc == rxCrc)
|
||||
{
|
||||
scopeCRCValid = true;
|
||||
QByteArray rxPacket((char *)m_bytes, byteCount - 2); // Don't include CRC
|
||||
//qDebug() << "RX: " << rxPacket.toHex();
|
||||
if (getMessageQueueToChannel())
|
||||
{
|
||||
// Calculate slot number based on time of start of transmission
|
||||
// This is unlikely to be accurate in absolute terms, given we don't know latency from SDR or buffering within SDRangel
|
||||
// But can be used to get an idea of congestion
|
||||
QDateTime currentTime = QDateTime::currentDateTime();
|
||||
if (m_settings.m_useFileTime)
|
||||
{
|
||||
QString hardwareId = m_aisDemod->getDeviceAPI()->getHardwareId();
|
||||
|
||||
if ((hardwareId == "FileInput") || (hardwareId == "SigMFFileInput"))
|
||||
{
|
||||
QString dateTimeStr;
|
||||
int deviceIdx = m_aisDemod->getDeviceSetIndex();
|
||||
|
||||
if (ChannelWebAPIUtils::getDeviceReportValue(deviceIdx, "absoluteTime", dateTimeStr)) {
|
||||
currentTime = QDateTime::fromString(dateTimeStr, Qt::ISODateWithMs);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int txTimeMs = (totalBitCount + 8 + 24 + 8) * (1000.0 / m_settings.m_baud); // Add ramp up, preamble and start-flag
|
||||
QDateTime startDateTime = currentTime.addMSecs(-txTimeMs);
|
||||
int ms = startDateTime.time().second() * 1000 + startDateTime.time().msec();
|
||||
float slotTime = 60.0f * 1000.0f / 2250.0f; // 2250 slots per minute, 26.6ms per slot
|
||||
int slot = ms / slotTime;
|
||||
int totalSlots = std::ceil(txTimeMs / slotTime);
|
||||
AISDemod::MsgMessage *msg = AISDemod::MsgMessage::create(rxPacket, currentTime, slot, totalSlots);
|
||||
getMessageQueueToChannel()->push(msg);
|
||||
}
|
||||
|
||||
// Skip over received packet, so we don't try to re-demodulate it
|
||||
m_rxBufCnt -= sampleIdx;
|
||||
}
|
||||
else
|
||||
{
|
||||
//qDebug() << QString("CRC mismatch: %1 %2").arg(calcCrc, 4, 16, QLatin1Char('0')).arg(rxCrc, 4, 16, QLatin1Char('0'));
|
||||
scopeCRCInvalid = true;
|
||||
}
|
||||
break;
|
||||
}
|
||||
else if (gotSOP)
|
||||
{
|
||||
// Repeated start flag without data or misalignment, something not right
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Start of packet
|
||||
gotSOP = true;
|
||||
bits = 0;
|
||||
bitCount = 0;
|
||||
byteCount = 0;
|
||||
totalBitCount = 0;
|
||||
}
|
||||
}
|
||||
onesCount = 0;
|
||||
}
|
||||
|
||||
if (gotSOP)
|
||||
{
|
||||
totalBitCount++;
|
||||
if (bitCount == 8)
|
||||
{
|
||||
// Could also check count according to message ID as that varies
|
||||
if (byteCount >= AISDEMOD_MAX_BYTES)
|
||||
{
|
||||
// Too many bytes
|
||||
break;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Got a complete byte
|
||||
m_bytes[byteCount] = bits;
|
||||
byteCount++;
|
||||
}
|
||||
bits = 0;
|
||||
bitCount = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Abort demod if we haven't found start flag within a couple of bytes of presumed preamble
|
||||
if (!gotSOP && (sampleIdx >= 16 * m_samplesPerSymbol)) {
|
||||
if (symbols >= maxSymbols) {
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (result == FrameGood)
|
||||
{
|
||||
// Skip over received packet, so we don't try to re-demodulate it
|
||||
m_rxBufCnt -= endSampleIdx;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Select signals to feed to scope
|
||||
sampleToScope(ci / SDR_RX_SCALEF, magsq, fmDemod, filt, m_rxBuf[m_rxBufIdx], corr / 100.0, thresholdMet, dcOffset, scopeCRCValid ? 1.0 : (scopeCRCInvalid ? -1.0 : 0));
|
||||
sampleToScope(ci / SDR_RX_SCALEF, magsq, fmDemod, filt, m_rxBuf[m_rxBufIdx],
|
||||
metric,
|
||||
thresholdMet, dcOffset, scopeCRCValid ? 1.0 : (scopeCRCInvalid ? -1.0 : 0));
|
||||
|
||||
// Send demod signal to Demod Analyzer feature
|
||||
m_demodBuffer[m_demodBufferFill++] = fmDemod * std::numeric_limits<int16_t>::max();
|
||||
@@ -373,6 +326,172 @@ void AISDemodSink::processOneSample(Complex &ci)
|
||||
}
|
||||
}
|
||||
|
||||
// HDLC deframing of the symbol decisions computeSymbols() has left in m_symSoft
|
||||
AISDemodSink::DemodResult AISDemodSink::deframe(int& endSampleIdx,
|
||||
bool& crcValid, bool& crcInvalid)
|
||||
{
|
||||
bool gotSOP = false;
|
||||
int bits = 0;
|
||||
int bitCount = 0;
|
||||
int onesCount = 0;
|
||||
int byteCount = 0;
|
||||
int symbolPrev = 0;
|
||||
int symbolIdx = 0;
|
||||
int totalBitCount = 0; // Count of bits after start flag, before bit stuffing removal, including stop flag
|
||||
|
||||
for (int sampleIdx = 0; sampleIdx < m_rxBufLength; sampleIdx += m_samplesPerSymbol)
|
||||
{
|
||||
if (symbolIdx >= (int) m_symSoft.size()) {
|
||||
return gotSOP ? FrameTruncated : FrameNone;
|
||||
}
|
||||
|
||||
// Sign of the Viterbi's decision for this symbol
|
||||
int symbol = m_symSoft[symbolIdx] >= 0.0f ? 1 : 0;
|
||||
|
||||
symbolIdx++;
|
||||
|
||||
// HDLC deframing
|
||||
|
||||
// NRZI decoding
|
||||
int bit;
|
||||
if (symbol != symbolPrev) {
|
||||
bit = 0;
|
||||
} else {
|
||||
bit = 1;
|
||||
}
|
||||
symbolPrev = symbol;
|
||||
|
||||
// Store in shift reg
|
||||
bits |= bit << bitCount;
|
||||
bitCount++;
|
||||
|
||||
if (bit == 1)
|
||||
{
|
||||
onesCount++;
|
||||
// Shouldn't ever get 7 1s in a row
|
||||
if ((onesCount == 7) && gotSOP) {
|
||||
return FrameNone;
|
||||
}
|
||||
}
|
||||
else if (bit == 0)
|
||||
{
|
||||
if (onesCount == 5)
|
||||
{
|
||||
// Remove bit-stuffing (5 1s followed by a 0)
|
||||
bitCount--;
|
||||
}
|
||||
else if (onesCount == 6)
|
||||
{
|
||||
// Start/end of packet
|
||||
if (gotSOP && (bitCount == 8) && (bits == 0x7e) && (byteCount >= 3))
|
||||
{
|
||||
// End of packet
|
||||
// Check CRC is valid
|
||||
m_crc.init();
|
||||
m_crc.calculate(m_bytes, byteCount - 2);
|
||||
uint16_t calcCrc = m_crc.get();
|
||||
uint16_t rxCrc = m_bytes[byteCount-2] | (m_bytes[byteCount-1] << 8);
|
||||
if (calcCrc == rxCrc)
|
||||
{
|
||||
crcValid = true;
|
||||
// Don't include CRC
|
||||
sendMessage(QByteArray((char *) m_bytes, byteCount - 2), totalBitCount);
|
||||
endSampleIdx = sampleIdx;
|
||||
return FrameGood;
|
||||
}
|
||||
else
|
||||
{
|
||||
//qDebug() << QString("CRC mismatch: %1 %2").arg(calcCrc, 4, 16, QLatin1Char('0')).arg(rxCrc, 4, 16, QLatin1Char('0'));
|
||||
crcInvalid = true;
|
||||
return FrameBadCrc;
|
||||
}
|
||||
}
|
||||
else if (gotSOP)
|
||||
{
|
||||
// Repeated start flag without data or misalignment, something not right
|
||||
return FrameNone;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Start of packet
|
||||
gotSOP = true;
|
||||
bits = 0;
|
||||
bitCount = 0;
|
||||
byteCount = 0;
|
||||
totalBitCount = 0;
|
||||
}
|
||||
}
|
||||
onesCount = 0;
|
||||
}
|
||||
|
||||
if (gotSOP)
|
||||
{
|
||||
totalBitCount++;
|
||||
if (bitCount == 8)
|
||||
{
|
||||
// Could also check count according to message ID as that varies
|
||||
if (byteCount >= AISDEMOD_MAX_BYTES)
|
||||
{
|
||||
// Too many bytes
|
||||
return FrameNone;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Got a complete byte
|
||||
m_bytes[byteCount] = bits;
|
||||
byteCount++;
|
||||
}
|
||||
bits = 0;
|
||||
bitCount = 0;
|
||||
}
|
||||
}
|
||||
|
||||
// Abort demod if we haven't found start flag within a couple of bytes of presumed preamble
|
||||
if (!gotSOP && (sampleIdx >= AISDEMOD_SOP_SYMBOLS * m_samplesPerSymbol)) {
|
||||
return FrameNone;
|
||||
}
|
||||
}
|
||||
|
||||
return gotSOP ? FrameTruncated : FrameNone;
|
||||
}
|
||||
|
||||
void AISDemodSink::sendMessage(const QByteArray& rxPacket, int totalBitCount)
|
||||
{
|
||||
if (!getMessageQueueToChannel()) {
|
||||
return;
|
||||
}
|
||||
|
||||
//qDebug() << "RX: " << rxPacket.toHex();
|
||||
|
||||
// Calculate slot number based on time of start of transmission
|
||||
// This is unlikely to be accurate in absolute terms, given we don't know latency from SDR or buffering within SDRangel
|
||||
// But can be used to get an idea of congestion
|
||||
QDateTime currentTime = QDateTime::currentDateTime();
|
||||
if (m_settings.m_useFileTime)
|
||||
{
|
||||
QString hardwareId = m_aisDemod->getDeviceAPI()->getHardwareId();
|
||||
|
||||
if ((hardwareId == "FileInput") || (hardwareId == "SigMFFileInput"))
|
||||
{
|
||||
QString dateTimeStr;
|
||||
int deviceIdx = m_aisDemod->getDeviceSetIndex();
|
||||
|
||||
if (ChannelWebAPIUtils::getDeviceReportValue(deviceIdx, "absoluteTime", dateTimeStr)) {
|
||||
currentTime = QDateTime::fromString(dateTimeStr, Qt::ISODateWithMs);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int txTimeMs = (totalBitCount + 8 + 24 + 8) * (1000.0 / m_settings.m_baud); // Add ramp up, preamble and start-flag
|
||||
QDateTime startDateTime = currentTime.addMSecs(-txTimeMs);
|
||||
int ms = startDateTime.time().second() * 1000 + startDateTime.time().msec();
|
||||
float slotTime = 60.0f * 1000.0f / 2250.0f; // 2250 slots per minute, 26.6ms per slot
|
||||
int slot = ms / slotTime;
|
||||
int totalSlots = std::ceil(txTimeMs / slotTime);
|
||||
AISDemod::MsgMessage *msg = AISDemod::MsgMessage::create(rxPacket, currentTime, slot, totalSlots);
|
||||
getMessageQueueToChannel()->push(msg);
|
||||
}
|
||||
|
||||
void AISDemodSink::applyChannelSettings(int channelSampleRate, int channelFrequencyOffset, bool force)
|
||||
{
|
||||
qDebug() << "AISDemodSink::applyChannelSettings:"
|
||||
@@ -426,6 +545,13 @@ void AISDemodSink::applySettings(const AISDemodSettings& settings, const QString
|
||||
m_rxBufIdx = 0;
|
||||
m_rxBufCnt = 0;
|
||||
|
||||
// The complex baseband the sequence detector works on, kept in step with m_rxBuf
|
||||
m_iqBuf.assign(m_rxBufLength, std::complex<double>(0.0, 0.0));
|
||||
|
||||
m_mlse.create(m_samplesPerSymbol, AISDemodSettings::AISDEMOD_MLSE_SPAN,
|
||||
AISDemodSettings::AISDEMOD_MLSE_BT);
|
||||
m_mlse.setLoopGains(AISDEMOD_MLSE_PHASE_GAIN, AISDEMOD_MLSE_FREQ_GAIN);
|
||||
|
||||
// Create 24-bit training sequence for correlation
|
||||
delete[] m_train;
|
||||
m_correlationLength = 24*m_samplesPerSymbol;
|
||||
@@ -444,6 +570,11 @@ void AISDemodSink::applySettings(const AISDemodSettings& settings, const QString
|
||||
m_train[i*m_samplesPerSymbol+j] = m_pulseShape.filter(trainNRZ[i] * 2.0f - 1.0f);
|
||||
}
|
||||
}
|
||||
|
||||
m_trainEnergy = 0.0f;
|
||||
for (int i = 0; i < m_correlationLength; i++) {
|
||||
m_trainEnergy += m_train[i] * m_train[i];
|
||||
}
|
||||
}
|
||||
|
||||
if (force) {
|
||||
|
||||
@@ -27,6 +27,7 @@
|
||||
#include "dsp/nco.h"
|
||||
#include "dsp/interpolator.h"
|
||||
#include "dsp/gaussian.h"
|
||||
#include "dsp/gmskmlse.h"
|
||||
#include "util/movingaverage.h"
|
||||
#include "util/messagequeue.h"
|
||||
#include "util/crc.h"
|
||||
@@ -39,6 +40,23 @@
|
||||
|
||||
#define AISDEMOD_MAX_BYTES 160
|
||||
|
||||
// Symbols searched for the HDLC start flag before giving up on a correlation trigger
|
||||
#define AISDEMOD_SOP_SYMBOLS 16
|
||||
|
||||
// Symbols the MLSE decodes once a start flag has been found. One AIS slot is about 256
|
||||
// symbols; longer messages occupy several slots and the span doubles until they fit
|
||||
#define AISDEMOD_MLSE_SYMBOLS 320
|
||||
|
||||
// Preamble symbols decoded to pull the per survivor phase loop in, then discarded. The
|
||||
// demodulator starts 18 symbols into the 24 bit training sequence, so there is room for at
|
||||
// most that many - beyond it the look back is clamped rather than reading off the start of
|
||||
// the receive buffer
|
||||
#define AISDEMOD_MLSE_WARMUP 12
|
||||
|
||||
// Per survivor phase tracking loop gains. The optimum is broad
|
||||
#define AISDEMOD_MLSE_PHASE_GAIN 0.3
|
||||
#define AISDEMOD_MLSE_FREQ_GAIN 0.05
|
||||
|
||||
class ChannelAPI;
|
||||
class AISDemod;
|
||||
class ScopeVis;
|
||||
@@ -119,6 +137,12 @@ private:
|
||||
int m_rxBufCnt; // Number of valid samples in buffer
|
||||
Real *m_train; // Training sequence to look for
|
||||
int m_correlationLength;
|
||||
Real m_trainEnergy; // Sum of squares of m_train, for the normalised correlation
|
||||
|
||||
GmskMlse m_mlse; // Coherent sequence detector, used instead of the slicer
|
||||
std::vector<std::complex<double>> m_iqBuf; // Complex baseband, in step with m_rxBuf
|
||||
std::vector<Real> m_mlseSoft; // Viterbi output, including the warm up symbols
|
||||
std::vector<Real> m_symSoft; // Symbol decisions handed to the deframer
|
||||
|
||||
unsigned char m_bytes[AISDEMOD_MAX_BYTES];
|
||||
crc16x25 m_crc;
|
||||
@@ -130,7 +154,19 @@ private:
|
||||
static const int m_sampleBufferSize = AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE / 20;
|
||||
int m_sampleBufferIndex;
|
||||
|
||||
enum DemodResult
|
||||
{
|
||||
FrameNone, //!< Nothing that looked like a frame
|
||||
FrameGood, //!< Complete frame, CRC passed
|
||||
FrameBadCrc, //!< Complete frame, CRC failed
|
||||
FrameTruncated //!< Started a frame but ran out of decoded symbols
|
||||
};
|
||||
|
||||
void processOneSample(Complex &ci);
|
||||
double estimateFrequency() const;
|
||||
void computeSymbols(int x, int n);
|
||||
DemodResult deframe(int& endSampleIdx, bool& crcValid, bool& crcInvalid);
|
||||
void sendMessage(const QByteArray& packet, int totalBitCount);
|
||||
MessageQueue *getMessageQueueToChannel() { return m_messageQueueToChannel; }
|
||||
void sampleToScope(Complex sample, Real magsq, Real fmDemod, Real filt, Real rxBuf, Real corr, Real thresholdMet, Real dcOffset, Real crcValid);
|
||||
};
|
||||
|
||||
@@ -40,11 +40,13 @@ This specifies the bandwidth of a LPF that is applied to the input signal to lim
|
||||
|
||||
<h3>5: Dev - Frequency deviation</h3>
|
||||
|
||||
Adjusts the expected peak frequency deviation in 0.1 kHz steps from 1 to 6 kHz. Typical values are 4.8 kHz, corresponding to a modulation index of 0.5 at 9,600 baud.
|
||||
Adjusts the expected peak frequency deviation in 0.1 kHz steps from 1 to 6 kHz. The default is 2.4 kHz: ITU-R M.1371-5 section 2.3.2 specifies a modulation index of 0.5, and for continuous phase modulation the modulation index is twice the peak deviation divided by the bit rate, so 0.5 at 9,600 bit/s gives a peak deviation of 2,400 Hz. (4,800 Hz is the separation between the mark and space frequencies, which is twice the deviation.)
|
||||
|
||||
This setting scales the output of the FM discriminator, which is used to detect the preamble and drive the scope traces. It has no direct effect on the demodulator's symbol decisions.
|
||||
|
||||
<h3>6: TH - Correlation Threshold</h3>
|
||||
|
||||
The correlation threshold between the received signal and the preamble (training sequence). A lower value should be able to demodulate weaker signals, but increases processor usage and may result in invalid messages if too low.
|
||||
The threshold for the normalised correlation between the received signal and the preamble (training sequence), from 0 to 1. Being normalised, it does not depend on signal level, so the same setting works for strong and weak signals. Real preambles correlate above 0.9 and the default of 0.6 sits well clear of the noise floor. A lower value may demodulate slightly weaker signals, but increases processor usage sharply, because every threshold crossing starts a sequence detection.
|
||||
|
||||
<h3>7: Find</h3>
|
||||
|
||||
|
||||
@@ -413,6 +413,7 @@ set(sdrbase_HEADERS
|
||||
dsp/kissengine.h
|
||||
dsp/firfilter.h
|
||||
dsp/gaussian.h
|
||||
dsp/gmskmlse.h
|
||||
dsp/mimochannel.h
|
||||
dsp/misc.h
|
||||
dsp/movingaverage.h
|
||||
|
||||
@@ -0,0 +1,315 @@
|
||||
///////////////////////////////////////////////////////////////////////////////////
|
||||
// Copyright (C) 2026 Jon Beniston, M7RCE <jon@beniston.com> //
|
||||
// Some code by AI //
|
||||
// //
|
||||
// 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 as 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 V3 for more details. //
|
||||
// //
|
||||
// You should have received a copy of the GNU General Public License //
|
||||
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
||||
///////////////////////////////////////////////////////////////////////////////////
|
||||
|
||||
#ifndef INCLUDE_GMSKMLSE_H
|
||||
#define INCLUDE_GMSKMLSE_H
|
||||
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
#include <vector>
|
||||
|
||||
#include "dsp/dsptypes.h"
|
||||
|
||||
// Maximum likelihood sequence estimator for GMSK with modulation index 1/2.
|
||||
//
|
||||
// GMSK is continuous phase modulation, not frequency shift keying that happens to be
|
||||
// filtered, and demodulating it with a phase discriminator and a slicer gives away
|
||||
// several dB. The transmitted phase is
|
||||
//
|
||||
// phi(t) = pi * sum_i a_i q(t - iT), a_i = +-1
|
||||
//
|
||||
// where q is the integral of the Gaussian frequency pulse, rising from 0 to 1/2 over L
|
||||
// symbol periods. Every symbol smears its phase contribution over L symbols, and only
|
||||
// after L symbols has it contributed its full +-pi/2.
|
||||
//
|
||||
// That makes the signal a finite state machine. Its state is the accumulated phase - four
|
||||
// values, since it advances in steps of pi/2 - together with the L-1 symbols still part
|
||||
// way through their pulse, so 4 * 2^(L-1) states, or 16 for L=3. What a slicer suffers as
|
||||
// inter symbol interference is structure rather than noise, and a sequence detector
|
||||
// resolves it. Unlike a discriminator it also never takes the argument of a noisy sample,
|
||||
// so it has no threshold effect to fall off.
|
||||
//
|
||||
// Branch metrics are formed by correlating the received complex baseband against the exact
|
||||
// waveform each state transition would have produced. The accumulated phase is applied as
|
||||
// a rotation after correlating, so only 2^L correlations are needed per symbol rather than
|
||||
// one per branch.
|
||||
//
|
||||
// Coherent detection needs a phase reference held over the whole burst, which for a signal
|
||||
// with any carrier offset a one shot estimate cannot provide - a few Hz of residual error
|
||||
// walks the phase away over a couple of hundred symbols. So each survivor carries its own
|
||||
// second order decision directed phase tracking loop, updated from its own decisions. A
|
||||
// survivor that is decoding correctly tracks the carrier; one that is not accumulates
|
||||
// phase error and dies, which is what is wanted. Without this the detector performs worse
|
||||
// than a discriminator, so the loop is not optional.
|
||||
//
|
||||
// decode() also returns a per symbol reliability - the metric difference between the best
|
||||
// surviving path asserting +1 and the best asserting -1 - which can drive soft decision
|
||||
// techniques such as Chase decoding of an outer CRC.
|
||||
|
||||
class GmskMlse
|
||||
{
|
||||
public:
|
||||
// samplesPerSymbol - the signal must be sampled at an integer multiple of the baud rate
|
||||
// span - assumed length of the Gaussian phase pulse in symbols, 1 to 4.
|
||||
// 3 is right for the BT values used in practice; 4 costs twice the
|
||||
// CPU for no measurable gain
|
||||
// bt - bandwidth symbol time product of the transmit filter
|
||||
void create(int samplesPerSymbol, int span, double bt)
|
||||
{
|
||||
m_sps = samplesPerSymbol;
|
||||
m_L = std::max(1, std::min(4, span));
|
||||
m_numCombos = 1 << m_L;
|
||||
m_numStates = 4 << (m_L - 1);
|
||||
m_histMask = (1 << (m_L - 1)) - 1;
|
||||
|
||||
// Gaussian frequency pulse, integrated to give the phase pulse q. Computed on a
|
||||
// fine grid and normalised numerically so that q(LT) is exactly 1/2, which is what
|
||||
// makes each symbol's total phase contribution pi/2.
|
||||
const int fine = 64;
|
||||
const int nFine = m_L * m_sps * fine;
|
||||
std::vector<double> g(nFine + 1);
|
||||
|
||||
double k = 2.0 * M_PI * bt / std::sqrt(std::log(2.0));
|
||||
double centre = m_L / 2.0;
|
||||
double sum = 0.0;
|
||||
|
||||
for (int i = 0; i <= nFine; i++)
|
||||
{
|
||||
double t = (double) i / (m_sps * fine) - centre;
|
||||
g[i] = 0.5 * (qfunc(k * (t - 0.5)) - qfunc(k * (t + 0.5)));
|
||||
sum += g[i];
|
||||
}
|
||||
|
||||
double total = sum - 0.5 * (g[0] + g[nFine]);
|
||||
double scale = (total > 0.0) ? (0.5 / total) : 0.0;
|
||||
|
||||
m_q.assign(m_L * m_sps, 0.0);
|
||||
|
||||
double acc = 0.0;
|
||||
|
||||
for (int i = 0, idx = 0; i < nFine; i++)
|
||||
{
|
||||
if ((i % fine) == 0) {
|
||||
m_q[idx++] = acc * scale;
|
||||
}
|
||||
acc += 0.5 * (g[i] + g[i+1]);
|
||||
}
|
||||
|
||||
// Waveform for every combination of the L symbols in flight, over one symbol
|
||||
// period, relative to the accumulated phase state
|
||||
m_w.assign((size_t) m_numCombos * m_sps, std::complex<double>(0.0, 0.0));
|
||||
|
||||
for (int combo = 0; combo < m_numCombos; combo++)
|
||||
{
|
||||
for (int m = 0; m < m_sps; m++)
|
||||
{
|
||||
double phase = 0.0;
|
||||
|
||||
for (int i = 0; i < m_L; i++)
|
||||
{
|
||||
double a = ((combo >> i) & 1) ? 1.0 : -1.0;
|
||||
phase += a * m_q[m + i*m_sps];
|
||||
}
|
||||
|
||||
phase *= M_PI;
|
||||
m_w[(size_t) combo*m_sps + m] = std::complex<double>(std::cos(phase), std::sin(phase));
|
||||
}
|
||||
}
|
||||
|
||||
for (int p = 0; p < 4; p++)
|
||||
{
|
||||
double a = p * M_PI / 2.0;
|
||||
m_rot[p] = std::complex<double>(std::cos(a), std::sin(a));
|
||||
}
|
||||
}
|
||||
|
||||
// Gains for the per survivor phase tracking loop. The optimum is broad; 0.3 and 0.05
|
||||
// measured best for AIS and anything from 0.2/0.01 to 0.6/0.1 is within a few percent.
|
||||
void setLoopGains(double phaseGain, double freqGain)
|
||||
{
|
||||
m_phaseGain = phaseGain;
|
||||
m_freqGain = freqGain;
|
||||
}
|
||||
|
||||
int getStates() const { return m_numStates; }
|
||||
|
||||
// Decode n symbols. samples(k, m) must return sample m of symbol period k as a
|
||||
// std::complex<double>. soft is filled with a signed reliability whose sign is the
|
||||
// decided symbol: positive for +1, negative for -1.
|
||||
//
|
||||
// The loop starts with no knowledge of the carrier phase, so the first several symbols
|
||||
// are unreliable - decode a preamble before the data and discard those symbols.
|
||||
template <class SampleFn>
|
||||
void decode(int n, SampleFn samples, std::vector<Real>& soft)
|
||||
{
|
||||
m_metric.assign(m_numStates, 0.0);
|
||||
m_next.assign(m_numStates, -1e30);
|
||||
m_from.assign((size_t) n * m_numStates, 0);
|
||||
m_delta.assign(n, 0.0);
|
||||
|
||||
m_ph.assign(m_numStates, std::complex<double>(1.0, 0.0));
|
||||
m_phNext.assign(m_numStates, std::complex<double>(1.0, 0.0));
|
||||
m_freq.assign(m_numStates, 0.0);
|
||||
m_freqNext.assign(m_numStates, 0.0);
|
||||
|
||||
soft.assign(n, 0.0f);
|
||||
|
||||
m_corr.resize(m_numCombos);
|
||||
|
||||
for (int kSym = 0; kSym < n; kSym++)
|
||||
{
|
||||
// Correlate this symbol period against every combination of the symbols in
|
||||
// flight. Done once per combination rather than once per branch, since the
|
||||
// accumulated phase state is only a rotation.
|
||||
for (int combo = 0; combo < m_numCombos; combo++)
|
||||
{
|
||||
std::complex<double> c(0.0, 0.0);
|
||||
|
||||
for (int m = 0; m < m_sps; m++) {
|
||||
c += samples(kSym, m) * std::conj(m_w[(size_t) combo*m_sps + m]);
|
||||
}
|
||||
|
||||
m_corr[combo] = c;
|
||||
}
|
||||
|
||||
for (int s = 0; s < m_numStates; s++) {
|
||||
m_next[s] = -1e30;
|
||||
}
|
||||
|
||||
double best1 = -1e30;
|
||||
double best0 = -1e30;
|
||||
|
||||
for (int s = 0; s < m_numStates; s++)
|
||||
{
|
||||
if (m_metric[s] <= -1e29) {
|
||||
continue;
|
||||
}
|
||||
|
||||
int p = s & 3;
|
||||
int hist = s >> 2;
|
||||
std::complex<double> ref = m_rot[p] * m_ph[s];
|
||||
|
||||
for (int u = 0; u < 2; u++)
|
||||
{
|
||||
int combo = u | (hist << 1);
|
||||
std::complex<double> c = m_corr[combo] * std::conj(ref);
|
||||
double metric = m_metric[s] + c.real();
|
||||
|
||||
// The oldest symbol has finished its pulse and joins the accumulated
|
||||
// phase
|
||||
int oldBit = (combo >> (m_L - 1)) & 1;
|
||||
int p2 = (p + (oldBit ? 1 : 3)) & 3;
|
||||
int s2 = p2 | ((combo & m_histMask) << 2);
|
||||
|
||||
if (metric > m_next[s2])
|
||||
{
|
||||
m_next[s2] = metric;
|
||||
m_from[(size_t) kSym * m_numStates + s2] = (s << 1) | u;
|
||||
|
||||
// Second order decision directed phase tracking, carried per
|
||||
// survivor. The error is the sine of the residual branch phase, so
|
||||
// it is bounded and needs no atan, and the rotation update uses a
|
||||
// small angle approximation to avoid a sincos per state.
|
||||
double mag = std::abs(c);
|
||||
double e = (mag > 1e-12) ? (c.imag() / mag) : 0.0;
|
||||
double freq = m_freq[s] + m_freqGain * e;
|
||||
double step = freq + m_phaseGain * e;
|
||||
|
||||
std::complex<double> rot(1.0, step);
|
||||
rot *= 1.0 / std::sqrt(1.0 + step*step);
|
||||
|
||||
m_phNext[s2] = m_ph[s] * rot;
|
||||
m_freqNext[s2] = freq;
|
||||
}
|
||||
|
||||
if (u) {
|
||||
best1 = std::max(best1, metric);
|
||||
} else {
|
||||
best0 = std::max(best0, metric);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
m_delta[kSym] = best1 - best0;
|
||||
|
||||
// Renormalise so the metrics cannot run away over a long burst
|
||||
double best = -1e30;
|
||||
|
||||
for (int s = 0; s < m_numStates; s++) {
|
||||
best = std::max(best, m_next[s]);
|
||||
}
|
||||
for (int s = 0; s < m_numStates; s++) {
|
||||
m_next[s] -= best;
|
||||
}
|
||||
|
||||
m_metric.swap(m_next);
|
||||
m_ph.swap(m_phNext);
|
||||
m_freq.swap(m_freqNext);
|
||||
}
|
||||
|
||||
// Traceback from the best final state
|
||||
int s = 0;
|
||||
double best = -1e30;
|
||||
|
||||
for (int i = 0; i < m_numStates; i++)
|
||||
{
|
||||
if (m_metric[i] > best)
|
||||
{
|
||||
best = m_metric[i];
|
||||
s = i;
|
||||
}
|
||||
}
|
||||
|
||||
for (int kSym = n - 1; kSym >= 0; kSym--)
|
||||
{
|
||||
int f = m_from[(size_t) kSym * m_numStates + s];
|
||||
int u = f & 1;
|
||||
|
||||
soft[kSym] = (Real) ((u ? 1.0 : -1.0) * std::fabs(m_delta[kSym]));
|
||||
s = f >> 1;
|
||||
}
|
||||
}
|
||||
|
||||
private:
|
||||
static double qfunc(double x) { return 0.5 * std::erfc(x / std::sqrt(2.0)); }
|
||||
|
||||
int m_sps = 6;
|
||||
int m_L = 3;
|
||||
int m_numCombos = 8;
|
||||
int m_numStates = 16;
|
||||
int m_histMask = 3;
|
||||
double m_phaseGain = 0.3;
|
||||
double m_freqGain = 0.05;
|
||||
|
||||
std::vector<double> m_q;
|
||||
std::vector<std::complex<double>> m_w;
|
||||
std::complex<double> m_rot[4];
|
||||
|
||||
std::vector<std::complex<double>> m_corr;
|
||||
std::vector<double> m_metric;
|
||||
std::vector<double> m_next;
|
||||
std::vector<int> m_from;
|
||||
std::vector<double> m_delta;
|
||||
std::vector<std::complex<double>> m_ph;
|
||||
std::vector<std::complex<double>> m_phNext;
|
||||
std::vector<double> m_freq;
|
||||
std::vector<double> m_freqNext;
|
||||
};
|
||||
|
||||
#endif // INCLUDE_GMSKMLSE_H
|
||||
@@ -14,6 +14,7 @@ AISDemodSettings:
|
||||
type: number
|
||||
format: float
|
||||
correlationThreshold:
|
||||
description: "Normalised correlation threshold with the preamble, from 0 to 1"
|
||||
type: number
|
||||
format: float
|
||||
udpEnabled:
|
||||
|
||||
@@ -14,6 +14,7 @@ AISDemodSettings:
|
||||
type: number
|
||||
format: float
|
||||
correlationThreshold:
|
||||
description: "Normalised correlation threshold with the preamble, from 0 to 1"
|
||||
type: number
|
||||
format: float
|
||||
udpEnabled:
|
||||
|
||||
Reference in New Issue
Block a user