Implement "matched filter on s3(64,63)" as the criterion for hinted

(aka experience-based) decoding.


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

lib/decode65a.f90

@@ -1,5 +1,5 @@
 subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
-     naggressive,ndepth,nexp_decode,sync2,a,dt,nsf,nhist,decoded)
+     naggressive,ndepth,nexp_decode,sync2,a,dt,nft,qual,nhist,decoded)
 
 ! Apply AFC corrections to a candidate JT65 signal, then decode it.
 
@@ -74,7 +74,7 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
 
   call timer('dec65b  ',0)
   call decode65b(s2,nflip,mode65,ntrials,naggressive,ndepth,nexp_decode,   &
-       nqd,nsf,nhist,decoded)
+       nqd,nft,qual,nhist,decoded)
   dt=dtbest !return new, improved estimate of dt
   call timer('dec65b  ',1)
 

lib/decode65b.f90

@@ -1,5 +1,5 @@
 subroutine decode65b(s2,nflip,mode65,ntrials,naggressive,ndepth,nexp_decode, &
-     nqd,nsf,nhist,decoded)
+     nqd,nft,qual,nhist,decoded)
 
   real s2(66,126)
   real s3(64,63)
@@ -23,15 +23,15 @@ subroutine decode65b(s2,nflip,mode65,ntrials,naggressive,ndepth,nexp_decode, &
   enddo
 
   nadd=mode65
-  call extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,     &
+  call extract(s3,nadd,nqd,ntrials,naggressive,ndepth,                  &
-       ncount,nhist,decoded,ltext,nsf)           !Extract the message
+       ncount,nhist,decoded,ltext,nft,qual)           !Extract the message
 
 ! Suppress "birdie messages" and other garbage decodes:
   if(decoded(1:7).eq.'000AAA ') ncount=-1
   if(decoded(1:7).eq.'0L6MWK ') ncount=-1
   if(nflip.lt.0 .and. ltext) ncount=-1
   if(ncount.lt.0) then 
-     nsf=0
+     nft=0
      decoded='                      '
   endif
 

lib/demod64a.f90

@@ -38,6 +38,7 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
            i1=i                              !Most reliable
         endif
      enddo
+     if(psum.eq.0.0) psum=1.e-6
 
      s2=-1.e30
      do i=1,64
@@ -47,7 +48,7 @@ subroutine demod64a(s3,nadd,afac1,mrsym,mrprob,mr2sym,mr2prob,ntest,nlow)
         endif
      enddo
      p1=s1/psum                !Use these for sfrsd
-     p2=s2/psum                !
+     p2=s2/psum                !...
      mrsym(j)=i1-1
      mr2sym(j)=i2-1
      mrprob(j)=scale*p1

lib/exp_decode65.f90

@@ -1,29 +1,27 @@
-subroutine exp_decode65(mrsym,mrprob,mr2sym,nexp_decode,nhard,nsoft,nbest,  &
+subroutine exp_decode65(s3,mrs,mrs2,mode65,flip,mycall,qual,decoded)
-     correct)
 
   use packjt
   use prog_args
   parameter (NMAX=10000)
+  real s3(64,63)
+  real pp(NMAX)
   integer*1 sym1(0:62,NMAX)
-  integer mrsym(0:62),mr2sym(0:62),mrprob(0:62)
+  integer*1 sym2(0:62,NMAX)
+  integer mrs(63),mrs2(63)
   integer dgen(12),sym(0:62),sym_rev(0:62)
-  integer test(0:62)  
-  integer correct(0:62)
   character*6 mycall,hiscall(NMAX)
   character*4 hisgrid(NMAX)
   character callsign*12,grid*4
   character*180 line
   character ceme*3,msg*22
+  character*22 msg0(1000),decoded
   logical*1 eme(NMAX)
   logical first
   data first/.true./
-  save first,sym1,nused
+  save first,sym1,nused,msg0,sym2
 
-!###
-  mycall='VK7MO'                   !### TEMPORARY ###
   if(first) then
      neme=1
-
      open(23,file=trim(data_dir)//'/CALL3.TXT',status='unknown')
      icall=0
      j=0
@@ -55,42 +53,84 @@ subroutine exp_decode65(mrsym,mrprob,mr2sym,nexp_decode,nhard,nsoft,nbest,  &
      j=0
      do i=1,ncalls
         if(neme.eq.1 .and. (.not.eme(i))) cycle
-        j=j+1
+
-        msg=mycall//' '//hiscall(i)//' '//hisgrid(i)
+!### Special for tests ###
-        call fmtmsg(msg,iz)
+        do isnr=-20,-30,-1
-        call packmsg(msg,dgen,itype)            !Pack message into 72 bits
+           j=j+1
-        call rs_encode(dgen,sym_rev)            !RS encode
+           msg=mycall//' '//hiscall(i)//' '//hisgrid(i)
-        sym(0:62)=sym_rev(62:0:-1)
+           write(msg(14:18),"(i3,'  ')") isnr
-        sym1(0:62,j)=sym
+           call fmtmsg(msg,iz)
+           call packmsg(msg,dgen,itype)            !Pack message into 72 bits
+           call rs_encode(dgen,sym_rev)            !RS encode
+           sym(0:62)=sym_rev(62:0:-1)
+           sym1(0:62,j)=sym
+
+           call interleave63(sym_rev,1)            !Interleave channel symbols
+           call graycode(sym_rev,63,1,sym_rev)     !Apply Gray code
+           sym2(0:62,j)=sym_rev(0:62)
+           msg0(j)=msg
+        enddo
+
      enddo
      nused=j
      first=.false.
   endif
-!###
 
-  nbest=999999
+  ref0=0.
-  do j=1,nused
+  do j=1,63
-     test=sym1(0:62,j)
+     ref0=ref0 + s3(mrs(j)+1,j)
-     nh=0
-     ns=0
-     do i=0,62
-        k=62-i
-        if(mrsym(k).ne.test(i)) then
-           nh=nh+1
-           if(mr2sym(k).ne.test(i)) ns=ns+mrprob(k)
-        endif
-     enddo
-     ns=63*ns/sum(mrprob)
-
-     if(nh+ns.lt.nbest) then
-        nhard=nh
-        nsoft=ns
-        nbest=nhard+nsoft
-        ncandidates=0
-        ntry=0
-        correct=test
-     endif
   enddo
 
+  p1=-1.e30
+  p2=-1.e30
+  ip1=1                                    !Silence compiler warning
+  do k=1,nused
+     pp(k)=0.
+     if(k.ge.2 .and. k.le.64 .and. flip.lt.0.0) cycle
+! Test all messages if flip=+1; skip the CQ messages if flip=-1.
+     if(flip.gt.0.0 .or. msg0(k)(1:3).ne.'CQ ') then
+        sum=0.
+        ref=ref0
+        do j=1,63
+!           i=ncode(j,k)+1
+           i=sym2(j-1,k)+1
+           sum=sum + s3(i,j)
+           if(i.eq.mrs(j)+1) ref=ref - s3(i,j) + s3(mrs2(j)+1,j)
+        enddo
+        p=sum/ref
+        pp(k)=p
+        if(p.gt.p1) then
+           p1=p
+           ip1=k
+        endif
+     endif
+!     write(*,3001) k,p,msg0(k)
+!3001 format(i5,f10.3,2x,a22)
+  enddo
+
+  do i=1,nused
+     if(pp(i).gt.p2 .and. pp(i).ne.p1) p2=pp(i)
+  enddo
+
+! ### DO NOT REMOVE ### 
+!  call cs_lock('deep65')
+  rewind 77
+  write(77,*) p1,p2
+  call flush(77)
+!  call cs_unlock
+! ### Works OK without it (in both Windows and Linux) if compiled 
+! ### without optimization.  However, in Windows this is a colossal 
+! ### pain because of the way F2PY wants to run the compile step.
+
+  bias=max(1.12*p2,0.335)
+  if(mode65.eq.2) bias=max(1.08*p2,0.405)
+  if(mode65.ge.4) bias=max(1.04*p2,0.505)
+
+  if(p2.eq.p1 .and. p1.ne.-1.e30) stop 'Error in deep65'
+  qual=100.0*(p1-bias)
+  if(qual.lt.0.0) qual=0.0
+  decoded='                      '
+  if(qual.ge.1.0) decoded=msg0(ip1)
+
   return
 end subroutine exp_decode65

lib/extract.f90

@@ -1,5 +1,5 @@
-subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,   &
+subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,              &
-     ncount,nhist,decoded,ltext,nsf)
+     ncount,nhist,decoded,ltext,nft,qual)
 
 ! Input:
 !   s3       64-point spectra for each of 63 data symbols
@@ -11,28 +11,32 @@ subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,   &
 !   nhist    maximum number of identical symbol values
 !   decoded  decoded message (if ncount >=0)
 !   ltext    true if decoded message is free text
-!   nsf    0=no decode; 1=BM decode; 2=KV decode
+!   nft      0=no decode; 1=FT decode; 2=hinted decode
 
   use prog_args                       !shm_key, exe_dir, data_dir
   use packjt
 
   real s3(64,63)
   character decoded*22
+  character*6 mycall
   integer dat4(12)
   integer mrsym(63),mr2sym(63),mrprob(63),mr2prob(63)
+  integer mrs(63),mrs2(63)
   integer correct(63),tmp(63)
   integer param(0:7)
   integer indx(0:62)
   real*8 tt
   logical nokv,ltext
   common/chansyms65/correct
+  common/test000/param                              !### TEST ONLY ###
   data nokv/.false./,nsec1/0/
   save
 
+  qual=0.
   nbirdie=20
   npct=50
   afac1=1.1
-  nsf=0
+  nft=0
   nfail=0
   decoded='                      '
   call pctile(s3,4032,npct,base)
@@ -54,6 +58,9 @@ subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,   &
      go to 1
   endif
 
+  mrs=mrsym
+  mrs2=mr2sym
+
   call graycode65(mrsym,63,-1)        !Remove gray code 
   call interleave63(mrsym,-1)         !Remove interleaving
   call interleave63(mrprob,-1)
@@ -72,22 +79,35 @@ subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,   &
   nhard=param(1)
   nsoft=param(2)
   nerased=param(3)
+  nsofter=param(4)
+  ntotal=param(5)
 
-  if(nhard.lt.0 .and. ndepth.ge.5) then
+  nhard_max=44
+  nd_a=72 + naggressive
+  if(nhard.le.nhard_max .and. ntotal.le.nd_a) nft=1
+
+!  print*,'AAA',ndepth
+  if(nft.eq.0 .and. ndepth.ge.5) then
+!     print*,'BBB',ndepth
      call timer('exp_deco',0)
-     call exp_decode65(mrsym,mrprob,mr2sym,nexp_decode,nhard,nsoft,nbest,  &
+     mode65=1
-          correct)
+     flip=1.0
-     if(nbest.gt.72+2*naggressive) then
+     mycall='K1ABC'                   !### TEMPORARY ###
-        nhard=-1
+     call exp_decode65(s3,mrs,mrs2,mode65,flip,mycall,qual,decoded)
+     if(qual.ge.1.0) then
+        nft=2
+     else
+        param=0
+        ntry=0
      endif
      call timer('exp_deco',1)
+     go to 900
   endif
 
-  ndone=ndone+1
   ncount=-1
   decoded='                      '
   ltext=.false.
-  if(nhard.ge.0) then
+  if(nft.gt.0) then
     !turn the corrected symbol array into channel symbols for subtraction
     !pass it back to jt65a via common block "chansyms65"
      do i=1,12
@@ -100,14 +120,13 @@ subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode,   &
      call interleave63(correct,63,1)
      call graycode65(correct,63,1)
      call unpackmsg(dat4,decoded)     !Unpack the user message
-!  open(43,file='stats.dat',access='append',status='unknown')
-!  write(43,*) nhard,nsoft,nerased,ntry
-!  close(43)
      ncount=0
      if(iand(dat4(10),8).ne.0) ltext=.true.
-     nsf=1
   endif
 900 continue
+  if(nft.eq.1 .and. nhard.lt.0) decoded='                      '
+!  write(*,3300) nft,nhard,ntotal,int(qual),decoded
+!3300 format(4i5,2x,a22)
 
   return
 end subroutine extract

lib/fillcom.f90

@@ -23,10 +23,11 @@ subroutine fillcom(nutc0,ndepth0,nrxfreq,mode,tx9,flow,fsplit,fhigh)
   nzhsym=181
   ndepth=ndepth0
   dttol=3.0
+  minsync=-1               !### TEST ONLY
 
   n2pass=1
-  nranera=8
+  nranera=10               !ntrials=100000
-  naggressive=10
+  naggressive=3
   nrobust=0
 
   if (tx9) then

lib/jt65a.f90

@@ -27,6 +27,8 @@ subroutine jt65a(dd0,npts,newdat,nutc,nf1,nf2,nfqso,ntol,nsubmode,   &
   type(decode) dec(30)
   common/decstats/ntry65a,ntry65b,n65a,n65b,num9,numfano
   common/steve/thresh0
+  common/test000/ncandidates,nhard_min,nsoft_min,nera_best,nsofter_min,   &
+       ntotal_min,ntry         !### TEST ONLY ###
   save
 
   dd=dd0
@@ -87,16 +89,13 @@ subroutine jt65a(dd0,npts,newdat,nutc,nf1,nf2,nfqso,ntol,nsubmode,   &
       freq=ca(icand)%freq
       dtx=ca(icand)%dt
       sync1=ca(icand)%sync
-
       if(ipass.eq.1) ntry65a=ntry65a + 1
       if(ipass.eq.2) ntry65b=ntry65b + 1
       call timer('decod65a',0)
       call decode65a(dd,npts,newdat,nqd,freq,nflip,mode65,nvec,     &
-           naggressive,ndepth,nexp_decode,sync2,a,dtx,nsf,nhist,decoded)
+           naggressive,ndepth,nexp_decode,sync2,a,dtx,nft,qual,nhist,decoded)
       call timer('decod65a',1)
-!write(*,*) icand,freq+a(1),dtx,sync1,sync2
+      if(decoded.eq.decoded0 .and. minsync.ge.0) cycle     !Don't display dupes
-      if(decoded.eq.decoded0) cycle            !Don't display dupes
-
       if(decoded.ne.'                      ' .or. minsync.lt.0) then
         if( nsubtract .eq. 1 ) then
            call timer('subtr65 ',0)
@@ -114,9 +113,12 @@ subroutine jt65a(dd0,npts,newdat,nutc,nf1,nf2,nfqso,ntol,nsubmode,   &
 !$omp critical(decode_results) 
         ndupe=0 ! de-dedupe
         do i=1, ndecoded
-          if( decoded==dec(i)%decoded ) ndupe=1
+          if(decoded==dec(i)%decoded) then
+             ndupe=1
+             exit
+          endif
         enddo
-        if(ndupe.ne.1) then
+        if(ndupe.ne.1 .or. minsync.lt.0) then
           if(ipass.eq.1) n65a=n65a + 1
           if(ipass.eq.2) n65b=n65b + 1
           ndecoded=ndecoded+1
@@ -124,14 +126,22 @@ subroutine jt65a(dd0,npts,newdat,nutc,nf1,nf2,nfqso,ntol,nsubmode,   &
           dec(ndecoded)%dt=dtx
           dec(ndecoded)%sync=sync2
           dec(ndecoded)%decoded=decoded
+          nqual=qual
+          if(nqual.gt.10) nqual=10
           write(*,1010) nutc,nsnr,dtx-1.0,nfreq,decoded
 1010      format(i4.4,i4,f5.1,i5,1x,'#',1x,a22)
           write(13,1012) nutc,nint(sync1),nsnr,dtx-1.0,float(nfreq),ndrift,  &
-             decoded,nsf
+             decoded,nft
 1012      format(i4.4,i4,i5,f6.1,f8.0,i4,3x,a22,' JT65',i4)
           call flush(6)
           call flush(13)
+          write(79,3001) nutc,nint(sync1),nsnr,dtx-1.0,nfreq,ncandidates,    &
+               nhard_min,ntotal_min,ntry,naggressive,nft,nqual,decoded
+3001      format(i4.4,i3,i4,f6.2,i5,i6,i3,i4,i8,i3,i2,i3,1x,a22)
+          flush(79)
         endif
+        decoded0=decoded
+        if(decoded0.eq.'                      ') decoded0='*'
 !$omp end critical(decode_results)
       endif
     enddo                                 !candidate loop