From 021336823af65fe6f17fe5c689930d82e32a8e90 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Sat, 26 Aug 2017 16:30:47 +0000 Subject: [PATCH] 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 | 1 + lib/fsk4hf/baseline.f90 | 48 +++++++++++++++++++++++++++++++++++++++++ lib/fsk4hf/ft8b.f90 | 10 +++++++-- lib/fsk4hf/sync8.f90 | 6 ++++-- lib/ft8_decode.f90 | 8 ++++--- 5 files changed, 66 insertions(+), 7 deletions(-) create mode 100644 lib/fsk4hf/baseline.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 5aa7a9932..694371d97 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/lib/fsk4hf/baseline.f90 b/lib/fsk4hf/baseline.f90 new file mode 100644 index 000000000..fac3f2ccb --- /dev/null +++ b/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 diff --git a/lib/fsk4hf/ft8b.f90 b/lib/fsk4hf/ft8b.f90 index 308f04206..8dc307f21 100644 --- a/lib/fsk4hf/ft8b.f90 +++ b/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 diff --git a/lib/fsk4hf/sync8.f90 b/lib/fsk4hf/sync8.f90 index 0813d7f88..a13297318 100644 --- a/lib/fsk4hf/sync8.f90 +++ b/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) diff --git a/lib/ft8_decode.f90 b/lib/ft8_decode.f90 index 20cbaf863..1ed94473b 100644 --- a/lib/ft8_decode.f90 +++ b/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