From d8f9f385b498d256ef3fbb535b7569a9d5f24bdd Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Wed, 12 Jul 2017 21:49:39 +0000 Subject: [PATCH] Use quarter-symbol steps for time sync. Lower sync threshold. Implement subtraction and two pass decoding. Use osd2 only near nfqso. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7861 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/fsk4hf/bpdecode174.f90 | 2 +- lib/fsk4hf/ft8_params.f90 | 3 +- lib/fsk4hf/ft8b.f90 | 52 ++++++++++++++------------ lib/fsk4hf/osd174.f90 | 38 +++++++++++-------- lib/fsk4hf/subtractft8.f90 | 2 +- lib/fsk4hf/sync8.f90 | 50 ++++++++++++++----------- lib/ft8_decode.f90 | 76 +++++++++++++++++++++----------------- lib/jt9.f90 | 2 +- 8 files changed, 127 insertions(+), 98 deletions(-) diff --git a/lib/fsk4hf/bpdecode174.f90 b/lib/fsk4hf/bpdecode174.f90 index aec0fcd3a..9442cc4a1 100644 --- a/lib/fsk4hf/bpdecode174.f90 +++ b/lib/fsk4hf/bpdecode174.f90 @@ -1,4 +1,4 @@ -subroutine bpdecode174(llr,apmask,maxiterations,decoded,cw,nharderror) +subroutine bpdecode174(llr,apmask,maxiterations,decoded,cw,nharderror,iter) ! ! A log-domain belief propagation decoder for the (174,87) code. ! diff --git a/lib/fsk4hf/ft8_params.f90 b/lib/fsk4hf/ft8_params.f90 index 8441fdb42..af063cc51 100644 --- a/lib/fsk4hf/ft8_params.f90 +++ b/lib/fsk4hf/ft8_params.f90 @@ -7,5 +7,6 @@ parameter (NSPS=1920) !Samples per symbol at 12000 S/s parameter (NZ=NSPS*NN) !Samples in full 15 s waveform (151,680) parameter (NMAX=15*12000) !Samples in iwave (180,000) parameter (NFFT1=2*NSPS, NH1=NFFT1/2) !Length of FFTs for symbol spectra -parameter (NHSYM=2*NMAX/NH1-1) !Number of symbol spectra (1/2-sym steps) +parameter (NSTEP=NSPS/4) !Rough time-sync step size +parameter (NHSYM=NMAX/NSTEP-1) !Number of symbol spectra (1/2-sym steps) parameter (NDOWN=60) !Downsample factor diff --git a/lib/fsk4hf/ft8b.f90 b/lib/fsk4hf/ft8b.f90 index b79840beb..adda05529 100644 --- a/lib/fsk4hf/ft8b.f90 +++ b/lib/fsk4hf/ft8b.f90 @@ -1,5 +1,5 @@ -subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, & - dmin,nbadcrc,message,xsnr) +subroutine ft8b(dd0,newdat,nfqso,ndepth,lsubtract,icand,sync0,f1,xdt,apsym,nharderrors, & + dmin,nbadcrc,iap,ipass,message,xsnr) use timer_module, only: timer include 'ft8_params.f90' @@ -18,12 +18,12 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, complex cd0(3200) complex ctwk(32) complex csymb(32) - logical newdat + logical newdat,lsubtract data rr73/-1,1,1,1,1,1,1,-1,1,1,-1/ - max_iterations=40 + max_iterations=30 norder=2 - if(ndepth.eq.3 .and. abs(nfqso-f1).lt.10.0) norder=3 + nharderrors=-1 fs2=12000.0/NDOWN dt2=1.0/fs2 twopi=8.0*atan(1.0) @@ -36,7 +36,7 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, i0=nint(xdt*fs2) !Initial guess for start of signal smax=0.0 - do idt=i0-16,i0+16 !Search over +/- half a symbol + do idt=i0-8,i0+8 !Search over +/- one quarter symbol call sync8d(cd0,idt,ctwk,0,sync) if(sync.gt.smax) then smax=sync @@ -109,30 +109,34 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, rxdata=rxdata/rxsig ss=0.84 llr=2.0*rxdata/(ss*ss) - -! do iap=0,3 + apmag=4.0 +! do iap=0,1 do iap=0,0 !### Temporary ### if(iap.eq.0) then apmask=0 apmask(160:162)=1 llrap=llr - llrap(160:162)=5.0*apsym(73:75)/ss + llrap(160:162)=apmag*apsym(73:75)/ss elseif(iap.eq.1) then apmask=0 apmask(88:115)=1 ! mycall + apmask(144)=1 ! not free text apmask(160:162)=1 ! 3 extra bits llrap=0.0 - llrap(88:115)=5.0*apsym(1:28)/ss - llrap(160:162)=5.0*apsym(73:75)/ss + llrap(88:115)=apmag*apsym(1:28)/ss + llrap(144)=-apmag/ss + llrap(160:162)=apmag*apsym(73:75)/ss where(apmask.eq.0) llrap=llr elseif(iap.eq.2) then apmask=0 apmask(88:115)=1 ! mycall apmask(116:143)=1 ! hiscall + apmask(144)=1 ! not free text apmask(160:162)=1 ! 3 extra bits llrap=0.0 - llrap(88:143)=5.0*apsym(1:56)/ss - llrap(160:162)=5.0*apsym(73:75)/ss + llrap(88:143)=apmag*apsym(1:56)/ss + llrap(144)=-apmag/ss + llrap(160:162)=apmag*apsym(73:75)/ss where(apmask.eq.0) llrap=llr elseif(iap.eq.3) then apmask=0 @@ -141,26 +145,28 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, apmask(144:154)=1 ! RRR or 73 apmask(160:162)=1 ! 3 extra bits llrap=0.0 - llrap(88:143)=5.0*apsym(1:56)/ss - llrap(144:154)=5.0*rr73/ss - llrap(160:162)=5.0*apsym(73:75)/ss + llrap(88:143)=apmag*apsym(1:56)/ss + llrap(144:154)=apmag*rr73/ss + llrap(160:162)=apmag*apsym(73:75)/ss where(apmask.eq.0) llrap=llr endif cw=0 call timer('bpd174 ',0) - call bpdecode174(llrap,apmask,max_iterations,decoded,cw,nharderrors) + call bpdecode174(llrap,apmask,max_iterations,decoded,cw,nharderrors,niterations) call timer('bpd174 ',1) dmin=0.0 - if(nharderrors.lt.0 .and. ndepth.ge.2) then + if(ndepth.eq.3 .and. abs(nfqso-f1).lt.10.0 .and. nharderrors.lt.0) then call timer('osd174 ',0) - call osd174(llrap,norder,decoded,cw,nharderrors,dmin) + call osd174(llrap,apmask,norder,decoded,cw,nharderrors,dmin) call timer('osd174 ',1) endif nbadcrc=1 message=' ' xsnr=-99.0 if(count(cw.eq.0).eq.174) cycle !Reject the all-zero codeword - if( nharderrors.ge.0 .and. dmin.le.30.0 .and. nharderrors .lt. 30) then +! if( nharderrors.ge.0 .and. dmin.le.30.0 .and. nharderrors .lt. 30) then +!*** These thresholds should probably be dependent on nap + if( nharderrors.ge.0 .and. dmin.le.50.0 .and. nharderrors .lt. 50) then call chkcrc12a(decoded,nbadcrc) else nharderrors=-1 @@ -169,7 +175,7 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, if(nbadcrc.eq.0) then call extractmessage174(decoded,message,ncrcflag,recent_calls,nrecent) call genft8(message,msgsent,msgbits,itone) -! call subtractft8(dd0,itone,f1,xdt2) + if(lsubtract) call subtractft8(dd0,itone,f1,xdt2) xsig=0.0 xnoi=0.0 do i=1,79 @@ -181,8 +187,8 @@ subroutine ft8b(dd0,newdat,nfqso,ndepth,icand,sync0,f1,xdt,apsym,nharderrors, if( xnoi.gt.0 .and. xnoi.lt.xsig ) xsnr=xsig/xnoi-1.0 xsnr=10.0*log10(xsnr)-27.0 if( xsnr .lt. -24.0 ) xsnr=-24.0 -! write(50,3050) icand,sync0,f1,xdt,nharderrors,dmin,message,iap -!3050 format(i3,3f10.3,i5,f10.3,2x,a22,i3) +! write(50,3050) icand,sync0,f1,xdt,xsnr,nharderrors,niterations,dmin,iap,ipass,message +!3050 format(i3,4f10.3,i5,i5,f10.3,i4,i4,2x,a22) return endif enddo diff --git a/lib/fsk4hf/osd174.f90 b/lib/fsk4hf/osd174.f90 index 1056cf61c..538da351d 100644 --- a/lib/fsk4hf/osd174.f90 +++ b/lib/fsk4hf/osd174.f90 @@ -1,9 +1,10 @@ -subroutine osd174(llr,norder,decoded,cw,nhardmin,dmin) +subroutine osd174(llr,apmask,norder,decoded,cw,nhardmin,dmin) ! ! An ordered-statistics decoder for the (174,87) code. ! include "ldpc_174_87_params.f90" +integer*1 apmask(N),apmaskr(N) integer*1 gen(K,N) integer*1 genmrb(K,N),g2(N,K) integer*1 temp(K),m0(K),me(K),mi(K) @@ -36,6 +37,8 @@ endif ! re-order received vector to place systematic msg bits at the end rx=llr(colorder+1) +apmaskr=apmask(colorder+1) + ! hard decode the received word hdec=0 @@ -91,6 +94,7 @@ hdec=hdec(indices) ! hard decisions from received symbols m0=hdec(1:K) ! zero'th order message absrx=absrx(indices) rx=rx(indices) +apmaskr=apmaskr(indices) s1=sum(absrx(1:K)) s2=sum(absrx(K+1:N)) @@ -110,22 +114,24 @@ do iorder=1,norder mi(K-iorder+1:K)=1 iflag=0 do while(iflag .ge. 0 ) - dpat=sum(mi*absrx(1:K)) - nt=nt+1 - if( dpat .lt. thresh ) then ! reject unlikely error patterns - me=ieor(m0,mi) - call mrbencode(me,ce,g2,N,K) - nxor=ieor(ce,hdec) - dd=sum(nxor*absrx) - if( dd .lt. dmin ) then - dmin=dd - cw=ce - nhardmin=sum(nxor) - thresh=rho*dmin + if(all(iand(apmaskr(1:K),mi).eq.0)) then ! reject patterns with ap bits + dpat=sum(mi*absrx(1:K)) + nt=nt+1 + if( dpat .lt. thresh ) then ! reject unlikely error patterns + me=ieor(m0,mi) + call mrbencode(me,ce,g2,N,K) + nxor=ieor(ce,hdec) + dd=sum(nxor*absrx) + if( dd .lt. dmin ) then + dmin=dd + cw=ce + nhardmin=sum(nxor) + thresh=rho*dmin + endif + else + nrejected=nrejected+1 endif - else - nrejected=nrejected+1 - endif + endif ! get the next test error pattern, iflag will go negative ! when the last pattern with weight iorder has been generated call nextpat(mi,k,iorder,iflag) diff --git a/lib/fsk4hf/subtractft8.f90 b/lib/fsk4hf/subtractft8.f90 index dcc089031..1031cf2b7 100644 --- a/lib/fsk4hf/subtractft8.f90 +++ b/lib/fsk4hf/subtractft8.f90 @@ -10,7 +10,7 @@ subroutine subtractft8(dd,itone,f0,dt) use timer_module, only: timer parameter (NMAX=15*12000,NFRAME=1920*79) - parameter (NFFT=NMAX,NFILT=400) + parameter (NFFT=NMAX,NFILT=1400) real*4 dd(NMAX), window(-NFILT/2:NFILT/2) complex cref,camp,cfilt,cw integer itone(79) diff --git a/lib/fsk4hf/sync8.f90 b/lib/fsk4hf/sync8.f90 index 35c402f58..c625191e2 100644 --- a/lib/fsk4hf/sync8.f90 +++ b/lib/fsk4hf/sync8.f90 @@ -1,7 +1,7 @@ -subroutine sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) +subroutine sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand) include 'ft8_params.f90' - parameter (JZ=31) !DT up to +/- 2.5 s + parameter (JZ=62) !DT up to +/- 1.25 s complex cx(0:NH1) real s(NH1,NHSYM) real savg(NH1) @@ -18,16 +18,13 @@ subroutine sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) data icos7/2,5,6,0,4,1,3/ !Costas 7x7 tone pattern equivalence (x,cx) -! Compute symbol spectra at half-symbol steps. +! Compute symbol spectra, stepping by NSTEP steps. savg=0. - istep=NSPS/2 !960 - tstep=istep/12000.0 !0.08 s + tstep=NSTEP/12000.0 df=12000.0/NFFT1 !3.125 Hz - -! Compute symbol spectra at half-symbol steps fac=1.0/300.0 do j=1,NHSYM - ia=(j-1)*istep + 1 + ia=(j-1)*NSTEP + 1 ib=ia+NSPS-1 x(1:NSPS)=fac*dd(ia:ib) x(NSPS+1:)=0. @@ -45,21 +42,23 @@ subroutine sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) ia=max(1,nint(nfa/df)) ib=nint(nfb/df) + nssy=NSPS/NSTEP ! # steps per symbol + nfos=NFFT1/NSPS ! # frequency bin oversampling factor do i=ia,ib do j=-JZ,JZ t=0. t0=0. do n=0,6 - k=j+2*n - if(k.ge.1) then - t=t + s(i+2*icos7(n),k) - t0=t0 + sum(s(i:i+12:2,k)) + k=j+nssy*n + if(k.ge.1.and.k.le.NHSYM) then + t=t + s(i+nfos*icos7(n),k) + t0=t0 + sum(s(i:i+nfos*6:nfos,k)) endif - t=t + s(i+2*icos7(n),k+72) - t0=t0 + sum(s(i:i+12:2,k+72)) + t=t + s(i+nfos*icos7(n),k+nssy*36) + t0=t0 + sum(s(i:i+nfos*6:nfos,k+nssy*36)) if(k+144.le.NHSYM) then - t=t + s(i+2*icos7(n),k+144) - t0=t0 + sum(s(i:i+12:2,k+144)) + t=t + s(i+nfos*icos7(n),k+nssy*72) + t0=t0 + sum(s(i:i+nfos*6:nfos,k+nssy*72)) endif enddo t0=(t0-t)/6.0 @@ -84,8 +83,7 @@ subroutine sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) candidate0=0. k=0 - syncmin=2.0 - do i=1,100 + do i=1,200 n=ia + indx(iz+1-i) - 1 if(red(n).lt.syncmin) exit if(k.lt.200) k=k+1 @@ -115,13 +113,21 @@ subroutine sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) fac=20.0/maxval(s) s=fac*s +! Sort by sync +! call indexx(candidate0(3,1:ncand),ncand,indx) +! Sort by frequency call indexx(candidate0(1,1:ncand),ncand,indx) + k=1 +! do i=ncand,1,-1 do i=1,ncand j=indx(i) - candidate(1,i)=abs(candidate0(1,j)) - candidate(2,i)=candidate0(2,j) - candidate(3,i)=candidate0(3,j) + if( candidate0(3,j) .ge. syncmin ) then + candidate(1,k)=abs(candidate0(1,j)) + candidate(2,k)=candidate0(2,j) + candidate(3,k)=candidate0(3,j) + k=k+1 + endif enddo - + ncand=k-1 return end subroutine sync8 diff --git a/lib/ft8_decode.f90 b/lib/ft8_decode.f90 index bb751a539..20a7325e2 100644 --- a/lib/ft8_decode.f90 +++ b/lib/ft8_decode.f90 @@ -24,17 +24,18 @@ contains subroutine decode(this,callback,iwave,nfqso,newdat,nutc,nfa, & nfb,nagain,ndepth,nsubmode,mycall12,hiscall12,hisgrid6) -!use wavhdr +! use wavhdr use timer_module, only: timer include 'fsk4hf/ft8_params.f90' -!type(hdr) h +! type(hdr) h class(ft8_decoder), intent(inout) :: this procedure(ft8_decode_callback) :: callback real s(NH1,NHSYM) real candidate(3,200) real dd(15*12000) - logical, intent(in) :: newdat, nagain + logical, intent(in) :: nagain + logical newdat,lsubtract character*12 mycall12, hiscall12 character*6 hisgrid6 integer*2 iwave(15*12000) @@ -50,38 +51,47 @@ contains call ft8apset(mycall12,hiscall12,hisgrid6,apsym) dd=iwave - call timer('sync8 ',0) - call sync8(dd,nfa,nfb,nfqso,s,candidate,ncand) - call timer('sync8 ',1) - syncmin=2.0 - do icand=1,ncand - sync=candidate(3,icand) - if(sync.lt.syncmin) cycle - f1=candidate(1,icand) - xdt=candidate(2,icand) - nsnr0=min(99,nint(10.0*log10(sync) - 25.5)) !### empirical ### - call timer('ft8b ',0) - call ft8b(dd,newdat,nfqso,ndepth,icand,sync,f1,xdt,apsym,nharderrors, & - dmin,nbadcrc,message,xsnr) - nsnr=xsnr - xdt=xdt-0.6 - call timer('ft8b ',1) - if (associated(this%callback)) call this%callback(sync,nsnr,xdt, & +! For now: +! ndepth=1: no subtraction, 1 pass, belief propagation only +! ndepth=2: subtraction, 2 passes, belief propagation only +! ndepth=3: subtraction, 2 passes, bp+osd2 at and near nfqso + if(ndepth.eq.1) npass=1 + if(ndepth.ge.2) npass=2 + do ipass=1,npass + newdat=.true. ! Is this a problem? I hijacked newdat. + if(ipass.eq.1) then + lsubtract=.true. + if(ndepth.eq.1) lsubtract=.false. + syncmin=1.3 + else + lsubtract=.false. + syncmin=1.3 + endif + call timer('sync8 ',0) + call sync8(dd,nfa,nfb,syncmin,nfqso,s,candidate,ncand) + call timer('sync8 ',1) + do icand=1,ncand + sync=candidate(3,icand) + f1=candidate(1,icand) + xdt=candidate(2,icand) + nsnr0=min(99,nint(10.0*log10(sync) - 25.5)) !### empirical ### + call timer('ft8b ',0) + call ft8b(dd,newdat,nfqso,ndepth,lsubtract,icand,sync,f1,xdt,apsym,nharderrors, & + dmin,nbadcrc,iap,ipass,message,xsnr) + nsnr=xsnr + xdt=xdt-0.6 + call timer('ft8b ',1) + if (associated(this%callback)) call this%callback(sync,nsnr,xdt, & f1,nbadcrc,message) -! write(*,'(f7.2,i5,f7.2,f9.1,i5,f7.2,2x,a22)') sync,nsnr,xdt,f1,nharderrors,dmin,message -! write(13,1110) datetime,0,nsnr,xdt,f1,nharderrors,dmin,message -!1110 format(a13,2i4,f6.2,f7.1,i4,' ~ ',f6.2,2x,a22,' FT8') -! write(51,3051) xdt,f1,sync,dmin,nsnr,nharderrors,nbadcrc,message -!3051 format(4f9.1,3i5,2x,a22) -! flush(51) - enddo -!h=default_header(12000,NMAX) -!open(10,file='subtract.wav',status='unknown',access='stream') -!iwave=nint(dd) -!write(10) h,iwave -!close(10) - return + enddo +! h=default_header(12000,NMAX) +! open(10,file='subtract.wav',status='unknown',access='stream') +! iwave=nint(dd) +! write(10) h,iwave +! close(10) + enddo + return end subroutine decode end module ft8_decode diff --git a/lib/jt9.f90 b/lib/jt9.f90 index 394dbdb64..0b4d32849 100644 --- a/lib/jt9.f90 +++ b/lib/jt9.f90 @@ -262,7 +262,7 @@ program jt9 shared_data%params%ndepth=ndepth shared_data%params%dttol=3. - shared_data%params%minsync=0 !### TEST ONLY +! shared_data%params%minsync=0 !### TEST ONLY ! shared_data%params%nfqso=1500 !### TEST ONLY ! mycall="G3WDG " !### TEST ONLY ! hiscall="VK7MO " !### TEST ONLY