Much improved detection of sync in JT4 decoder.

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

CMakeLists.txt

@@ -327,6 +327,7 @@ set (wsjt_FSRCS
   lib/filbig.f90
   lib/flat1.f90
   lib/flat1a.f90
+  lib/flat1b.f90
   lib/flat2.f90
   lib/flat4.f90
   lib/flat65.f90
@@ -352,7 +353,6 @@ set (wsjt_FSRCS
   lib/hashing.f90
   lib/hint65.f90
   lib/hspec.f90
-  lib/image.f90
   lib/indexx.f90
   lib/init_random_seed.f90
   lib/interleave4.f90
@@ -430,7 +430,6 @@ set (wsjt_FSRCS
   lib/wavhdr.f90
   lib/xcor.f90
   lib/xcor4.f90
-  lib/zplt.f90
   lib/wavhdr.f90
   lib/wqencode.f90
   lib/wspr_downsample.f90

lib/decoder.f90

@@ -162,9 +162,8 @@ contains
 
     if (have_sync) then
        decoded=decoded0
-!       write(*,3001) 'A',is_deep,is_average,int(qual),ave,decoded
+       cflags='   '
-!3001   format(a1,2L2,2i4,1x,a22)
+       if(decoded.ne.'                      ') cflags='f  '
-       cflags='f  '
        if(is_deep) then
           cflags(1:2)='d1'
           write(cflags(3:3),'(i1)') min(int(qual),9)

lib/flat1b.f90

@@ -0,0 +1,29 @@
+subroutine flat1b(psavg,nsmo,s2,nh,nsteps,nhmax,nsmax)
+
+  real psavg(nh)
+  real s2(nhmax,nsmax)
+  real x(8192)
+
+  ia=nsmo/2 + 1
+  ib=nh - nsmo/2 - 1
+  do i=ia,ib
+     call pctile(psavg(i-nsmo/2),nsmo,50,x(i))
+  enddo
+  do i=1,ia-1
+     x(i)=x(ia)
+  enddo
+  do i=ib+1,nh
+     x(i)=x(ib)
+  enddo
+
+  do i=1,nh
+     psavg(i)=psavg(i)/x(i)
+     do j=1,nsteps
+        s2(i,j)=s2(i,j)/x(i)
+     enddo
+  enddo
+
+  return
+end subroutine flat1b
+
+      

lib/jt4_decode.f90

@@ -110,9 +110,13 @@ contains
     logical, intent(in) :: NAgain,NClearAve
     character(len=12), intent(in) :: mycall,hiscall
     character(len=6), intent(in) :: hisgrid
-
     real, intent(in) :: dat(npts) !Raw data
-    real z(458,65)
+
+    real ccfblue(-5:540)                             !CCF in time
+    real ccfred(-224:224)                            !CCF in frequency
+    real ps0(450)
+
+!    real z(458,65)
     logical first,prtavg
     character decoded*22,special*5
     character*22 avemsg,deepmsg,deepave,blank,deepmsg0,deepave1
@@ -129,7 +133,8 @@ contains
     endif
 
     zz=0.
-    syncmin=3.0 + minsync
+!    syncmin=3.0 + minsync
+    syncmin=1.0+minsync
     naggressive=0
     if(ndepth.ge.2) naggressive=1
     nq1=3
@@ -150,29 +155,14 @@ contains
 
 ! Attempt to synchronize: look for sync pattern, get DF and DT.
     call timer('sync4   ',0)
-    call sync4(dat,npts,mode4,minw)
+    mousedf=nint(nfqso + 1.5*4.375*mode4 - 1270.46)
+    call sync4(dat,npts,ntol,1,MouseDF,4,mode4,minw+1,dtx,dfx,    &
+         snrx,snrsync,ccfblue,ccfred,flip,width,ps0)
+    sync=snrsync
+    dtxz=dtx-0.8
+    nfreqz=dfx + 1270.46 - 1.5*4.375*mode4
     call timer('sync4   ',1)
 
-    call timer('zplt    ',0)
-    do ich=1,7
-       z(1:458,1:65)=zz(274:731,1:65,ich)
-       call zplt(z,ich-4,syncz,dtxz,nfreqz,flipz,sync2z,0,emedelay,dttol,     &
-            nfqso,ntol)
-    enddo
-    call timer('zplt    ',1)
-
-! Use results from zplt
-!### NB: JT4 is severely "sync limited" at present...  (Maybe not still true???)
-
-!###  TESTS ONLY! ###
-    nfreqz=1000
-    dtxz=0.0
-    flipz=1.0
-    syncz=5.0
-!###
-
-    flip=flipz
-    sync=syncz
     snrx=db(sync) - 26.
     nsnr=nint(snrx)
     if(sync.lt.syncmin) then
@@ -368,7 +358,7 @@ contains
        if(flipsave(i).lt.0.0) csync='#'
        if (associated (this%average_callback)) then
           call this%average_callback(cused(i) .eq. '$',iutc(i),               &
-               syncsave(i) - 5.,dtsave(i),nfsave(i),flipsave(i) .lt.0.)
+               syncsave(i),dtsave(i),nfsave(i),flipsave(i).lt.0.)
        end if
     enddo
 

lib/sync4.f90

@@ -1,22 +1,33 @@
-subroutine sync4(dat,jz,mode4,minw)
+subroutine sync4(dat,jz,ntol,NFreeze,MouseDF,mode,mode4,minwidth,    &
+     dtx,dfx,snrx,snrsync,ccfblue,ccfred1,flip,width,ps0)
 
 ! Synchronizes JT4 data, finding the best-fit DT and DF.  
 
-  use jt4
-  use timer_module, only: timer
-
   parameter (NFFTMAX=2520)         !Max length of FFTs
   parameter (NHMAX=NFFTMAX/2)      !Max length of power spectra
   parameter (NSMAX=525)            !Max number of half-symbol steps
+  integer ntol                     !Range of DF search
   real dat(jz)
   real psavg(NHMAX)                !Average spectrum of whole record
+  real ps0(450)                    !Avg spectrum for plotting
   real s2(NHMAX,NSMAX)             !2d spectrum, stepped by half-symbols
+  real ccfblue(-5:540)             !CCF with pseudorandom sequence
+  real ccfred(-450:450)            !Peak of ccfblue, as function of freq
+  real red(-450:450)               !Peak of ccfblue, as function of freq
+  real ccfred1(-224:224)           !Peak of ccfblue, as function of freq
   real tmp(1260)
+  integer ipk1(1)
+  integer nch(7)
+  logical savered
+  equivalence (ipk1,ipk1a)
+  data nch/1,2,4,9,18,36,72/
   save
 
+!  write(*,3001) 'A',ntol,nfreeze,mousedf,mode,mode4,minwidth
+!3001 format(a1,6i6)
+
 ! Do FFTs of twice symbol length, stepped by half symbols.  Note that 
 ! we have already downsampled the data by factor of 2.
-
   nsym=207
   nfft=2520
   nh=nfft/2
@@ -24,38 +35,145 @@ subroutine sync4(dat,jz,mode4,minw)
   nsteps=jz/nq - 1
   df=0.5*11025.0/nfft
   psavg(1:nh)=0.
+  if(mode.eq.-999) width=0.                        !Silence compiler warning
 
-  call timer('ps4     ',0)
   do j=1,nsteps                     !Compute spectrum for each step, get average
      k=(j-1)*nq + 1
      call ps4(dat(k),nfft,s2(1,j))
      psavg(1:nh)=psavg(1:nh) + s2(1:nh,j)
   enddo
-  call timer('ps4     ',1)
 
-  call timer('flat1a  ',0)
   nsmo=min(10*mode4,150)
-  call flat1a(psavg,nsmo,s2,nh,nsteps,NHMAX,NSMAX)        !Flatten spectra
+  call flat1b(psavg,nsmo,s2,nh,nsteps,NHMAX,NSMAX)        !Flatten spectra
-  call timer('flat1a  ',1)
 
-  call timer('smo     ',0)
   if(mode4.ge.9) call smo(psavg,nh,tmp,mode4/4)
-  call timer('smo     ',1)
 
-  ia=600.0/df
+  i0=132
-  ib=1600.0/df
+  do i=1,450
+     ps0(i)=5.0*(psavg(i0+2*i) + psavg(i0+2*i+1) - 2.0)
+  enddo
 
-!  ichmax=1.0+log(float(mode4))/log(2.0)
+! Set freq and lag ranges
-  do ich=minw+1,7                     !Find best width
+  famin=200.0 + 3*mode4*df
+  fbmax=2700.0 - 3*mode4*df
+  fa=famin
+  fb=fbmax
+  if(NFreeze.eq.1) then
+     fa=max(famin,1270.46+MouseDF-ntol)
+     fb=min(fbmax,1270.46+MouseDF+ntol)
+  else
+     fa=max(famin,1270.46+MouseDF-600)
+     fb=min(fbmax,1270.46+MouseDF+600)
+  endif
+  ia=fa/df - 3*mode4                   !Index of lowest tone, bottom of range
+  ib=fb/df - 3*mode4                   !Index of lowest tone, top of range
+  i0=nint(1270.46/df)
+  irange=450
+  if(ia-i0.lt.-irange) ia=i0-irange
+  if(ib-i0.gt.irange)  ib=i0+irange
+  lag1=-5
+  lag2=59
+  syncbest=-1.e30
+  ccfred=0.
+  jmax=-1000
+  jmin=1000
+
+  do ich=minwidth,7                       !Find best width
      kz=nch(ich)/2
-! Set istep>1 for wide submodes?
+     savered=.false.
      do i=ia+kz,ib-kz                     !Find best frequency channel for CCF
-        call timer('xcor4   ',0)
+        call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0,   &
-        call xcor4(s2,i,nsteps,nsym,ich,mode4)
+             lagpk0,flip)
-        call timer('xcor4   ',1)
+        j=i-i0 + 3*mode4
+        if(j.ge.-372 .and. j.le.372) then
+           ccfred(j)=ccf0
+           jmax=max(j,jmax)
+           jmin=min(j,jmin)
+        endif
+
+! Find rms of the CCF, without main peak
+        call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0)
+        sync=abs(ccfblue(lagpk0))
+
+! Find best sync value
+        if(sync.gt.syncbest) then
+           ipk=i
+           lagpk=lagpk0
+           ichpk=ich
+           syncbest=sync
+           savered=.true.
+        endif
      enddo
+     if(savered) red=ccfred
   enddo
 
+  ccfred=red
+!  width=df*nch(ichpk)
+  dfx=(ipk-i0 + 3*mode4)*df
+
+! Peak up in time, at best whole-channel frequency
+  call xcor4(s2,ipk,nsteps,nsym,lag1,lag2,ichpk,mode4,ccfblue,ccfmax,   &
+       lagpk,flip)
+  xlag=lagpk
+  if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then
+     call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
+     xlag=lagpk+dx2
+  endif
+
+! Find rms of the CCF, without the main peak
+  call slope(ccfblue(lag1),lag2-lag1+1,xlag-lag1+1.0)
+  sq=0.
+  nsq=0
+  do lag=lag1,lag2
+     if(abs(lag-xlag).gt.2.0) then
+        sq=sq+ccfblue(lag)**2
+        nsq=nsq+1
+     endif
+  enddo
+  rms=sqrt(sq/nsq)
+  snrsync=max(0.0,db(abs(ccfblue(lagpk)/rms - 1.0)) - 4.5)
+  snrx=-26.
+  if(mode4.eq.2)  snrx=-25.
+  if(mode4.eq.4)  snrx=-24.
+  if(mode4.eq.9)  snrx=-23.
+  if(mode4.eq.18) snrx=-22.
+  if(mode4.eq.36) snrx=-21.
+  if(mode4.eq.72) snrx=-20.
+  snrx=snrx + snrsync
+
+  dt=2.0/11025.0
+  istart=xlag*nq
+  dtx=istart*dt
+  ccfred1=0.
+  jmin=max(jmin,-224)
+  jmax=min(jmax,224)
+  do i=jmin,jmax
+     ccfred1(i)=ccfred(i)
+  enddo
+
+  ipk1=maxloc(ccfred1) - 225
+  ns=0
+  s=0.
+  iw=min(mode4,(ib-ia)/4)
+  do i=jmin,jmax
+     if(abs(i-ipk1a).gt.iw) then
+        s=s+ccfred1(i)
+        ns=ns+1
+     endif
+  enddo
+  base=s/ns
+  ccfred1=ccfred1-base
+  ccf10=0.5*maxval(ccfred1)
+  do i=ipk1a,jmin,-1
+     if(ccfred1(i).le.ccf10) exit
+  enddo
+  i1=i
+  do i=ipk1a,jmax
+     if(ccfred1(i).le.ccf10) exit
+  enddo
+  width=(i-i1)*df
+
   return
 end subroutine sync4
 
+include 'flat1b.f90'

lib/xcor4.f90

@@ -1,18 +1,44 @@
-subroutine xcor4(s2,ipk,nsteps,nsym,ich,mode4)
+subroutine xcor4(s2,ipk,nsteps,nsym,lag1,lag2,ich,mode4,ccf,ccf0,   &
+     lagpk,flip)
 
-! Computes ccf of the 4-FSK spectral array s2 and the pseudo-random 
+! Computes ccf of the 4_FSK spectral array s2 and the pseudo-random 
 ! array pr2.  Returns peak of CCF and the lag at which peak occurs.  
 ! The CCF peak may be either positive or negative, with negative
-! implying a message with report.
+! implying the "OOO" message.
 
-  use jt4
   parameter (NHMAX=1260)           !Max length of power spectra
   parameter (NSMAX=525)            !Max number of half-symbol steps
   real s2(NHMAX,NSMAX)             !2d spectrum, stepped by half-symbols
   real a(NSMAX)
+  real ccf(-5:540)
+  integer nch(7)
+  integer npr2(207)
+  real pr2(207)
+  logical first
+  data lagmin/0/                    !Silence compiler warning
+  data first/.true./
+  data npr2/                                                        &
+       0,0,0,0,1,1,0,0,0,1,1,0,1,1,0,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0, &
+       0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,1,1,0,1,0,1,1,1,1,1,0,1,0,0,0, &
+       1,0,0,1,0,0,1,1,1,1,1,0,0,0,1,0,1,0,0,0,1,1,1,1,0,1,1,0,0,1, &
+       0,0,0,1,1,0,1,0,1,0,1,0,1,0,1,1,1,1,1,0,1,0,1,0,1,1,0,1,0,1, &
+       0,1,1,1,0,0,1,0,1,1,0,1,1,1,1,0,0,0,0,1,1,0,1,1,0,0,0,1,1,1, &
+       0,1,1,1,0,1,1,1,0,0,1,0,0,0,1,1,0,1,1,0,0,1,0,0,0,1,1,1,1,1, &
+       1,0,0,1,1,0,0,0,0,1,1,0,0,0,1,0,1,1,0,1,1,1,1,0,1,0,1/
+  data nch/1,2,4,9,18,36,72/
   save
 
+  if(first) then
+     do i=1,207
+        pr2(i)=2*npr2(i)-1
+     enddo
+     first=.false.
+  endif
+
+  ccfmax=0.
+  ccfmin=0.
   nw=nch(ich)
+
   do j=1,nsteps
      n=2*mode4
      if(mode4.eq.1) then
@@ -36,14 +62,34 @@ subroutine xcor4(s2,ipk,nsteps,nsym,ich,mode4)
      endif
   enddo
 
-  do lag=1,65
+  do lag=lag1,lag2
      x=0.
      do i=1,nsym
         j=2*i-1+lag
-        if(j.ge.1 .and. j.le.nsteps) x=x+a(j)*float(2*npr(i)-1)
+        if(j.ge.1 .and. j.le.nsteps) x=x+a(j)*pr2(i)
      enddo
-     zz(ipk,lag,ich)=x
+     ccf(lag)=2*x                        !The 2 is for plotting scale
+     if(ccf(lag).gt.ccfmax) then
+        ccfmax=ccf(lag)
+        lagpk=lag
+     endif
+
+     if(ccf(lag).lt.ccfmin) then
+        ccfmin=ccf(lag)
+        lagmin=lag
+     endif
   enddo
 
+  ccf0=ccfmax
+  flip=1.0
+  if(-ccfmin.gt.ccfmax) then
+     do lag=lag1,lag2
+        ccf(lag)=-ccf(lag)
+     enddo
+     lagpk=lagmin
+     ccf0=-ccfmin
+     flip=-1.0
+  endif
+
   return
 end subroutine xcor4