Use 3rd order polynomial fit to estimate the noise baseline. The

polynomial fit is done over 400 Hz bandwidth for T/R periods longer
than 15s, and over approx. 600 Hz (10 times the signal bandwidth) for
T/R period of 15s.

CMakeLists.txt

@@ -618,6 +618,7 @@ set (wsjt_FSRCS
   lib/fst4/osd240_101.f90
   lib/fst4/osd240_74.f90
   lib/fst4/get_crc24.f90
+  lib/fst4/fst4_baseline.f90
   )
 
 # temporary workaround for a gfortran v7.3 ICE on Fedora 27 64-bit

lib/fst4/fst4_baseline.f90

@@ -0,0 +1,48 @@
+subroutine fst4_baseline(s,np,ia,ib,npct,sbase)
+
+! Fit baseline to spectrum (for FST4)
+! Input:  s(npts)         Linear scale in power
+! Output: sbase(npts)    Baseline
+
+  implicit real*8 (a-h,o-z)
+  real*4 s(np),sw(np)
+  real*4 sbase(np)
+  real*4 base
+  real*8 x(1000),y(1000),a(5)
+  data nseg/8/
+
+  do i=ia,ib
+     sw(i)=10.0*log10(s(i))            !Convert to dB scale
+  enddo
+
+  nterms=3
+  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(sw(ja),nlen,npct,base) !Find lowest npct of points
+     do i=ja,jb
+        if(sw(i).le.base) then
+           if (k.lt.1000) k=k+1       !Save all "lower envelope" points
+           x(k)=i-i0
+           y(k)=sw(i)
+        endif
+     enddo
+  enddo
+  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))) + 0.2
+!     write(51,3051) i,sw(i),sbase(i)
+!3051 format(i8,2f12.3)
+  enddo
+
+  sbase=10**(sbase/10.0)
+
+  return
+end subroutine fst4_baseline

lib/fst4_decode.f90

@@ -49,7 +49,7 @@ contains
       complex, allocatable :: cframe(:)
       complex, allocatable :: c_bigfft(:)          !Complex waveform
       real llr(240),llrs(240,4)
-      real candidates(200,4)
+      real candidates(200,5)
       real bitmetrics(320,4)
       real s4(0:3,NN)
       real minsync
@@ -254,14 +254,19 @@ contains
 
       nhicoh=1
       nsyncoh=8
-      fa=max(100,nint(nfqso+1.5*hmod*baud-ntol))
+      if(iwspr.eq.1) then
-      fb=min(4800,nint(nfqso+1.5*hmod*baud+ntol))
+         fa=1400.0
-      minsync=1.2
+         fb=1600.0
+      else
+         fa=max(100,nint(nfqso+1.5*hmod*baud-ntol))
+         fb=min(4800,nint(nfqso+1.5*hmod*baud+ntol))
+      endif
+      minsync=1.20
       if(ntrperiod.eq.15) minsync=1.15
 
 ! Get first approximation of candidate frequencies
       call get_candidates_fst4(c_bigfft,nfft1,nsps,hmod,fs,fa,fb,     &
-         minsync,ncand,candidates,base)
+         minsync,ncand,candidates)
 
       ndecodes=0
       decodes=' '
@@ -317,7 +322,7 @@ contains
       enddo
       ncand=ic
       xsnr=0.
-
+!write(*,*) 'ncand ',ncand
       do icand=1,ncand
          sync=candidates(icand,2)
          fc_synced=candidates(icand,3)
@@ -465,6 +470,7 @@ contains
                      do i=1,NN
                         xsig=xsig+s4(itone(i),i)
                      enddo
+                     base=candidates(icand,5)
                      arg=600.0*(xsig/base)-1.0
                      if(arg.gt.0.0) then
                         xsnr=10*log10(arg)-35.5-12.5*log10(nsps/8200.0)
@@ -645,14 +651,15 @@ contains
    end subroutine fst4_downsample
 
    subroutine get_candidates_fst4(c_bigfft,nfft1,nsps,hmod,fs,fa,fb,   &
-      minsync,ncand,candidates,base)
+      minsync,ncand,candidates)
 
       complex c_bigfft(0:nfft1/2)              !Full length FFT of raw data
       integer hmod                             !Modulation index (submode)
       integer im(1)                            !For maxloc
-      real candidates(200,4)                   !Candidate list
+      real candidates(200,5)                   !Candidate list
       real, allocatable :: s(:)                !Low resolution power spectrum
       real, allocatable :: s2(:)               !CCF of s() with 4 tones
+      real, allocatable :: sbase(:)            !noise baseline estimate
       real xdb(-3:3)                           !Model 4-tone CCF peaks
       real minsync
       data xdb/0.25,0.50,0.75,1.0,0.75,0.50,0.25/
@@ -668,17 +675,17 @@ contains
       signal_bw=4*(12000.0/nsps)*hmod
       analysis_bw=min(4800.0,fb)-max(100.0,fa)
       xnoise_bw=10.0*signal_bw                  !Is this a good compromise?
-      if(analysis_bw.gt.xnoise_bw) then
+      if(xnoise_bw .lt. 400.0) xnoise_bw=400.0  
-         ina=ia
+      if(analysis_bw.gt.xnoise_bw) then         !Estimate noise baseline over analysis bw
-         inb=ib
+         ina=0.9*ia
-      else
+         inb=min(int(1.1*ib),nfft1/2)
-         fcenter=(fa+fb)/2.0                      !If noise_bw > analysis_bw,
+      else                                      !Estimate noise baseline over noise bw
-         fl = max(100.0,fcenter-xnoise_bw/2.)/df2  !we'll search over noise_bw
+         fcenter=(fa+fb)/2.0                    
+         fl = max(100.0,fcenter-xnoise_bw/2.)/df2  
          fh = min(4800.0,fcenter+xnoise_bw/2.)/df2
          ina=nint(fl)
          inb=nint(fh)
       endif
-
       nnw=nint(48000.*nsps*2./fs)
       allocate (s(nnw))
       s=0.                                  !Compute low-resolution power spectrum
@@ -692,12 +699,16 @@ contains
       ina=max(ina,1+3*hmod)                       !Don't run off the ends
       inb=min(inb,nnw-3*hmod)
       allocate (s2(nnw))
+      allocate (sbase(nnw))
       s2=0.
       do i=ina,inb                                !Compute CCF of s() and 4 tones
          s2(i)=s(i-hmod*3) + s(i-hmod) +s(i+hmod) +s(i+hmod*3)
       enddo
-      call pctile(s2(ina+hmod*3:inb-hmod*3),inb-ina+1-hmod*6,30,base)
+      npct=30
-      s2=s2/base                                  !Normalize wrt noise level
+      call fst4_baseline(s2,nnw,ina+hmod*3,inb-hmod*3,npct,sbase)
+      if(any(sbase(ina:inb).le.0.0)) return
+      s2(ina:inb)=s2(ina:inb)/sbase(ina:inb)             !Normalize wrt noise level
+
       ncand=0
       candidates=0
       if(ia.lt.3) ia=3
@@ -717,12 +728,11 @@ contains
                s2(k)=max(0.,s2(k)-0.9*pval*xdb(i))
             endif
          enddo
-!         s2(max(1,iploc-2*hmod*3):min(nnw,iploc+2*hmod*3))=0.0
          ncand=ncand+1
          candidates(ncand,1)=df2*iploc         !Candidate frequency
          candidates(ncand,2)=pval              !Rough estimate of SNR
+         candidates(ncand,5)=sbase(iploc)
       enddo
-
       return
    end subroutine get_candidates_fst4