Improve AIS acquisition.

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)));

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);

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();

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<double>(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<double> 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<double>(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<double>(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];

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<std::complex<double>> 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
 

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.
 

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)
     {