Final (?) version of test_snr, with echo_snr() moved into wsjt Fortran library.

CMakeLists.txt

@@ -382,6 +382,7 @@ set (wsjt_FSRCS
   lib/demod64a.f90
   lib/determ.f90
   lib/downsam9.f90
+  lib/echo_snr.f90
   lib/encode232.f90
   lib/encode4.f90
   lib/encode_msk40.f90

lib/echo_snr.f90

@@ -0,0 +1,56 @@
+subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
+
+  parameter (NZ=4096)
+  real sa(NZ)
+  real sb(NZ)
+  real blue(NZ)
+  real red(NZ)
+  integer ipkv(1)
+  equivalence (ipk,ipkv)
+
+  df=12000.0/32768.0
+  wh=0.5*fspread0+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(i3+1),i4-i3,50,r1)
+  scale=(r1-r0)/(0.5*(i3+i4) - 0.5*(i1+i2))
+  do i=1,NZ
+     x=i
+     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)
+!  blue=sa/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))
+  do i=i2,i3
+     ssum=sum(red(i-mh:i+mh))
+     if(ssum.gt.smax) then
+        smax=ssum
+        ipk=i
+     endif
+  enddo
+  fpeak=ipk*df - 750.0
+
+  call averms(red(i1:i2-1),i2-i1,-1,ave1,rms1)
+  call averms(red(i3+1:i4),i4-i3,-1,ave2,rms2)
+  perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3))
+  snr_detect=psig/perr
+  db_err=99.0
+  if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500)
+
+  return
+end subroutine echo_snr

lib/test_snr.f90

@@ -35,13 +35,12 @@ program test_snr
   nfsample=h%nsamrate
   df=12000.0/NFFT
   s=0.
+  fac=1.0/NMAX
   
   do iping=1,npings
      read(10) id2(1:NMAX)
-     x(1:NMAX)=id2(1:NMAX)
-     x(NMAX+1:)=0.
-     rms=sqrt(dot_product(x(1:NMAX),x(1:NMAX))/NMAX)
-     
+     x(1:NMAX)=fac*id2(1:NMAX)
+     x(NMAX+1:)=0.     
      call four2a(x,NFFT,1,-1,0)
      do i=1,8192                             !Accumulate spectrum 0 - 3 kHz
         s(i)=s(i) + real(c(i))**2 + aimag(c(i))**2
@@ -53,14 +52,7 @@ program test_snr
 
   call echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
   
-  nsum=npings
-  nqual=0
-  if(nsum.ge.2 .and. nsum.lt.4)  nqual=(snr_detect-4)/5
-  if(nsum.ge.4 .and. nsum.lt.8)  nqual=(snr_detect-3)/4
-  if(nsum.ge.8 .and. nsum.lt.12) nqual=(snr_detect-3)/3
-  if(nsum.ge.12) nqual=(snr_detect-2.5)/2.5
-  if(nqual.lt.0)  nqual=0
-  if(nqual.gt.10) nqual=10
+  nqual=min(10,int(snr_detect-4.0))
 
   write(*,1010) fspread0,snrdb0,snrdb,snrdb-snrdb0,db_err,fpeak,   &
        snr_detect,nqual
@@ -69,52 +61,7 @@ program test_snr
 
   do i=1,8192
      write(12,1100) i*df,s(i)
-1100 format(f10.3,e12.3)
+1100 format(f10.3,e15.6)
   enddo
 
 999 end program test_snr
-
-subroutine echo_snr(sa,sb,fspread0,blue,red,snrdb,db_err,fpeak,snr_detect)
-
-  parameter (NZ=4096)
-  real sa(NZ)
-  real sb(NZ)
-  real blue(NZ)
-  real red(NZ)
-  integer ipkv(1)
-  equivalence (ipk,ipkv)
-
-  df=12000.0/32768.0
-  wh=0.5*fspread0+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
-
-  baseline=(sum(sb(i1:i2-1)) + sum(sb(i3+1:i4)))/(i2+i4-i1-i3)
-  blue=sa/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))
-  do i=i2,i3
-     ssum=sum(red(i-mh:i+mh))
-     if(ssum.gt.smax) then
-        smax=ssum
-        ipk=i
-     endif
-  enddo
-  fpeak=ipk*df - 750.0
-
-  call averms(red(i1:i2-1),i2-i1,-1,ave1,rms1)
-  call averms(red(i3+1:i4),i4-i3,-1,ave2,rms2)
-  perr=0.707*(rms1+rms2)*sqrt(float(i2-i1+i4-i3))
-  snr_detect=psig/perr
-  db_err=99.0
-  if(psig.gt.perr) db_err=snrdb - db((psig-perr)/pnoise_2500)
-
-  return
-end subroutine echo_snr