Integrate the new Echo SNR algorithm into WSJT-X.

lib/avecho.f90

@@ -1,4 +1,4 @@
-subroutine avecho(id2,ndop,nfrit,nqual,f1,xlevel,sigdb,snr,dfreq,width)
+subroutine avecho(id2,ndop,nfrit,nqual,f1,xlevel,snrdb,db_err,dfreq,width)
 
   integer TXLENGTH
   parameter (TXLENGTH=27648)           !27*1024
@@ -13,11 +13,22 @@ subroutine avecho(id2,ndop,nfrit,nqual,f1,xlevel,sigdb,snr,dfreq,width)
   real s(8192)
   real x(NFFT)
   integer ipkv(1)
+  logical ex
   complex c(0:NH)
   equivalence (x,c),(ipk,ipkv)
   common/echocom/nclearave,nsum,blue(NZ),red(NZ)
+  common/echocom2/echo_spread
   save dop0,sa,sb
 
+  fspread=echo_spread                !### Use the predicted Doppler spread ###
+  inquire(file='fspread.txt',exist=ex)
+  if(ex) then
+     open(39,file='fspread.txt',status='old')
+     read(39,*) fspread
+     close(39)
+  endif
+  width=fspread
+
   dop=ndop
   sq=0.
   do i=1,TXLENGTH
@@ -48,69 +59,28 @@ subroutine avecho(id2,ndop,nfrit,nqual,f1,xlevel,sigdb,snr,dfreq,width)
   if(ia.gt.7590 .or. ib.gt.7590) go to 900
 
   nsum=nsum+1
-
   do i=1,NZ
      sa(i)=sa(i) + s(ia+i-2048)    !Center at initial doppler freq
      sb(i)=sb(i) + s(ib+i-2048)    !Center at expected echo freq
   enddo
-
+  call echo_snr(sa,sb,fspread,blue,red,snrdb,db_err,dfreq,snr_detect)
-  call pctile(sb,200,50,r0)
+  nqual=snr_detect-2
-  call pctile(sb(1800),200,50,r1)
+  if(nqual.lt.0) nqual=0
-
-  sum=0.
-  sq=0.
-  do i=1,NZ
-     y=r0 + (r1-r0)*(i-100.0)/1800.0
-     blue(i)=sa(i)/y
-     red(i)=sb(i)/y
-     if(i.le.500 .or. i.ge.3597) then
-        sum=sum+red(i)
-        sq=sq + (red(i)-1.0)**2
-     endif
-  enddo
-  ave=sum/1000.0
-  rms=sqrt(sq/1000.0)
-
-  redmax=maxval(red)
-  ipkv=maxloc(red)
-  fac=10.0/max(redmax,10.0)
-  dfreq=(ipk-2048)*df
-  snr=(redmax-ave)/rms
-
-  sigdb=-99.0
-  if(ave.gt.0.0) sigdb=10.0*log10(redmax/ave - 1.0) - 35.7
-
-  nqual=0
-  if(nsum.ge.2 .and. nsum.lt.4)  nqual=(snr-4)/5
-  if(nsum.ge.4 .and. nsum.lt.8)  nqual=(snr-3)/4
-  if(nsum.ge.8 .and. nsum.lt.12) nqual=(snr-3)/3
-  if(nsum.ge.12) nqual=(snr-2.5)/2.5
-  if(nqual.lt.0)  nqual=0
   if(nqual.gt.10) nqual=10
 
 ! Scale for plotting
+  redmax=maxval(red)
+  fac=10.0/max(redmax,10.0)
   blue=fac*blue
   red=fac*red
-
+  nsmo=max(0.0,0.25*width/df)
-  sum=0.
-  do i=ipk,ipk+300
-     if(i.gt.NZ) exit
-     if(red(i).lt.1.0) exit
-     sum=sum+(red(i)-1.0)
-  enddo
-  do i=ipk-1,ipk-300,-1
-     if(i.lt.1) exit
-     if(red(i).lt.1.0) exit
-     sum=sum+(red(i)-1.0)
-  enddo
-  bins=sum/(red(ipk)-1.0)
-  width=df*bins
-  nsmo=max(0.0,0.25*bins)
-
   do i=1,nsmo
      call smo121(red,NZ)
      call smo121(blue,NZ)
   enddo
 
+!  write(*,3001) snrdb,db_err,dfreq,snr_detect,redmax,nqual,nsmo,nclearave,nsum
+!3001 format('A',5f10.1,4i4)
+
 900  return
 end subroutine avecho

lib/echo_snr.f90

@@ -1,4 +1,4 @@
-subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
+subroutine echo_snr(sa,sb,fspread,blue,red,snrdb,db_err,fpeak,snr_detect)
 
   parameter (NZ=4096)
   real sa(NZ)
@@ -9,33 +9,27 @@ subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
   equivalence (ipk,ipkv)
 
   df=12000.0/32768.0
-  wh=0.5*fspread0+10.0
+  wh=0.5*fspread+10.0
   i1=nint((1500.0 - 2.0*wh)/df) - 2048
   i2=nint((1500.0 - wh)/df) - 2048
   i3=nint((1500.0 + wh)/df) - 2048
   i4=nint((1500.0 + 2.0*wh)/df) - 2048
 
-  call pctile(sb(i1),i2-i1,50,r0)
+!  call pctile(sb(i1),i2-i1,50,r0)
-  call pctile(sb(i3+1),i4-i3,50,r1)
+!  call pctile(sb(i3+1),i4-i3,50,r1)
-  scale=(r1-r0)/(0.5*(i3+i4) - 0.5*(i1+i2))
+!  ave=0.5*(r0+r1)
-  do i=1,NZ
+!  blue=sa/ave
-     x=i
+!  red=sb/ave
-     y=r0 + scale*(i-0.5*(i1+i2))
-     blue(i)=sa(i)/y
-     red(i)=sb(i)/y
-!     write(13,1100) i,i*df+750.0,sb(i),y,red(i)
-!1100 format(i6,f10.2,2f10.2,f10.6)
-enddo
 
-!  baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3)
+  baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3)
-!  blue=sa/baseline
+  blue=sa/baseline
-!  red=sb/baseline
+  red=sb/baseline
   psig=sum(red(i2:i3)-1.0)
   pnoise_2500 = 2500.0/df
   snrdb=db(psig/pnoise_2500)
 
   smax=0.
-  mh=max(1,nint(0.2*fspread0/df))
+  mh=max(1,nint(0.2*fspread/df))
   do i=i2,i3
      ssum=sum(red(i-mh:i+mh))
      if(ssum.gt.smax) then
@@ -51,6 +45,7 @@ enddo
   snr_detect=psig/perr
   db_err=99.0
   if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500)
+  if(db_err.lt.0.5) db_err=0.5
 
   return
 end subroutine echo_snr

widgets/mainwindow.cpp

@@ -1529,7 +1529,7 @@ void MainWindow::dataSink(qint64 frames)
   }
   if(bCallDecoder) {
     if(m_mode=="Echo") {
-      float snr=0;
+      float dBerr=0.0;
       int nfrit=0;
       int nqual=0;
       float f1=1500.0;
@@ -1540,10 +1540,10 @@ void MainWindow::dataSink(qint64 frames)
       echocom_.nclearave=m_nclearave;
       int nDop=0;
       avecho_(dec_data.d2,&nDop,&nfrit,&nqual,&f1,&xlevel,&sigdb,
-          &snr,&dfreq,&width);
+          &dBerr,&dfreq,&width);
       QString t;
-      t = t.asprintf("%3d %7.1f %7.1f %7.1f %7.1f %3d",echocom_.nsum,xlevel,sigdb,
+      t = t.asprintf("%3d %7.1f %7.1f %7.1f %7.1f %7.1f %3d",echocom_.nsum,xlevel,sigdb,
-                dfreq,width,nqual);
+                dBerr,dfreq,width,nqual);
       t=QDateTime::currentDateTimeUtc().toString("hh:mm:ss  ") + t;
       if (ui) ui->decodedTextBrowser->appendText(t);
       if(m_echoGraph->isVisible()) m_echoGraph->plotSpec();
@@ -6645,7 +6645,7 @@ void MainWindow::on_actionEcho_triggered()
   m_bFastMode=false;
   m_bFast9=false;
   WSPR_config(true);
-  ui->lh_decodes_headings_label->setText("   UTC      N   Level    Sig      DF    Width   Q");
+  ui->lh_decodes_headings_label->setText("   UTC      N   Level    SNR     dBerr    DF    Width   Q");
   //                       01234567890123456789012345678901234567
   displayWidgets(nWidgets("00000000000000000000001000000000000000"));
   fast_config(false);