From 56a58fe6ab45e5bfc46ba494b7c1bc6a6b3c5b64 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Wed, 9 Dec 2015 21:02:37 +0000 Subject: [PATCH] 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 | 4 +- lib/decode65b.f90 | 8 +-- lib/demod64a.f90 | 3 +- lib/exp_decode65.f90 | 118 +++++++++++++++++++++++++++++-------------- lib/extract.f90 | 49 ++++++++++++------ lib/fillcom.f90 | 5 +- lib/jt65a.f90 | 26 +++++++--- 7 files changed, 142 insertions(+), 71 deletions(-) diff --git a/lib/decode65a.f90 b/lib/decode65a.f90 index 71f1b35b9..11d8a9618 100644 --- a/lib/decode65a.f90 +++ b/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) diff --git a/lib/decode65b.f90 b/lib/decode65b.f90 index f33d07109..4bbcaa78f 100644 --- a/lib/decode65b.f90 +++ b/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, & - ncount,nhist,decoded,ltext,nsf) !Extract the message + call extract(s3,nadd,nqd,ntrials,naggressive,ndepth, & + 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 diff --git a/lib/demod64a.f90 b/lib/demod64a.f90 index e7521d8ab..034b8e3ed 100644 --- a/lib/demod64a.f90 +++ b/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 diff --git a/lib/exp_decode65.f90 b/lib/exp_decode65.f90 index 6de452c19..d17a21262 100644 --- a/lib/exp_decode65.f90 +++ b/lib/exp_decode65.f90 @@ -1,29 +1,27 @@ -subroutine exp_decode65(mrsym,mrprob,mr2sym,nexp_decode,nhard,nsoft,nbest, & - correct) +subroutine exp_decode65(s3,mrs,mrs2,mode65,flip,mycall,qual,decoded) 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) - 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 + +!### Special for tests ### + do isnr=-20,-30,-1 + j=j+1 + msg=mycall//' '//hiscall(i)//' '//hisgrid(i) + write(msg(14:18),"(i3,' ')") isnr + 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 - do j=1,nused - test=sym1(0:62,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 + ref0=0. + do j=1,63 + ref0=ref0 + s3(mrs(j)+1,j) 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 diff --git a/lib/extract.f90 b/lib/extract.f90 index 28881e5ee..3e72d2fcb 100644 --- a/lib/extract.f90 +++ b/lib/extract.f90 @@ -1,5 +1,5 @@ -subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth,nexp_decode, & - ncount,nhist,decoded,ltext,nsf) +subroutine extract(s3,nadd,nqd,ntrials,naggressive,ndepth, & + 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, & - correct) - if(nbest.gt.72+2*naggressive) then - nhard=-1 + mode65=1 + flip=1.0 + mycall='K1ABC' !### TEMPORARY ### + 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 diff --git a/lib/fillcom.f90 b/lib/fillcom.f90 index 4e9b017d5..af131fad6 100644 --- a/lib/fillcom.f90 +++ b/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 - naggressive=10 + nranera=10 !ntrials=100000 + naggressive=3 nrobust=0 if (tx9) then diff --git a/lib/jt65a.f90 b/lib/jt65a.f90 index 24c998c97..aa0e4887e 100644 --- a/lib/jt65a.f90 +++ b/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) cycle !Don't display dupes - + if(decoded.eq.decoded0 .and. minsync.ge.0) 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