From 4bcc4f35a1e4d79076988265124f13585cd90230 Mon Sep 17 00:00:00 2001 From: Joe Taylor Date: Thu, 19 May 2016 19:19:47 +0000 Subject: [PATCH] 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 | 3 +- lib/decoder.f90 | 5 +- lib/flat1b.f90 | 29 ++++++ lib/jt4_decode.f90 | 40 +++----- lib/sync4.f90 | 240 +++++++++++++++++++++++++++++++++------------ lib/xcor4.f90 | 144 ++++++++++++++++++--------- 6 files changed, 321 insertions(+), 140 deletions(-) create mode 100644 lib/flat1b.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 50d8827c7..b9816f0a9 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/lib/decoder.f90 b/lib/decoder.f90 index 48cb1f901..2a3145595 100644 --- a/lib/decoder.f90 +++ b/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 -!3001 format(a1,2L2,2i4,1x,a22) - cflags='f ' + cflags=' ' + if(decoded.ne.' ') cflags='f ' if(is_deep) then cflags(1:2)='d1' write(cflags(3:3),'(i1)') min(int(qual),9) diff --git a/lib/flat1b.f90 b/lib/flat1b.f90 new file mode 100644 index 000000000..bc22b8557 --- /dev/null +++ b/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 + + diff --git a/lib/jt4_decode.f90 b/lib/jt4_decode.f90 index 413e3ca89..a1cb9a0b4 100644 --- a/lib/jt4_decode.f90 +++ b/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 diff --git a/lib/sync4.f90 b/lib/sync4.f90 index dd92ea123..c7f33b470 100644 --- a/lib/sync4.f90 +++ b/lib/sync4.f90 @@ -1,61 +1,179 @@ -subroutine sync4(dat,jz,mode4,minw) - -! 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 - real dat(jz) - real psavg(NHMAX) !Average spectrum of whole record - real s2(NHMAX,NSMAX) !2d spectrum, stepped by half-symbols - real tmp(1260) - save - -! 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 - nq=nfft/4 - nsteps=jz/nq - 1 - df=0.5*11025.0/nfft - psavg(1:nh)=0. - - 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 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 - ib=1600.0/df - -! ichmax=1.0+log(float(mode4))/log(2.0) - do ich=minw+1,7 !Find best width - kz=nch(ich)/2 -! Set istep>1 for wide submodes? - do i=ia+kz,ib-kz !Find best frequency channel for CCF - call timer('xcor4 ',0) - call xcor4(s2,i,nsteps,nsym,ich,mode4) - call timer('xcor4 ',1) - enddo - enddo - - return -end subroutine sync4 - +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. + + 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 + nq=nfft/4 + nsteps=jz/nq - 1 + df=0.5*11025.0/nfft + psavg(1:nh)=0. + if(mode.eq.-999) width=0. !Silence compiler warning + + 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 + + nsmo=min(10*mode4,150) + call flat1b(psavg,nsmo,s2,nh,nsteps,NHMAX,NSMAX) !Flatten spectra + + if(mode4.ge.9) call smo(psavg,nh,tmp,mode4/4) + + i0=132 + do i=1,450 + ps0(i)=5.0*(psavg(i0+2*i) + psavg(i0+2*i+1) - 2.0) + enddo + +! Set freq and lag ranges + 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 + savered=.false. + do i=ia+kz,ib-kz !Find best frequency channel for CCF + call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0, & + lagpk0,flip) + 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' diff --git a/lib/xcor4.f90 b/lib/xcor4.f90 index 296c65cc2..9e1f5a529 100644 --- a/lib/xcor4.f90 +++ b/lib/xcor4.f90 @@ -1,49 +1,95 @@ -subroutine xcor4(s2,ipk,nsteps,nsym,ich,mode4) - -! 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. - - 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) - save - - nw=nch(ich) - do j=1,nsteps - n=2*mode4 - if(mode4.eq.1) then - a(j)=max(s2(ipk+n,j),s2(ipk+3*n,j)) - max(s2(ipk ,j),s2(ipk+2*n,j)) - else - kz=max(1,nw/2) - ss0=0. - ss1=0. - ss2=0. - ss3=0. - wsum=0. - do k=-kz+1,kz-1 - w=float(kz-iabs(k))/nw - wsum=wsum+w - ss0=ss0 + w*s2(ipk +k,j) - ss1=ss1 + w*s2(ipk+ n+k,j) - ss2=ss2 + w*s2(ipk+2*n+k,j) - ss3=ss3 + w*s2(ipk+3*n+k,j) - enddo - a(j)=(max(ss1,ss3) - max(ss0,ss2))/sqrt(wsum) - endif - enddo - - do lag=1,65 - 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) - enddo - zz(ipk,lag,ich)=x - enddo - - return -end subroutine xcor4 +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 +! 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 the "OOO" message. + + 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 + a(j)=max(s2(ipk+n,j),s2(ipk+3*n,j)) - max(s2(ipk ,j),s2(ipk+2*n,j)) + else + kz=max(1,nw/2) + ss0=0. + ss1=0. + ss2=0. + ss3=0. + wsum=0. + do k=-kz+1,kz-1 + w=float(kz-iabs(k))/nw + wsum=wsum+w + ss0=ss0 + w*s2(ipk +k,j) + ss1=ss1 + w*s2(ipk+ n+k,j) + ss2=ss2 + w*s2(ipk+2*n+k,j) + ss3=ss3 + w*s2(ipk+3*n+k,j) + enddo + a(j)=(max(ss1,ss3) - max(ss0,ss2))/sqrt(wsum) + endif + enddo + + 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)*pr2(i) + enddo + 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