Add submode to argument list of JT9 subroutine decode().

Add test program fersum; several tweaks to fer65.
Code cleanup in lorentzian, and return chisqr in a(5).


git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6567 ab8295b8-cf94-4d9e-aec4-7959e3be5d79

lib/decoder.f90

@@ -93,9 +93,10 @@ subroutine multimode_decoder(ss,id2,params,nfsample)
   else if(params%nmode.eq.9 .or. (params%nmode.eq.(65+9) .and. params%ntxmode.eq.9)) then
      ! We're in JT9 mode, or should do JT9 first
      call timer('decjt9  ',0)
-     call my_jt9%decode(jt9_decoded,ss,id2,params%nutc,params%nfqso,newdat9,params%npts8,   &
+     call my_jt9%decode(jt9_decoded,ss,id2,params%nutc,params%nfqso,       &
-          params%nfa,params%nfsplit,params%nfb,params%ntol,params%nzhsym,                   &
+          newdat9,params%npts8,params%nfa,params%nfsplit,params%nfb,       &
-          logical(params%nagain),params%ndepth,params%nmode)
+          params%ntol,params%nzhsym,logical(params%nagain),params%ndepth,  &
+          params%nmode,params%nsubmode)
      call timer('decjt9  ',1)
   endif
 
@@ -114,9 +115,10 @@ subroutine multimode_decoder(ss,id2,params,nfsample)
         call timer('jt65a   ',1)
      else
         call timer('decjt9  ',0)
-        call my_jt9%decode(jt9_decoded,ss,id2,params%nutc,params%nfqso,newdat9,params%npts8,&
+        call my_jt9%decode(jt9_decoded,ss,id2,params%nutc,params%nfqso,    &
-             params%nfa,params%nfsplit,params%nfb,params%ntol,params%nzhsym,                &
+             newdat9,params%npts8,params%nfa,params%nfsplit,params%nfb,    &
-             logical(params%nagain),params%ndepth,params%nmode)
+             params%ntol,params%nzhsym,logical(params%nagain),             &
+             params%ndepth,params%nmode,params%nsubmode)
         call timer('decjt9  ',1)
      end if
   endif

lib/fchisq0.f90

@@ -2,7 +2,7 @@ real function fchisq0(y,npts,a)
 
   real y(npts),a(4)
   
-  rewind 51
+!  rewind 51
   chisq = 0.
   do i=1,npts
      x=i

lib/fer65.f90

@@ -20,7 +20,7 @@ program fer65
   implicit real*8 (a-h,o-z)
   real*8 s(7),sq(7)
   character arg*12,cmnd*100,decoded*22,submode*1,csync*1,f1*15,f2*15
-  character*10 outfile
+  character*12 outfile
   logical syncok
 
   nargs=iargc()
@@ -44,19 +44,18 @@ program fer65
   call getarg(7,arg)
   read(arg,*) iters
 
-
   write(outfile,1001) submode,d,navg,nds
-1001 format(a1,f4.1,'_',i2.2,'_',i1)
+1001 format(a1,f6.2,'_',i2.2,'_',i1)
   if(outfile(2:2).eq.' ') outfile(2:2)='0'
+  if(outfile(3:3).eq.' ') outfile(3:3)='0'
 
   ndepth=3
   if(navg.gt.1) ndepth=ndepth+16
   if(nds.ne.0) ndepth=ndepth+32
 
-  dfmax=min(d,0.5*2.69)
+  dfmax=3
-  if(submode.eq.'b' .or. submode.eq.'B') dfmax=min(d,2.69)
+  if(submode.eq.'b' .or. submode.eq.'B') dfmax=6
-  if(submode.eq.'c' .or. submode.eq.'C') dfmax=min(d,2.0*2.69)
+  if(submode.eq.'c' .or. submode.eq.'C') dfmax=11
-  if(dfmax.lt.0.5*2.69) dfmax=0.5*2.69
 
   ntrials=1000
   naggressive=10
@@ -64,9 +63,9 @@ program fer65
   open(20,file=outfile,status='unknown')
   open(21,file='fer65.21',status='unknown')
 
-  write(20,1000) submode,iters,ntrials,naggressive,d,iand(ndepth,3),navg,nds
+  write(20,1000) submode,iters,ntrials,naggressive,d,ndepth,navg,nds
 1000 format(/'JT65',a1,'  Iters:',i5,'  T:',i6,'  Aggr:',i3,  &
-          '  Dop:',f5.1,'  Depth:',i2,'  Navg:',i3,'  DS:',i2)
+          '  Dop:',f6.2,'  Depth:',i2,'  Navg:',i3,'  DS:',i2)
   write(20,1002) 
 1002 format(/'  dB  nsync ngood nbad     sync       dsnr        ',     &
             'DT       Freq      Nsum     Width'/85('-'))
@@ -75,7 +74,6 @@ program fer65
   do isnr=0,20
      snr=snr1+isnr
      if(snr.gt.snr2) exit
-
      nsync=0
      ngood=0
      nbad=0
@@ -83,7 +81,7 @@ program fer65
      sq=0.
      do iter=1,iters
         write(cmnd,1010) submode,d,snr,navg
-1010    format('./jt65sim -n 1 -m ',a1,' -d',f6.1,' -s \\',f5.1,' -f',i3,' >devnull')
+1010    format('./jt65sim -n 1 -m ',a1,' -d',f7.2,' -s \\',f5.1,' -f',i3,' >devnull')
         call unlink('000000_????.wav')
         call system(cmnd)
         if(navg.gt.1) then

lib/fersum.f90

@@ -0,0 +1,72 @@
+program fersum
+
+  character mode*5
+  character infile*40
+  real dop(0:9)
+  real thresh(0:9,12),threshsync(0:9,12)
+  data dop/0.25,0.5,1.0,2.0,4.0,8.0,16.0,32.0,64.0,128.0/
+
+  nargs=iargc()
+  if(nargs.ne.1) then
+     print*,'Usage: fersum <infile>'
+     go to 999
+  endif
+  call getarg(1,infile)
+  open(10,file=infile,status='old')
+  thresh=0.
+  threshsync=0.
+
+  do iblk=1,999
+1    read(10,1002,end=100) mode,iters,ntot,naggr,d,navg,nds
+1002 format(a5,8x,i5,4x,i6,7x,i3,6x,f6.2,17x,i3,5x,i2)
+     write(33,*) iblk,mode
+     if(mode.eq.'     ') go to 1
+     read(10,1002)
+     read(10,1002)
+     read(10,1002)
+
+     nsync0=0
+     ngood0=0
+     xsum0=0.
+     do n=1,99
+        read(10,1010,end=100) snr,nsync,ngood,nbad,xsync,esync,dsnr,esnr,   &
+             xdt,edt,dfreq,efreq,xsum,esum,xwidth,ewidth
+1010    format(f5.1,2i6i4,2f6.1,f6.1,f5.1,f6.2,f5.2,6f5.1)
+        if(snr.eq.0.0) exit
+        if(mode(5:5).eq.'A') nmode=1
+        if(mode(5:5).eq.'B') nmode=2
+        if(mode(5:5).eq.'C') nmode=3
+        j=nint(log(d)/log(2.0) + 2.0)
+        if(navg.eq.1 .and. nds.eq.0) k=nmode
+        if(navg.eq.1 .and. nds.eq.1) k=nmode+3
+        if(navg.gt.1 .and. nds.eq.0) k=nmode+6
+        if(navg.gt.1 .and. nds.eq.1) k=nmode+9
+        if(nsync0.le.iters/2 .and. nsync.ge.iters/2) then
+           threshsync(j,k)=snr-float(nsync-iters/2)/(nsync-nsync0)
+        endif
+        if(ngood0.le.iters/2 .and. ngood.ge.iters/2) then
+           threshold=snr-float(ngood-iters/2)/(ngood-ngood0)
+           xsumavg=max(1.0,0.5*(xsum0+xsum))
+!           write(*,1020) mode,iters,ntot,naggr,d,navg,nds,threshold,xsumavg
+!1020       format(a5,i7,i7,i3,f7.2,i3,i3,f7.1,f6.1)
+           thresh(j,k)=threshold
+        endif
+        nsync0=nsync
+        ngood0=ngood
+        xsum0=xsum
+     enddo
+  enddo
+
+100 write(12,1100)
+1100 format(' ')
+  do i=0,9
+     write(12,1110) dop(i),thresh(i,1:12)
+1110 format(f6.2,13f6.1)
+  enddo
+
+  write(12,1110)
+  do i=0,9
+     write(12,1110) dop(i),threshsync(i,1:12)
+  enddo
+
+999 end program fersum

lib/jt9_decode.f90

@@ -7,7 +7,8 @@ module jt9_decode
   end type jt9_decoder
 
   abstract interface
-     subroutine jt9_decode_callback (this, utc, sync, snr, dt, freq, drift, decoded)
+     subroutine jt9_decode_callback (this, utc, sync, snr, dt, freq, drift, &
+          decoded)
        import jt9_decoder
        implicit none
        class(jt9_decoder), intent(inout) :: this
@@ -23,8 +24,8 @@ module jt9_decode
 
 contains
 
-  subroutine decode(this,callback,ss,id2,nutc,nfqso,newdat,npts8,nfa,nfsplit,nfb,ntol,  &
+  subroutine decode(this,callback,ss,id2,nutc,nfqso,newdat,npts8,nfa,    &
-       nzhsym,nagain,ndepth,nmode)
+       nfsplit,nfb,ntol,nzhsym,nagain,ndepth,nmode,nsubmode)
     use timer_module, only: timer
 
     include 'constants.f90'
@@ -42,6 +43,8 @@ contains
     common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
     save ccfred,red2
 
+    write(60) nutc,nfqso,ntol,ndepth,nmode,nsubmode,nzhsym,ss,id2(1:60*12000)
+
     this%callback => callback
     nsynced=0
     ndecoded=0

lib/lorentzian.f90

@@ -8,15 +8,15 @@ subroutine lorentzian(y,npts,a)
 !         a(4) = width
 
   real y(npts)
-  real a(4)
+  real a(5)
   real deltaa(4)
-  real x(2000)
 
   if(npts.gt.2000) stop 'Error in lorentzian'
 
   a=0.
   df=12000.0/8192.0                               !df = 1.465 Hz
   width=0.
+  ipk=0
   ymax=-1.e30
   do i=1,npts
      if(y(i).gt.ymax) then
@@ -93,12 +93,10 @@ subroutine lorentzian(y,npts,a)
      enddo
      chisqr=fchisq0(y,npts,a)
 !       write(*,4000) 0,0,a,chisqr
-     if(chisqr/chisqr0.gt.0.99) go to 30
+     if(chisqr/chisqr0.gt.0.99) exit
      chisqr0=chisqr
   enddo
-    
+  a(5)=chisqr
-30 ccfbest=ccfmax * (1378.125/fsample)**2
-  dtbest=dtmax
 
   return
 end subroutine lorentzian