wip integrate SNR est into RADE V1 tools

radae/dsp.py

@@ -372,7 +372,9 @@ def __init__(self,latent_dim,Fs,M,Ncp,Wfwd,Nc,Ns,w,P,bottleneck,pilot_gain,time_
          a = local_path_delay_s*self.Fs
          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
+
    def est_pilots(self, rx_sym_pilots, num_modem_frames, Nc, Ns):
       rx_pilots = torch.zeros(num_modem_frames+1, Nc, dtype=torch.complex64)
       # 3-pilot least squares fit across frequency, ref: freedv_low.pdf
@@ -391,15 +393,31 @@ def est_pilots(self, rx_sym_pilots, num_modem_frames, Nc, Ns):
                rx_pilots[i,c] = g[0] + g[1]*torch.exp(-1j*self.w[c]*a)
 
       return rx_pilots
-
+
+   # update SNR estimate
+   def update_snr_est(self, rx_sym_pilots, rx_pilots):
+      Pcn_hat = rx_sym_pilots[0,0,0,:]
+      rx_phase = torch.angle(rx_pilots[0,:])
+      Rcn_hat = Pcn_hat * torch.exp(-1j*rx_phase)
+      S1 = torch.sum(torch.abs(Pcn_hat)**2)
+      S2 = torch.sum(torch.abs(Rcn_hat.imag)**2)   
+      snr_est = S1/(2*S2) - 1
+      # remove occasional illegal values
+      if snr_est <= 0:
+         snr_est = 0.1
+      snrdB_est = 10*np.log10(snr_est)
+      # moving average smoothing, roughly 1 second time constant
+      self.snrdB_est = 0.9*self.snrdB_est + 0.1*snrdB_est
+
    # One frame version of do_pilot_eq() for streaming implementation
    def do_pilot_eq_one(self, num_modem_frames, rx_sym_pilots):
       Nc = self.Nc 
       Ns = self.Ns + 1
 
-      # First, estimate the (complex) value of each received pilot symbol
+      # First, estimate the (complex) value of each received pilot symbol, and update SNR est
       rx_pilots = self.est_pilots(rx_sym_pilots, num_modem_frames, Nc, Ns)
-
+      self.update_snr_est(rx_sym_pilots, rx_pilots)
+
       # Linearly interpolate between two pilots to EQ data symbol phase
       for i in torch.arange(num_modem_frames):
          for c in torch.arange(0,Nc):

radae_rxe.py

@@ -230,6 +230,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)
             if auxdata and self.state == "sync":
                print(f" uw_err: {self.uw_errors:d}", file=sys.stderr)
             else: