Starting to implement decoding of JT65 shorthand messages. Needs testing!

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

CMakeLists.txt

@@ -393,6 +393,8 @@ set (wsjt_FSRCS
   lib/savec2.f90
   lib/sec_midn.f90
   lib/setup65.f90
+  lib/sh65.f90
+  lib/sh65snr.f90
   lib/slasubs.f
   lib/sleep_msec.f90
   lib/slope.f90

lib/Makefile.msk

@@ -54,9 +54,13 @@ OBJS1 = t1.o four2a.o db.o
 t1:   $(OBJS1)
 	$(FC) -o t1 $(OBJS1) -lfftw3f
 
-OBJS2 = t6.o four2a.o db.o
-t6:   $(OBJS2)
-	$(FC) -o t6 $(OBJS2) -lfftw3f
+OBJS2 = t2.o four2a.o db.o
+t2:   $(OBJS2)
+	$(FC) -o t2 $(OBJS2) -lfftw3f
+
+OBJS6 = t6.o four2a.o db.o
+t6:   $(OBJS6)
+	$(FC) -o t6 $(OBJS6) -lfftw3f
 
 nhash.o: wsprd/nhash.h wsprd/nhash.c
 	 $(CC) -c -O2 wsprd/nhash.c

lib/Makefile.mskWin

@@ -61,9 +61,13 @@ OBJS1 = fixwav.o wavhdr.o
 fixwav.exe: $(OBJS1)
 	$(FC) -o fixwav.exe $(OBJS1)
 
-OBJS2 = t6.o four2a.o db.o
-t6:   $(OBJS2)
-	$(FC) -o t6 $(OBJS2) C:\JTSDK\fftw3f\libfftw3f-3.dll
+OBJS2 = t2.o four2a.o db.o
+t2:   $(OBJS2)
+	$(FC) -o t2 $(OBJS2) C:\JTSDK\fftw3f\libfftw3f-3.dll
+
+OBJS6 = t6.o four2a.o db.o
+t6:   $(OBJS6)
+	$(FC) -o t6 $(OBJS6) C:\JTSDK\fftw3f\libfftw3f-3.dll
 
 .PHONY : clean
 

lib/decode65a.f90

@@ -1,6 +1,6 @@
 subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
-     naggressive,ndepth,mycall,hiscall,hisgrid,nexp_decode,sync2,    &
-     a,dt,nft,qual,nhist,nsmo,decoded)
+     naggressive,ndepth,ntol,mycall,hiscall,hisgrid,nexp_decode,     &
+     single_decode,sync2,a,dt,nft,qual,nhist,nsmo,decoded)
 
 ! Apply AFC corrections to a candidate JT65 signal, then decode it.
 
@@ -15,7 +15,7 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
   complex c5a(512)
   real s2(66,126)
   real a(5)
-  logical first
+  logical single_decode,first
   character decoded*22,decoded_best*22
   character mycall*12,hiscall*12,hisgrid*6
   data first/.true./,jjjmin/1000/,jjjmax/-1000/
@@ -27,9 +27,23 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
   call filbig(dd,npts,f0,newdat,cx,n5,sq0)
   if(mode65.eq.4) call filbig(dd,npts,f0+355.297852,newdat,cx1,n5,sq0)
   call timer('filbig  ',1)
-
 ! NB: cx has sample rate 12000*77125/672000 = 1378.125 Hz
 
+! Check for a shorthand message
+  if(single_decode) then
+     call sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
+     if(nspecial.gt.0) then
+        a=0.
+        a(1)=xdf
+        if(nspecial.eq.2) decoded='RO'
+        if(nspecial.eq.3) decoded='RRR'
+        if(nspecial.eq.4) decoded='73'
+        nflip=0
+        sync2=snrdb
+        go to 900
+     endif
+  endif
+
 ! Find best DF, drift, curvature, and DT.  Start by downsampling to 344.53125 Hz
   call timer('fil6521 ',0)
   call fil6521(cx,n5,c5x,n6)
@@ -144,5 +158,5 @@ subroutine decode65a(dd,npts,newdat,nqd,f0,nflip,mode65,ntrials,     &
 
   call timer('dec65b  ',1)
 
-  return
+900 return
 end subroutine decode65a

lib/decoder.f90

@@ -257,7 +257,7 @@ contains
           endif
        endif
        csync=' '
-       if(sync.ge.float(minsync)) then
+       if(nflip.ne.0 .and. sync.ge.max(0.0,float(minsync))) then
           csync='*'
           if(nflip.eq.-1) then
              csync='#'

lib/jt65_decode.f90

@@ -190,15 +190,17 @@ contains
           if(single_decode) then
              flip=ca(icand)%flip
              nflip=flip
-             if(sync1.lt.float(minsync)) cycle
           endif
           if(ipass.eq.1) ntry65a=ntry65a + 1
           if(ipass.eq.2) ntry65b=ntry65b + 1
           call timer('decod65a',0)
           call decode65a(dd,npts,first_time,nqd,freq,nflip,mode65,nvec,     &
-               naggressive,ndepth,mycall,hiscall,hisgrid,nexp_decode,       &
-               sync2,a,dtx,nft,qual,nhist,nsmo,decoded)
+               naggressive,ndepth,ntol,mycall,hiscall,hisgrid,              &
+               nexp_decode,single_decode,sync2,a,dtx,nft,qual,nhist,        &
+               nsmo,decoded)
           call timer('decod65a',1)
+          if(sync1.lt.float(minsync) .and.                                  &
+               decoded.eq.'                      ') nflip=0
           if(nft.ne.0) nsum=1
 
 !          ncandidates=param(0)

lib/sh65.f90

@@ -0,0 +1,110 @@
+subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
+  parameter(NFFT=2048,NH=NFFT/2,MAXSTEPS=150)
+  complex cx(90000)
+  complex c(0:NFFT-1)
+  real s(-NH+1:NH)
+  real s2(-NH+1:NH,MAXSTEPS)
+  real ss(-NH+1:NH,8)
+  real sigmax(8)
+  integer ipk(8)
+
+  s=0.
+  ss=0.
+
+  jstep=NFFT/4
+  nblks=n5/jstep - 3
+  ia=-jstep+1
+  do iblk=1,nblks
+     ia=ia+jstep
+     ib=ia+NFFT-1
+     c=cx(ia:ib)
+     call four2a(c,nfft,1,1,1)            !c2c FFT
+     do i=0,NFFT-1
+        j=i
+        if(j.gt.NH) j=j-NFFT
+        p=real(c(i))**2 + aimag(c(i))**2
+        s(j)=s(j) + p
+        s2(j,iblk)=p
+     enddo
+     n=mod(iblk-1,8) +1
+     ss(-NH+1:NH,n)=ss(-NH+1:NH,n) + s2(-NH+1:NH,iblk)
+  enddo
+
+  s=1.e-6*s
+  ss=1.e-6*ss
+  df=1378.1285/NFFT
+  do i=-NH+1,NH
+     f=i*df
+     write(13,1010) f,s(i),ss(i,1:8)
+1010 format(10f10.3)
+  enddo
+
+  nfac=40*mode65
+  dtstep=0.25/df
+
+! Define range of frequencies to be searched
+  fa=-ntol
+  fb=ntol
+  ia2=max(-NH+1,nint(fa/df))
+! Upper tone is above sync tone by 4*nfac*df Hz
+  ib2=min(NH,nint(fb/df + 4.1*nfac))  
+
+! Find strongest line in each of the 4 phases, repeating for each drift rate.
+  sbest=0.
+  snrbest=0.
+  nbest=1
+  ipk=0
+
+  do n=1,8
+     sigmax(n)=0.
+     do i=ia2,ib2
+        sig=ss(i,n)
+        if(sig.ge.sigmax(n)) then
+           ipk(n)=i
+           sigmax(n)=sig
+           if(sig.ge.sbest) then
+              sbest=sig
+              nbest=n
+           endif
+        endif
+     enddo
+  enddo
+  n2best=nbest+4
+  if(n2best.gt.8) n2best=nbest-4
+  xdf=min(ipk(nbest),ipk(n2best))*df
+  nspecial=0
+  if(abs(xdf).gt.ntol) go to 10
+
+  idiff=abs(ipk(nbest)-ipk(n2best))
+  xk=float(idiff)/nfac
+  k=nint(xk)
+  iderr=nint((xk-k)*nfac)
+  maxerr=nint(0.008*abs(idiff) + 0.51)
+  if(abs(iderr).le.maxerr .and. k.ge.2 .and. k.le.4) nspecial=k
+  nstest=0
+  if(nspecial.gt.0) then
+     call sh65snr(ss(ia2,nbest),ib2-ia2+1,snr1)
+     call sh65snr(ss(ia2,n2best),ib2-ia2+1,snr2)
+     snr=0.5*(snr1+snr2)
+     if(snr.gt.snrbest) then
+        snrbest=snr
+        nspecialbest=nspecial
+        nstest=snr/2.0 - 2.0             !Threshold set here
+        if(nstest.lt.0) nstest=0
+        if(nstest.gt.10) nstest=10
+        dfsh=nint(xdf)
+        iderrbest=iderr
+        snrdb=db(snr) - db(2500.0/df) - db(sqrt(nblks/4.0))+1.8
+        n1=nbest
+        n2=n2best
+        ipk1=ipk(n1)
+        ipk2=ipk(n2)
+     endif
+  endif
+  if(nstest.eq.0) nspecial=0
+10 continue
+
+!  print*,'a',ia2,ib2,snrdb,xdf,nspecial
+
+  return
+end subroutine sh65

lib/sh65snr.f90

@@ -0,0 +1,36 @@
+subroutine sh65snr(x,nz,snr)
+
+  real x(nz)
+
+  ipk=0 !Shut up compiler warnings. -db
+  smax=-1.e30
+  do i=1,nz
+     if(x(i).gt.smax) then
+        ipk=i
+        smax=x(i)
+     endif
+     s=s+x(i)
+  enddo
+
+  s=0.
+  ns=0
+  do i=1,nz
+     if(abs(i-ipk).ge.3) then
+        s=s+x(i)
+        ns=ns+1
+     endif
+  enddo
+  ave=s/ns
+
+  sq=0.
+  do i=1,nz
+     if(abs(i-ipk).ge.3) then
+        sq=sq+(x(i)-ave)**2
+        ns=ns+1
+     endif
+  enddo
+  rms=sqrt(sq/(nz-2))
+  snr=(smax-ave)/rms
+
+  return
+end subroutine sh65snr