Use a baseline-fitting procedure to improve S/N estimates for FT8.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@8036 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

CMakeLists.txt

@@ -366,6 +366,7 @@ set (wsjt_FSRCS
   lib/averms.f90
   lib/azdist.f90
   lib/badmsg.f90
+  lib/fsk4hf/baseline.f90
   lib/bpdecode40.f90
   lib/bpdecode144.f90
   lib/fsk4hf/bpdecode120.f90

lib/fsk4hf/baseline.f90

@@ -0,0 +1,48 @@
+subroutine baseline(s,nfa,nfb,sbase)
+
+! Fit baseline to spectrum (for FT8)
+! Input:  s(npts)         Linear scale in power
+! Output: sbase(npts)    Baseline
+
+  implicit real*8 (a-h,o-z)
+  real*4 s(1920)
+  real*4 sbase(1920)
+  real*4 base
+  real*8 x(1000),y(1000),a(5)
+  data nseg/10/,npct/10/
+
+  df=12000.0/3840.0                    !3.125 Hz
+  ia=nint(nfa/df)
+  ib=nint(nfb/df)
+  do i=ia,ib
+     s(i)=10.0*log10(s(i))            !Convert to dB scale
+  enddo
+
+  nterms=5
+  nlen=(ib-ia+1)/nseg                 !Length of test segment
+  i0=(ib-ia+1)/2                      !Midpoint
+  k=0
+  do n=1,nseg                         !Loop over all segments
+     ja=ia + (n-1)*nlen
+     jb=ja+nlen-1
+     call pctile(s(ja),nlen,npct,base) !Find lowest npct of points
+     do i=ja,jb
+        if(s(i).le.base) then
+           if (k.lt.1000) k=k+1       !Save all "lower envelope" points
+           x(k)=i-i0
+           y(k)=s(i)
+        endif
+     enddo
+  enddo
+  kz=k
+  a=0.
+  call polyfit(x,y,y,kz,nterms,0,a,chisqr)  !Fit a low-order polynomial
+  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
+!     write(51,3051) i*df,s(i),sbase(i)
+!3051 format(3f12.3)
+  enddo
+
+  return
+end subroutine baseline

lib/fsk4hf/ft8b.f90

@@ -1,6 +1,6 @@
 subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid,       &
-     lsubtract,nagain,iaptype,mygrid6,bcontest,sync0,f1,xdt,apsym,nharderrors,&
-     dmin,nbadcrc,ipass,iera,message,xsnr)  
+     lsubtract,nagain,iaptype,mygrid6,bcontest,sync0,f1,xdt,xbase,apsym,      &
+     nharderrors,dmin,nbadcrc,ipass,iera,message,xsnr)  
 
   use timer_module, only: timer
   include 'ft8_params.f90'
@@ -337,6 +337,12 @@ subroutine ft8b(dd0,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid,       &
         xsnr=0.001
         if(xnoi.gt.0 .and. xnoi.lt.xsig) xsnr=xsig/xnoi-1.0
         xsnr=10.0*log10(xsnr)-27.0
+!###
+        xsnr2=db(xsig/xbase - 1.0) - 32.0
+!        write(52,3052) f1,xdt,xsig,xnoi,xbase,xsnr,xsnr2
+!3052    format(7f10.2)
+        xsnr=xsnr2
+!###
         if(xsnr .lt. -24.0) xsnr=-24.0
         return
      endif

lib/fsk4hf/sync8.f90

@@ -1,4 +1,4 @@
-subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand)
+subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand,sbase)
 
   include 'ft8_params.f90'
 ! Search over +/- 1.5s relative to 0.5s TX start time. 
@@ -6,6 +6,7 @@ subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand)
   complex cx(0:NH1)
   real s(NH1,NHSYM)
   real savg(NH1)
+  real sbase(NH1)
   real x(NFFT1)
   real sync2d(NH1,-JZ:JZ)
   real red(NH1)
@@ -35,7 +36,8 @@ subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand)
      enddo
      savg=savg + s(1:NH1,j)                   !Average spectrum
   enddo
-  savg=savg/NHSYM
+  call baseline(savg,nfa,nfb,sbase)
+!  savg=savg/NHSYM
 !  do i=1,NH1
 !     write(51,3051) i*df,savg(i),db(savg(i))
 !3051 format(f10.3,e12.3,f12.3)

lib/ft8_decode.f90

@@ -34,6 +34,7 @@ contains
     class(ft8_decoder), intent(inout) :: this
     procedure(ft8_decode_callback) :: callback
     real s(NH1,NHSYM)
+    real sbase(NH1)
     real candidate(3,200)
     real dd(15*12000)
     logical, intent(in) :: lapon,nagain
@@ -86,17 +87,18 @@ contains
       endif 
 
       call timer('sync8   ',0)
-      call sync8(dd,ifa,ifb,syncmin,nfqso,s,candidate,ncand)
+      call sync8(dd,ifa,ifb,syncmin,nfqso,s,candidate,ncand,sbase)
       call timer('sync8   ',1)
       do icand=1,ncand
         sync=candidate(3,icand)
         f1=candidate(1,icand)
         xdt=candidate(2,icand)
+        xbase=10.0**(0.1*(sbase(nint(f1/3.125))-40.0))
         nsnr0=min(99,nint(10.0*log10(sync) - 25.5))    !### empirical ###
         call timer('ft8b    ',0)
         call ft8b(dd,newdat,nQSOProgress,nfqso,nftx,ndepth,lapon,napwid,    &
-             lsubtract,nagain,iaptype,mygrid6,bcontest,sync,f1,xdt,apsym,   &
-             nharderrors,dmin,nbadcrc,iappass,iera,message,xsnr)
+             lsubtract,nagain,iaptype,mygrid6,bcontest,sync,f1,xdt,xbase,   &
+             apsym,nharderrors,dmin,nbadcrc,iappass,iera,message,xsnr)
         nsnr=nint(xsnr) 
         xdt=xdt-0.5
         hd=nharderrors+dmin