diff --git a/plugins/channelrx/demodais/aisdemodgui.cpp b/plugins/channelrx/demodais/aisdemodgui.cpp index 420e5d7ac..3e1336e80 100644 --- a/plugins/channelrx/demodais/aisdemodgui.cpp +++ b/plugins/channelrx/demodais/aisdemodgui.cpp @@ -769,7 +769,7 @@ AISDemodGUI::AISDemodGUI(PluginAPI* pluginAPI, DeviceUISet *deviceUISet, Baseban ui->scopeGUI->changeTrigger(0, triggerData); ui->scopeGUI->focusOnTrigger(0); // re-focus to take changes into account in the GUI - m_scopeVis->setLiveRate(9600*6); + m_scopeVis->setLiveRate(AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE); //m_scopeVis->setFreeRun(false); // FIXME: add method rather than call m_scopeVis->configure() ui->deltaFrequencyLabel->setText(QString("%1f").arg(QChar(0x94, 0x03))); diff --git a/plugins/channelrx/demodais/aisdemodsettings.cpp b/plugins/channelrx/demodais/aisdemodsettings.cpp index 11d27cead..00342423c 100644 --- a/plugins/channelrx/demodais/aisdemodsettings.cpp +++ b/plugins/channelrx/demodais/aisdemodsettings.cpp @@ -38,7 +38,7 @@ void AISDemodSettings::resetToDefaults() m_inputFrequencyOffset = 0; m_rfBandwidth = 16000.0f; m_fmDeviation = 2400.0f; - m_correlationThreshold = 0.6f; + m_correlationThreshold = 0.7f; m_filterMMSI = ""; m_udpEnabled = false; m_udpAddress = "127.0.0.1"; @@ -107,7 +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); + s.writeFloat(37, m_correlationThreshold); for (int i = 0; i < AISDEMOD_MESSAGE_COLUMNS; i++) s.writeS32(100 + i, m_messageColumnIndexes[i]); @@ -141,7 +141,7 @@ bool AISDemodSettings::deserialize(const QByteArray& data) d.readString(7, &m_udpAddress, "127.0.0.1"); d.readU32(8, &utmp); - if ((utmp > 1023) && (utmp < 65535)) { + if ((utmp > 1023) && (utmp <= 65535)) { m_udpPort = utmp; } else { m_udpPort = 9999; @@ -195,7 +195,7 @@ 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); + d.readFloat(37, &m_correlationThreshold, 0.7f); for (int i = 0; i < AISDEMOD_MESSAGE_COLUMNS; i++) { d.readS32(100 + i, &m_messageColumnIndexes[i], i); diff --git a/plugins/channelrx/demodais/aisdemodsettings.h b/plugins/channelrx/demodais/aisdemodsettings.h index b6d477a03..dfaf4e8a8 100644 --- a/plugins/channelrx/demodais/aisdemodsettings.h +++ b/plugins/channelrx/demodais/aisdemodsettings.h @@ -40,7 +40,10 @@ struct AISDemodSettings //!< 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 + Real m_correlationThreshold; //!< Normalised matched filter output for the preamble, + //!< 0 to 1. Noise sits near 0.886/sqrt(block length), + //!< so this cannot be retuned without regard to + //!< AISDEMOD_IQ_BLOCKS - the two are coupled QString m_filterMMSI; bool m_udpEnabled; QString m_udpAddress; @@ -74,15 +77,39 @@ struct AISDemodSettings int m_messageColumnSizes[AISDEMOD_MESSAGE_COLUMNS]; //!< Size of the columns in the table static const int AISDEMOD_BAUD_RATE = 9600; - 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 AISDEMOD_CHANNEL_SAMPLE_RATE = 96000; //!< 10x 9600 baud rate (use even multiple so Gaussian filter has odd number of taps). + //!< 6x was enough for a discriminator, but the sequence detector is sensitive to + //!< symbol timing and the finer grid is worth measurable sensitivity 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; + + //! Preamble detection runs a matched filter on the complex baseband rather than + //! correlating the discriminator output, which has an FM threshold and so gives up + //! exactly where the sequence detector still works. Measured at the noise floor it + //! reaches the same message count on 629 triggers where the old statistic needed + //! 173,000. + //! + //! Carrier phase is unknown so the blocks are combined non-coherently, and carrier + //! frequency is unknown too - a single coherent correlation over all 24 symbols would + //! be cancelled by a few hundred Hz. 6 blocks tolerates about +/-1200 Hz; 4 was tried + //! and lost real messages on a recording with larger offsets. + static const int AISDEMOD_IQ_BLOCKS = 6; + + //! The correlation peak is about a symbol wide, so it does not need evaluating on + //! every one of 10 samples per symbol. Sub sample alignment is recovered by the + //! timing phase retries, which search +/-2 samples anyway. + static const int AISDEMOD_IQ_DECIM = 2; + + static const int AISDEMOD_MLSE_SPAN = 4; //!< Symbols of Gaussian pulse the trellis resolves. 4 gives a 32 state + //!< trellis against 16 for 3, and is the largest single sensitivity + //!< win available - measured +25% messages on a weak recording. An + //!< earlier note here said 3 was enough and 4 measured no better; + //!< that was at 6 samples per symbol with the symbol timing a third + //!< of a symbol out, so the trellis could not use the extra span AISDemodSettings(); void resetToDefaults(); diff --git a/plugins/channelrx/demodais/aisdemodsink.cpp b/plugins/channelrx/demodais/aisdemodsink.cpp index cc965dc6b..597c8e293 100644 --- a/plugins/channelrx/demodais/aisdemodsink.cpp +++ b/plugins/channelrx/demodais/aisdemodsink.cpp @@ -33,6 +33,7 @@ AISDemodSink::AISDemodSink(AISDemod *aisDemod) : m_scopeSink(nullptr), m_aisDemod(aisDemod), + m_channel(nullptr), m_channelSampleRate(AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE), m_channelFrequencyOffset(0), m_magsqSum(0.0f), @@ -45,6 +46,9 @@ AISDemodSink::AISDemodSink(AISDemod *aisDemod) : m_sampleBufferIndex(0) { m_magsq = 0.0; + m_sampleCounter = 0; + m_lastAttemptPos = 0; + m_haveAttempted = false; m_demodBuffer.resize(1<<12); m_demodBufferFill = 0; @@ -206,6 +210,7 @@ void AISDemodSink::processOneSample(Complex &ci) m_rxBuf[m_rxBufIdx] = filtClipped; m_iqBuf[m_rxBufIdx] = std::complex(ci.real() / SDR_RX_SCALEF, ci.imag() / SDR_RX_SCALEF); m_rxBufIdx = (m_rxBufIdx + 1) % m_rxBufLength; + m_sampleCounter++; m_rxBufCnt = std::min(m_rxBufCnt + 1, m_rxBufLength); Real corr = 0.0f; @@ -217,28 +222,41 @@ 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 - // Calculate sum to estimate DC offset - for (int i = 0; i < m_correlationLength; i++) + // Matched filter for the preamble, on the complex baseband. Correlating the + // discriminator output instead - as this used to - means detecting with a statistic + // that has already fallen off the FM threshold at the levels the sequence detector + // can still demodulate, which is why noise used to trigger it for several percent of + // all samples. Blocks are combined non-coherently because neither carrier phase nor + // frequency is known; see AISDEMOD_IQ_BLOCKS. + if ((m_sampleCounter % AISDemodSettings::AISDEMOD_IQ_DECIM) == 0) { - 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]; - } + const int blocks = AISDemodSettings::AISDEMOD_IQ_BLOCKS; + const int blockLen = m_correlationLength / blocks; + double sumMag = 0.0; + double iqEnergy = 0.0; - // If we meet threshold, try to demod - // 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; + for (int b = 0; b < blocks; b++) + { + std::complex acc(0.0, 0.0); + + for (int i = b*blockLen; i < (b+1)*blockLen; i++) + { + int j = (m_rxBufIdx + i) % m_rxBufLength; + acc += m_iqBuf[j] * std::conj(m_trainIQ[i]); + iqEnergy += std::norm(m_iqBuf[j]); + trainingSum += m_rxBuf[j]; + } + + sumMag += std::abs(acc); + } + + // |m_trainIQ| is 1, so a clean match of amplitude A gives sumMag = A*N against + // sqrt(iqEnergy*N) = A*N, i.e. 1. Noise lands near 0.886/sqrt(blockLen). + corr = (Real) sumMag; + metric = (Real) (sumMag / (sqrt(iqEnergy * (double) (blockLen*blocks)) + 1e-12)); + thresholdMet = metric >= m_settings.m_correlationThreshold; + } if (thresholdMet) { @@ -247,7 +265,46 @@ void AISDemodSink::processOneSample(Complex &ci) dcOffset = trainingSum/m_correlationLength; // Start demod after (most of) preamble - int x = (m_rxBufIdx + m_correlationLength*3/4 + 4) % m_rxBufLength; + // Symbol centres sit at m_samplesPerSymbol/2 - 1 past the training symbol + // boundary: the transmit pulse shaping delay and the receive one cancel, + // leaving the half symbol offset and the one sample the phase discriminator + // consumes. Do not write this as a constant - it was +4, which is only + // correct at 10 samples per symbol, and the MLSE will not tolerate a quarter + // symbol error the way the old slicer did. + int base = (m_rxBufIdx + m_correlationLength*3/4 + + m_samplesPerSymbol/2 - 1) % m_rxBufLength; + + // Try the nominal timing first, then progressively larger offsets either side + for (int phase = 0; phase <= 2*AISDEMOD_TIMING_PHASES; phase++) + { + // Only the alignment that is kept should mark the scope trace. deframe() + // only ever sets these, so a bad CRC from an alignment that is about to be + // retried would otherwise stay set even when a later phase succeeds. + scopeCRCValid = false; + scopeCRCInvalid = false; + + int d = (phase + 1) / 2; + + if (phase & 1) { + d = -d; + } + + // Consecutive triggers overlap by all but one alignment, and re-demodulating a + // position almost always gives the same answer, so skipping the repeats is worth + // about 2.5x the CPU of the whole retry stage. Not quite free: the look ahead + // available from a position grows as samples arrive, so a multi slot frame that + // ran out of buffer on the first attempt can succeed on a later one. Measured at + // one message in 1550 on a clean recording and none at the noise floor. + qint64 pos = (qint64) m_sampleCounter + d; + + if (m_haveAttempted && (pos <= m_lastAttemptPos)) { + continue; + } + + m_lastAttemptPos = pos; + m_haveAttempted = true; + + int x = (base + d + m_rxBufLength) % m_rxBufLength; int endSampleIdx = 0; @@ -291,6 +348,8 @@ void AISDemodSink::processOneSample(Complex &ci) { // Skip over received packet, so we don't try to re-demodulate it m_rxBufCnt -= endSampleIdx; + break; + } } } } @@ -476,8 +535,18 @@ void AISDemodSink::sendMessage(const QByteArray& rxPacket, int totalBitCount) QString dateTimeStr; int deviceIdx = m_aisDemod->getDeviceSetIndex(); - if (ChannelWebAPIUtils::getDeviceReportValue(deviceIdx, "absoluteTime", dateTimeStr)) { - currentTime = QDateTime::fromString(dateTimeStr, Qt::ISODateWithMs); + if (ChannelWebAPIUtils::getDeviceReportValue(deviceIdx, "absoluteTime", dateTimeStr)) + { + QDateTime fileTime = QDateTime::fromString(dateTimeStr, Qt::ISODateWithMs); + + // An unparseable timestamp gives an invalid QDateTime, whose time().second() + // is -1. That makes ms and then the slot number negative, which the slot map + // indexes with and the log reports. Fall back to the wall clock instead. + if (fileTime.isValid()) { + currentTime = fileTime; + } else { + qDebug() << "AISDemodSink::sendMessage: could not parse absoluteTime" << dateTimeStr; + } } } } @@ -486,6 +555,11 @@ void AISDemodSink::sendMessage(const QByteArray& rxPacket, int totalBitCount) 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 + + if (ms < 0) { + ms = 0; + } + int slot = ms / slotTime; int totalSlots = std::ceil(txTimeMs / slotTime); AISDemod::MsgMessage *msg = AISDemod::MsgMessage::create(rxPacket, currentTime, slot, totalSlots); @@ -534,7 +608,16 @@ void AISDemodSink::applySettings(const AISDemodSettings& settings, const QString if ((settingsKeys.contains("baud")) || force) { - m_samplesPerSymbol = AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE / settings.m_baud; + // Clamp before dividing: baud is written by the web API and read from the config + // blob without validation, and 0 divides while anything over the channel rate leaves + // m_samplesPerSymbol at 0, which makes m_rxBufLength 0 and the next sample a % 0. + int baud = settings.m_baud; + + if ((baud < 1) || (baud > AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE)) { + baud = 9600; + } + + m_samplesPerSymbol = AISDemodSettings::AISDEMOD_CHANNEL_SAMPLE_RATE / baud; qDebug() << "AISDemodSink::applySettings: m_samplesPerSymbol: " << m_samplesPerSymbol << " baud " << settings.m_baud; m_pulseShape.create(0.5, 3, m_samplesPerSymbol); @@ -571,6 +654,22 @@ void AISDemodSink::applySettings(const AISDemodSettings& settings, const QString } } + // The transmitted preamble, for matched filtering on the complex baseband. + // m_train is the Gaussian filtered NRZ, i.e. the expected instantaneous + // frequency; integrating gives phase, and GMSK with h=1/2 puts +-pi/2 in each + // symbol while m_train sums to +-samplesPerSymbol over one. + m_trainIQ.assign(m_correlationLength, std::complex(0.0, 0.0)); + { + double phase = 0.0; + double k = M_PI / (2.0 * m_samplesPerSymbol); + + for (int i = 0; i < m_correlationLength; i++) + { + phase += k * m_train[i]; + m_trainIQ[i] = std::complex(std::cos(phase), std::sin(phase)); + } + } + m_trainEnergy = 0.0f; for (int i = 0; i < m_correlationLength; i++) { m_trainEnergy += m_train[i] * m_train[i]; diff --git a/plugins/channelrx/demodais/aisdemodsink.h b/plugins/channelrx/demodais/aisdemodsink.h index 79abdef68..bb2a90bb6 100644 --- a/plugins/channelrx/demodais/aisdemodsink.h +++ b/plugins/channelrx/demodais/aisdemodsink.h @@ -53,6 +53,12 @@ // the receive buffer #define AISDEMOD_MLSE_WARMUP 12 +// Symbol timing offsets tried either side of the nominal start, in samples. The correlator +// locates the burst to about a sample, but the sequence detector wants better than that, and +// one sample is a tenth of a symbol at 96 kHz. Retries only happen when the nominal timing +// fails to produce a good frame, so the cost falls on bursts that would otherwise be lost. +#define AISDEMOD_TIMING_PHASES 2 + // Per survivor phase tracking loop gains. The optimum is broad #define AISDEMOD_MLSE_PHASE_GAIN 0.3 #define AISDEMOD_MLSE_FREQ_GAIN 0.05 @@ -112,7 +118,10 @@ private: ChannelAPI *m_channel; int m_channelSampleRate; int m_channelFrequencyOffset; - int m_samplesPerSymbol; // Number of samples per symbol + int m_samplesPerSymbol; + quint64 m_sampleCounter; //!< Monotonic, so already tried alignments can be skipped + qint64 m_lastAttemptPos; + bool m_haveAttempted; // Number of samples per symbol NCO m_nco; Interpolator m_interpolator; @@ -135,7 +144,8 @@ private: int m_rxBufLength; // Size in elements in m_rxBuf int m_rxBufIdx; // Index in to circular buffer int m_rxBufCnt; // Number of valid samples in buffer - Real *m_train; // Training sequence to look for + Real *m_train; + std::vector> m_trainIQ; //!< Preamble as transmitted, for matched filtering // Training sequence to look for int m_correlationLength; Real m_trainEnergy; // Sum of squares of m_train, for the normalised correlation diff --git a/plugins/channelrx/demodais/readme.md b/plugins/channelrx/demodais/readme.md index d3ab7b662..2be282416 100644 --- a/plugins/channelrx/demodais/readme.md +++ b/plugins/channelrx/demodais/readme.md @@ -8,7 +8,7 @@ AIS is broadcast globally on 25kHz channels at 161.975MHz and 162.025MHz, with o The AIS demodulators can send received messages to the [AIS feature](../../feature/ais/readme.md), which displays a table combining the latest data for vessels amalgamated from multiple demodulators and sends their positions to the [Map Feature](../../feature/map/readme.md) for display in 2D or 3D. -AIS uses GMSK/FM modulation at a baud rate of 9,600, with a modulation index of 0.5. The demodulator works at a sample rate of 57,600Sa/s. +AIS uses GMSK/FM modulation at a baud rate of 9,600, with a modulation index of 0.5. The demodulator works at a sample rate of 96,000Sa/s. Received AIS messages can be NMEA encoded and forwarded via UDP to 3rd party applications. diff --git a/sdrbase/dsp/gmskmlse.h b/sdrbase/dsp/gmskmlse.h index 55403df97..853db9b26 100644 --- a/sdrbase/dsp/gmskmlse.h +++ b/sdrbase/dsp/gmskmlse.h @@ -67,8 +67,12 @@ 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 + // 4 doubles the state count over 3 and is worth it: measured 25% more + // AIS messages on a weak recording, and a fifth of a dB. An earlier + // measurement said 3 was enough, but that was taken at 6 samples per + // symbol with the symbol timing a third of a symbol out - the trellis + // could not use the extra span it was being given. Worth re-measuring + // rather than assuming either way, and cheap to try // bt - bandwidth symbol time product of the transmit filter void create(int samplesPerSymbol, int span, double bt) {