first pass C API integration of SNR estimtor

doc/radae.tex

@@ -979,16 +979,26 @@ \section{SNR Estimator Attempt 2}
 
 This algorithm was tested in the script \emph{est\_snr.py}, which supports single point testing and generating curves over evolving multiplath channels. There was some error in the $S/N$ estimates (comnpared to a genie estimator) due to (\ref{eq:snr2_noise}) being sensitive to the noise in the phase estimator $\mathcal{N}(0,\alpha)$.  The algorithm works reasonably well for AWGN, MPG and MPP channels, however the error increases with faster fading.  A straight line was fitted to the estimator output:
 \begin{equation}
-SNR_{est} = m*SNR + c
+SNR_{est} = mSNR + c
 \end{equation}
 where $SNR$ is the genie estimate in dB for the time window $T$. To correct for the error the equation can be re-arranged:
 \begin{equation}
 SNR = (SNR_{est} - c)/m
 \end{equation}
 The straight line fit is different for each channel tested.  So the mean of $m$ and $c$ for each channel was found and averaged.
 
-Figure \ref{fig:est_snr2a} shows the results for an AWGN channel, and Figure \ref{fig:est_snr2b} shows a straight line fit for the estimator over the three channels tested (in both figures the correction has been applied). The results were calculated over 60 second runs using the \emph{est\_snr\_curves.sh} script. As the algorithm works at the symbol rate $S/N=E_s/N_0$.  To obtain the SNR in a different noise bandwidth (e.g 3000Hz) a constant offset will need to be applied. Between -5 and 10dB, the estimator is within +/-1dB for all channels except MPP (1 Hz fading) which reads several dB low.
- 
+Figure \ref{fig:est_snr2a} shows the results for an AWGN channel, and Figure \ref{fig:est_snr2b} shows a straight line fit for the estimator over the three channels tested (in both figures the correction has been applied). The results were calculated over 60 second runs using the \emph{est\_snr\_curves.sh} script.
+
+As the algorithm works at the symbol rate $S/N=E_s/N_0$.  Between -5 and 10dB $E_s/N_0$, the estimator is within +/-1dB for all channels except MPP (1 Hz fading) which reads several dB low. To obtain the SNR in a different noise bandwidth (e.g 3000Hz), and using the RADE V1 waveform parameters:
+\begin{equation}
+\begin{split}
+SNR_{3k} &= 10log_{10}(E_s/N_0) + 10log_{10}(RsNc/B) + 10log_{10}((M+N_{cp})/M) \\
+         &= 10log_{10}(E_s/N_0) + 10log_{10}((50)(30)/3000) + 10log_{10}((160+32)/160) \\
+         &= 10log_{10}(E_s/N_0) - 2.22
+\end{split}
+\end{equation}
+The last term accounts for the carrier power in the cyclic prefix, which is not included in $Es$.
+
 \begin{figure}[H]
 \caption{Straight line fit to estimator output for AWGN, MPG, and MPP channels.}
 \label{fig:est_snr2b}

radae/dsp.py

@@ -373,7 +373,7 @@ def __init__(self,latent_dim,Fs,M,Ncp,Wfwd,Nc,Ns,w,P,bottleneck,pilot_gain,time_
          A = torch.tensor([[1, torch.exp(-1j*self.w[c_mid-1]*a)], [1, torch.exp(-1j*self.w[c_mid]*a)], [1, torch.exp(-1j*self.w[c_mid+1]*a)]])
          self.Pmat[c] = torch.matmul(torch.inverse(torch.matmul(torch.transpose(A,0,1),A)),torch.transpose(A,0,1))
 
-      self.snrdB_est = 0
+      self.snrdB_3k_est = 0
       self.m = 0.8070
       self.c = 2.513
 
@@ -409,10 +409,13 @@ def update_snr_est(self, rx_sym_pilots, rx_pilots):
          snr_est = 0.1
       snrdB_est = 10*np.log10(snr_est)
       # correction based on average of straight line fit to AWGN/MPG/MPP
-      snrdB_est = (snrdB_est - self.c)/self.match
-
+      snrdB_est = (snrdB_est - self.c)/self.m
+      # convert to 3000Hz noise badnwidth, and account for carrier power in cyclic prefix
+      Rs = self.Fs/self.M
+      snrdB_3k_est = snrdB_est + 10*math.log10(Rs*self.Nc/3000) + 10*math.log10((self.M+self.Ncp)/self.M)
+
       # moving average smoothing, roughly 1 second time constant
-      self.snrdB_est = 0.9*self.snrdB_est + 0.1*snrdB_est
+      self.snrdB_3k_est = 0.9*self.snrdB_3k_est + 0.1*snrdB_3k_est
 
    # One frame version of do_pilot_eq() for streaming implementation
    def do_pilot_eq_one(self, num_modem_frames, rx_sym_pilots):

radae_rxe.py

@@ -154,6 +154,9 @@ def get_nin(self):
    def get_sync(self):
       return self.state == "sync"
 
+   def get_snrdB_3k_est(self):
+      return int(self.receiver.snrdB_3k_est)
+
    def sum_uw_errors(self,new_uw_errors):
       self.uw_errors += new_uw_errors
 
@@ -230,7 +233,7 @@ def do_radae_rx(self, buffer_complex, floats_out):
          if v == 2 or (v == 1 and (self.state == "search" or self.state == "candidate" or prev_state == "candidate")):
             print(f"{self.mf:3d} state: {self.state:10s} valid: {candidate:d} {endofover:d} {self.valid_count:2d} Dthresh: {acq.Dthresh:8.2f} ", end='', file=sys.stderr)
             print(f"Dtmax12: {acq.Dtmax12:8.2f} {acq.Dtmax12_eoo:8.2f} tmax: {self.tmax:4d} fmax: {self.fmax:6.2f}", end='', file=sys.stderr)
-            print(f" SNRdB: {receiver.snrdB_est:5.2f}", end='', file=sys.stderr)
+            print(f" SNRdB: {receiver.snrdB_3k_est:5.2f}", end='', file=sys.stderr)
             if auxdata and self.state == "sync":
                print(f" uw_err: {self.uw_errors:d}", file=sys.stderr)
             else:

src/radae_rx.c

@@ -42,6 +42,7 @@ int main(int argc, char *argv[])
             fflush(stdout);
         }
         nin = rade_nin(r);
+        fprintf(stderr, "SNR3k (dB): %d\n", rade_snrdB_3k_est(r));
     }
 
     rade_close(r);

src/rade_api.c

@@ -514,3 +514,16 @@ float rade_freq_offset(struct rade *r) {
   return 0;
 }
 
+RADE_EXPORT int rade_snrdB_3k_est(struct rade *r) {
+  assert(r != NULL);
+
+  // Acquire the Python GIL (needed for multithreaded use)
+  PyGILState_STATE gstate = PyGILState_Ensure();
+
+  int result = (int)call_getter(r->pInst_radae_rx, "get_snrdB_3k_est");
+
+  // Release Python GIL
+  PyGILState_Release(gstate);
+
+  return result;
+}

src/rade_api.h

@@ -118,6 +118,9 @@ RADE_EXPORT int rade_sync(struct rade *r);
 // returns the current frequency offset of the Rx signal ( when rade_sync()!=0 )
 RADE_EXPORT float rade_freq_offset(struct rade *r);
 
+// returns the current SNR estimate (in dB) of the Rx signal ( when rade_sync()!=0 )
+RADE_EXPORT int rade_snrdB_3k_est(struct rade *r);
+
 #ifdef __cplusplus
 }
 #endif