From a0d471cb0bd5b12a36c777d8eb0f49e6d35a906f Mon Sep 17 00:00:00 2001
From: Steve Franke <s.j.franke@icloud.com>
Date: Wed, 27 Nov 2019 15:58:52 -0600
Subject: [PATCH] Improve FT8 SNR estimates in two ways: (i) SNR no longer
 saturates at +20 dB (ii) a large signal in the passband no longer causes the
 SNR of weaker signals to be biased low.

---
 CMakeLists.txt                    |  1 +
 lib/ft8/baseline.f90              |  1 +
 lib/ft8/get_spectrum_baseline.f90 | 43 +++++++++++++++++++++++++++++++
 lib/ft8/sync8.f90                 |  3 ++-
 4 files changed, 47 insertions(+), 1 deletion(-)
 create mode 100644 lib/ft8/get_spectrum_baseline.f90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index f578ac270..4f2a987d5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -487,6 +487,7 @@ set (wsjt_FSRCS
   lib/geodist.f90
   lib/getlags.f90
   lib/getmet4.f90
+  lib/ft8/get_spectrum_baseline.f90
   lib/ft2/gfsk_pulse.f90
   lib/graycode.f90
   lib/graycode65.f90
diff --git a/lib/ft8/baseline.f90 b/lib/ft8/baseline.f90
index 035f9ab4b..496795108 100644
--- a/lib/ft8/baseline.f90
+++ b/lib/ft8/baseline.f90
@@ -37,6 +37,7 @@ subroutine baseline(s,nfa,nfb,sbase)
   kz=k
   a=0.
   call polyfit(x,y,y,kz,nterms,0,a,chisqr)  !Fit a low-order polynomial
+  sbase=0.0
   do i=ia,ib
      t=i-i0
      sbase(i)=a(1)+t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5))))) + 0.65
diff --git a/lib/ft8/get_spectrum_baseline.f90 b/lib/ft8/get_spectrum_baseline.f90
new file mode 100644
index 000000000..843e08c2c
--- /dev/null
+++ b/lib/ft8/get_spectrum_baseline.f90
@@ -0,0 +1,43 @@
+subroutine get_spectrum_baseline(dd,nfa,nfb,sbase)
+
+  include 'ft8_params.f90'
+  parameter(NST=NFFT1/2,NF=NMAX/NST-1)
+  real s(NH1,NF)
+  real savg(NH1)
+  real sbase(NH1)
+  real x(NFFT1)
+  real window(NFFT1)
+  complex cx(0:NH1)
+  real dd(NMAX)
+  equivalence (x,cx)
+  logical first
+  data first/.true./
+  save first,window
+
+  if(first) then
+    first=.false.
+    pi=4.0*atan(1.)
+    window=0.
+    call nuttal_window(window,NFFT1)
+    window=window/sum(window)*NSPS*2/300.0
+  endif
+
+! Compute symbol spectra, stepping by NSTEP steps.  
+  savg=0.
+  df=12000.0/NFFT1  
+  do j=1,NF
+     ia=(j-1)*NST + 1
+     ib=ia+NFFT1-1
+     if(ib.gt.NMAX) exit
+     x=dd(ia:ib)*window
+     call four2a(x,NFFT1,1,-1,0)              !r2c FFT
+     s(1:NH1,j)=abs(cx(1:NH1))**2
+     savg=savg + s(1:NH1,j)                   !Average spectrum
+  enddo
+
+  if(nfa.lt.nint(200.0)) nfa=nint(200.0)
+  if(nfb.gt.nint(4910.0)) nfb=nint(4910.0)
+  call baseline(savg,nfa,nfb,sbase)
+
+return
+end subroutine get_spectrum_baseline
diff --git a/lib/ft8/sync8.f90 b/lib/ft8/sync8.f90
index 219c6ebf1..3d2141c02 100644
--- a/lib/ft8/sync8.f90
+++ b/lib/ft8/sync8.f90
@@ -37,7 +37,8 @@ subroutine sync8(dd,nfa,nfb,syncmin,nfqso,maxcand,s,candidate,   &
      enddo
      savg=savg + s(1:NH1,j)                   !Average spectrum
   enddo
-  call baseline(savg,nfa,nfb,sbase)
+!  call baseline(savg,nfa,nfb,sbase)
+  call get_spectrum_baseline(dd,nfa,nfb,sbase)
 
   ia=max(1,nint(nfa/df))
   ib=nint(nfb/df)