FT4: Implement polynomial baseline fit.

CMakeLists.txt

@@ -387,6 +387,7 @@ set (wsjt_FSRCS
   lib/azdist.f90
   lib/badmsg.f90
   lib/ft8/baseline.f90
+  lib/ft4/ft4_baseline.f90
   lib/bpdecode40.f90
   lib/bpdecode128_90.f90
   lib/ft8/bpdecode174_91.f90

lib/ft4/ft4_baseline.f90

@@ -0,0 +1,49 @@
+subroutine ft4_baseline(s,nfa,nfb,sbase)
+
+! Fit baseline to spectrum (for FT8)
+! Input:  s(npts)         Linear scale in power
+! Output: sbase(npts)    Baseline
+
+  include 'ft4_params.f90'
+  implicit real*8 (a-h,o-z)
+  real*4 s(NH1)
+  real*4 sbase(NH1)
+  real*4 base
+  real*8 x(1000),y(1000),a(5)
+  data nseg/10/,npct/10/
+
+  df=12000.0/NFFT1                    !5.21 Hz
+  ia=max(nint(200.0/df),nfa)
+  ib=min(NH1,nfb)
+  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)
+     sbase(i)=10**(sbase(i)/10.0)
+  enddo
+  return
+end subroutine ft4_baseline

lib/ft4/ft4sim_mult.f90

@@ -17,8 +17,6 @@ program ft4sim_mult
   integer itone(NN)
   integer*1 msgbits(77)
   integer*2 iwave(NZZ)                  !Generated full-length waveform
-  integer icos4(4)
-  data icos4/0,1,3,2/
   
 ! Get command-line argument(s)
   nargs=iargc()

lib/ft4/getcandidates4.f90

@@ -10,7 +10,6 @@ subroutine getcandidates4(dd,fa,fb,syncmin,nfqso,maxcand,savg,candidate,   &
   complex cx(0:NH1)
   real candidate(3,maxcand)
   real dd(NMAX)
-  integer indx(NH1)
   integer ipk(1)
   equivalence (x,cx)
   logical first
@@ -26,7 +25,6 @@ subroutine getcandidates4(dd,fa,fb,syncmin,nfqso,maxcand,savg,candidate,   &
 
 ! Compute symbol spectra, stepping by NSTEP steps.  
   savg=0.
-  tstep=NSTEP/12000.0                         
   df=12000.0/NFFT1                            
   fac=1.0/300.0
   do j=1,NHSYM
@@ -40,27 +38,20 @@ subroutine getcandidates4(dd,fa,fb,syncmin,nfqso,maxcand,savg,candidate,   &
      enddo
      savg=savg + s(1:NH1,j)                   !Average spectrum
   enddo
+  savg=savg/NHSYM
   savsm=0.
   do i=8,NH1-7
     savsm(i)=sum(savg(i-7:i+7))/15.
   enddo
+
   nfa=fa/df
-  if(nfa.lt.1) nfa=1
+  if(nfa.lt.8) nfa=8
   nfb=fb/df
   if(nfb.gt.nint(5000.0/df)) nfb=nint(5000.0/df)
-  n300=300/df
-  n2500=2500/df
-!  np=nfb-nfa+1
-  np=n2500-n300+1
-  indx=0
-  call indexx(savsm(n300:n2500),np,indx)
-  xn=savsm(n300+indx(nint(0.3*np)))
   ncand=0
-  if(xn.le.1.e-8) return 
-  savsm=savsm/xn
-!  call ft4_baseline(savg,nfa,nfb,sbase)
-!  savsm=savsm/sbase
-
+  call ft4_baseline(savg,nfa,nfb,sbase)
+  if(any(sbase(nfa:nfb).le.0)) return 
+  savsm(nfa:nfb)=savsm(nfa:nfb)/sbase(nfa:nfb)
   f_offset = -1.5*12000.0/NSPS
   do i=nfa+1,nfb-1
      if(savsm(i).ge.savsm(i-1) .and. savsm(i).ge.savsm(i+1) .and.      &