From 76c958b03ef8a6d541a3ab983d292edd95ee6710 Mon Sep 17 00:00:00 2001
From: Oliver Jowett <oliver@mutability.co.uk>
Date: Mon, 22 Sep 2014 14:53:06 +0100
Subject: [PATCH] Use a lookup table for SNR calculation.

---
 dump1090.c |  8 ++++++++
 dump1090.h |  3 ++-
 mode_ac.c  | 12 +++++++++---
 mode_s.c   | 19 +++++++++++--------
 4 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/dump1090.c b/dump1090.c
index b3bef28..3db8ee5 100644
--- a/dump1090.c
+++ b/dump1090.c
@@ -98,6 +98,7 @@ void modesInit(void) {
          ((Modes.pFileData  = (uint16_t *) malloc(MODES_ASYNC_BUF_SIZE)                                         ) == NULL) ||
          ((Modes.magnitude  = (uint16_t *) malloc(MODES_ASYNC_BUF_SIZE+MODES_PREAMBLE_SIZE+MODES_LONG_MSG_SIZE) ) == NULL) ||
          ((Modes.maglut     = (uint16_t *) malloc(sizeof(uint16_t) * 256 * 256)                                 ) == NULL) ||
+         ((Modes.log10lut   = (uint16_t *) malloc(sizeof(uint16_t) * 256 * 256)                                 ) == NULL) ||
          ((Modes.beastOut   = (char     *) malloc(MODES_RAWOUT_BUF_SIZE)                                        ) == NULL) ||
          ((Modes.rawOut     = (char     *) malloc(MODES_RAWOUT_BUF_SIZE)                                        ) == NULL) ) 
     {
@@ -185,6 +186,13 @@ void modesInit(void) {
         }
     }
 
+    // Prepare the log10 lookup table.
+    // This maps from a magnitude value x (scaled as above) to 100log10(x)
+    for (i = 0; i <= 65535; i++) {
+        int l10 = (int) round(100 * log10( (i + 365.4798) * 258.433254) );
+        Modes.log10lut[i] = (uint16_t) ((l10 < 65535 ? l10 : 65535));
+    }
+
     // Prepare error correction tables
     modesInitErrorInfo();
 }
diff --git a/dump1090.h b/dump1090.h
index fe607d8..bcb3e80 100644
--- a/dump1090.h
+++ b/dump1090.h
@@ -88,7 +88,7 @@
 #define MODES_ASYNC_BUF_SAMPLES    (MODES_ASYNC_BUF_SIZE / 2) // Each sample is 2 bytes
 #define MODES_AUTO_GAIN            -100                       // Use automatic gain
 #define MODES_MAX_GAIN             999999                     // Use max available gain
-#define MODES_MSG_SQUELCH_FACTOR   16                         // Min SNR expressed as an amplitude ratio, scaled by 10. 20log(16/10) = 4.1dB
+#define MODES_MSG_SQUELCH_DB       4.0                        // Minimum SNR, in dB
 #define MODES_MSG_ENCODER_ERRS     3                          // Maximum number of encoding errors
 
 // When changing, change also fixBitErrors() and modesInitErrorTable() !!
@@ -257,6 +257,7 @@ struct {                             // Internal state
     int             fd;              // --ifile option file descriptor
     uint32_t       *icao_cache;      // Recently seen ICAO addresses cache
     uint16_t       *maglut;          // I/Q -> Magnitude lookup table
+    uint16_t       *log10lut;        // Magnitude -> log10 lookup table
     int             exit;            // Exit from the main loop when true
 
     // RTLSDR
diff --git a/mode_ac.c b/mode_ac.c
index d3d8d47..d96b7b3 100644
--- a/mode_ac.c
+++ b/mode_ac.c
@@ -144,7 +144,7 @@ int detectModeA(uint16_t *m, struct modesMessage *mm)
   int F1_sig, F1_noise;
   int F2_sig, F2_noise;
   int fSig, fNoise, fLevel, fLoLo;
-  float snr;
+  int snr;
 
   // m[0] contains the energy from    0 ->  499 nS
   // m[1] contains the energy from  500 ->  999 nS
@@ -307,8 +307,14 @@ int detectModeA(uint16_t *m, struct modesMessage *mm)
   if ((ModeABits < 3) || (ModeABits & 0xFFFF8808) || (ModeAErrs) )
     {return (ModeABits = 0);}
 
-  snr = 5 * 20.0 * log10f((float)(fSig + fNoise + 365) / (fNoise + 365)); // 365 to adjust for magnitude value offset
-  mm->signalLevel = ((snr < 255) ? (uint8_t)round(snr) : 255);
+  // snr = 5 * 20log10(fSig / (fSig+fNoise))         (in units of 0.2dB)
+  //     = 100log10(fSig) - 100log10(fSig+fNoise)
+  while (fSig > 65535 || (fSig + fNoise) > 65535) {
+      fSig >>= 1;
+      fNoise >>= 1;
+  }
+  snr = Modes.log10lut[fSig] - Modes.log10lut[fSig + fNoise];
+  mm->signalLevel = ((snr < 255) ? (uint8_t)snr : 255);
 
   return ModeABits;
   }
diff --git a/mode_s.c b/mode_s.c
index e730003..da09f52 100644
--- a/mode_s.c
+++ b/mode_s.c
@@ -1575,6 +1575,7 @@ void detectModeS(uint16_t *m, uint32_t mlen) {
         uint8_t  theByte, theErrs;
         int msglen, scanlen;
         uint32_t sigLevel, noiseLevel;
+        uint16_t snr;
 
         pPreamble = &m[j];
         pPayload  = &m[j+MODES_PREAMBLE_SAMPLES];
@@ -1779,22 +1780,24 @@ void detectModeS(uint16_t *m, uint32_t mlen) {
             }
         }
 
-        // adjust for magnitude zero offset
-        sigLevel += 365*56;
-        noiseLevel += 365*56;
+        // snr = 5 * 20log10(sigLevel / noiseLevel)         (in units of 0.2dB)
+        //     = 100log10(sigLevel) - 100log10(noiseLevel)
+
+        while (sigLevel > 65535 || noiseLevel > 65535) {
+            sigLevel >>= 1;
+            noiseLevel >>= 1;
+        }
+        snr = Modes.log10lut[sigLevel] - Modes.log10lut[noiseLevel];
 
         // When we reach this point, if error is small, and the signal strength is large enough
         // we may have a Mode S message on our hands. It may still be broken and the CRC may not 
         // be correct, but this can be handled by the next layer.
         if ( (msglen) 
-          && ((sigLevel * 10) > (noiseLevel * MODES_MSG_SQUELCH_FACTOR))    // (sigLevel/noiseLevel) > (MODES_MSG_SQUELCH_FACTOR/10)
+          && ((2 * snr) > (int) (MODES_MSG_SQUELCH_DB * 10))
           && (errors      <= MODES_MSG_ENCODER_ERRS) ) {
-            float snr;
-
             // Set initial mm structure details
             mm.timestampMsg = Modes.timestampBlk + (j*6);
-            snr = 5.0 * 20.0 * log10f( (float)sigLevel / noiseLevel ); // sig/noise levels are amplitudes, so square them when computing SNR
-            mm.signalLevel = (snr > 255 ? 255 : (uint8_t)round(snr));
+            mm.signalLevel = (snr > 255 ? 255 : (uint8_t)snr);
             mm.phase_corrected = use_correction;
 
             // Decode the received message