From 304cbaa672b9d0bb72b9facdf4d45691a92eca7b Mon Sep 17 00:00:00 2001 From: Steven Franke Date: Sat, 1 Jul 2017 13:08:30 +0000 Subject: [PATCH] Use complex signal for fine sync and symbol demod. Lower sync threshold. Improved SNR estimate. git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7757 ab8295b8-cf94-4d9e-aec4-7959e3be5d79 --- lib/fsk4hf/ft8_params.f90 | 2 +- lib/fsk4hf/ft8b.f90 | 185 +++++++++++++++++++++++++++++++++----- lib/fsk4hf/ft8sim.f90 | 7 +- lib/fsk4hf/sync8.f90 | 2 +- lib/ft8_decode.f90 | 12 ++- 5 files changed, 175 insertions(+), 33 deletions(-) diff --git a/lib/fsk4hf/ft8_params.f90 b/lib/fsk4hf/ft8_params.f90 index d8dc73969..b139d6913 100644 --- a/lib/fsk4hf/ft8_params.f90 +++ b/lib/fsk4hf/ft8_params.f90 @@ -7,4 +7,4 @@ parameter (NSPS=2048) !Samples per symbol at 12000 S/s parameter (NZ=NSPS*NN) !Samples in full 15 s waveform (161,792) 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/NSPS-1) !Number of symbol spectra (half-symbol steps) +parameter (NHSYM=2*NMAX/NH1-1) !Number of symbol spectra (1/2-symbol steps) diff --git a/lib/fsk4hf/ft8b.f90 b/lib/fsk4hf/ft8b.f90 index b232547c0..ef3e3d630 100644 --- a/lib/fsk4hf/ft8b.f90 +++ b/lib/fsk4hf/ft8b.f90 @@ -1,37 +1,134 @@ -subroutine ft8b(s,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message) +subroutine ft8b(dd0,s,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message,xsnr) use timer_module, only: timer include 'ft8_params.f90' - parameter(NRECENT=10) + parameter(NRECENT=10,NP2=2812) character*12 recent_calls(NRECENT) - character message*22 + character message*22,msgsent*22 + real a(5) real s(NH1,NHSYM) - real s1(0:7,ND) + real s1(0:7,ND),s2(0:7,NN) real ps(0:7) real rxdata(3*ND),llr(3*ND) !Soft symbols + real dd0(15*12000) integer*1 decoded(KK),apmask(3*ND),cw(3*ND) + integer itone(NN) + integer icos7(0:6) + complex cd0(NP2) + complex csync(0:6,32) + complex ctwk(32) + complex csymb(32) + logical first + data icos7/2,5,6,0,4,1,3/ + data first/.true./ + save first,twopi,fs2,dt2,taus,baud,csync + + if( first ) then + twopi=8.0*atan(1.0) + fs2=12000.0/64.0 + dt2=1/fs2 + taus=32*dt2 + baud=1/taus + do i=0,6 + phi=0.0 + dphi=twopi*icos7(i)*baud*dt2 + do j=1,32 + csync(i,j)=cmplx(cos(phi),sin(phi)) + phi=mod(phi+dphi,twopi) + enddo + enddo + first=.false. + endif max_iterations=40 norder=2 ! if(abs(nfqso-f1).lt.10.0) norder=3 tstep=0.5*NSPS/12000.0 df=12000.0/NFFT1 + call ft8_downsample(dd0,f1,cd0) - i0=max(1,nint(f1/df)) - j0=nint(xdt/tstep) + i0=xdt*fs2 + smax=0.0 + do idt=i0-16,i0+16 + sync=0 + do i=0,6 + i1=idt+i*32 + i2=i1+36*32 + i3=i1+72*32 + term1=0.0 ! this needs to be fixed... + term2=0.0 + term3=0.0 + if( i1.ge.1 .and. i1+31.le.NP2 ) & + term1=abs(sum(cd0(i1:i1+31)*conjg(csync(i,1:32)))) + if( i2.ge.1 .and. i2+31.le.NP2 ) & + term2=abs(sum(cd0(i2:i2+31)*conjg(csync(i,1:32)))) + if( i3.ge.1 .and. i3+31.le.NP2 ) & + term3=abs(sum(cd0(i3:i3+31)*conjg(csync(i,1:32)))) + sync=sync+term1+term2+term3 + enddo + if( sync .gt. smax ) then + smax=sync + ibest=idt + delfbest=delf + ifbest=if + endif + enddo + xdt2=ibest*dt2 + +! peak up the frequency + i0=xdt2*fs2 + smax=0.0 + do ifr=-5,5 + delf=ifr*0.5 + dphi=twopi*delf*dt2 + phi=0.0 + do i=1,32 + ctwk(i)=cmplx(cos(phi),sin(phi)) + phi=mod(phi+dphi,twopi) + enddo + sync=0.0 + do i=0,6 + i1=i0+i*32 + i2=i1+36*32 + i3=i1+72*32 + term1=0.0 ! this needs to be fixed... + term2=0.0 + term3=0.0 + if( i1.ge.1 .and. i1+31.le.NP2 ) & + term1=abs(sum(cd0(i1:i1+31)*conjg(ctwk*csync(i,1:32)))) + if( i2.ge.1 .and. i2+31.le.NP2 ) & + term2=abs(sum(cd0(i2:i2+31)*conjg(ctwk*csync(i,1:32)))) + if( i3.ge.1 .and. i3+31.le.NP2 ) & + term3=abs(sum(cd0(i3:i3+31)*conjg(ctwk*csync(i,1:32)))) + sync=sync+term1+term2+term3 + enddo + if( sync .gt. smax ) then + smax=sync + delfbest=delf + endif + enddo + a=0.0 + a(1)=-delfbest + call twkfreq1(cd0,NP2,fs2,a,cd0) + xdt=xdt2 + f1=f1+delfbest j=0 - ia=i0 - ib=i0+14 do k=1,NN - if(k.le.7) cycle - if(k.ge.37 .and. k.le.43) cycle - if(k.gt.72) cycle - n=j0+2*(k-1)+1 - if(n.lt.1) cycle - j=j+1 - s1(0:7,j)=s(ia:ib:2,n) - enddo + i1=ibest+(k-1)*32 + csymb=cd0(i1:i1+31) + call four2a(csymb,32,1,-1,1) + s2(0:7,k)=abs(csymb(1:8)) + enddo + j=0 + do k=1,NN + if(k.le.7) cycle + if(k.ge.37 .and. k.le.43) cycle + if(k.gt.72) cycle + j=j+1 + s1(0:7,j)=s2(0:7,k) + enddo + do j=1,ND ps=s1(0:7,j) where (ps.gt.0.0) ps=log(ps) @@ -56,9 +153,6 @@ subroutine ft8b(s,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message) llr=2.0*rxdata/(ss*ss) apmask=0 cw=0 -! cw will be needed for subtraction. -! dmin is the correlation discrepancy of a returned codeword - it is -! used to select the best codeword within osd174. call timer('bpd174 ',0) call bpdecode174(llr,apmask,max_iterations,decoded,cw,nharderrors) call timer('bpd174 ',1) @@ -67,17 +161,62 @@ subroutine ft8b(s,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message) call timer('osd174 ',0) call osd174(llr,norder,decoded,cw,nharderrors,dmin) call timer('osd174 ',1) -! This threshold needs to be tuned. 99.0 should pass everything. - if( dmin .gt. 99.0 ) nharderrors=-1 endif nbadcrc=1 message=' ' + xsnr=-99.0 if(count(cw.eq.0).eq.174) go to 900 !Reject the all-zero codeword - if(nharderrors.ge.0) call chkcrc12a(decoded,nbadcrc) + if( nharderrors.ge.0 .and. dmin.le.30.0 .and. nharderrors .lt. 30) then + call chkcrc12a(decoded,nbadcrc) + else + nharderrors=-1 + go to 900 + endif if(nbadcrc.eq.0) then call extractmessage174(decoded,message,ncrcflag,recent_calls,nrecent) + call genft8(message,msgsent,itone) + xsig=0.0 + xnoi=0.0 + do i=1,79 + xsig=xsig+s2(itone(i),i)**2 + ios=mod(itone(i)+4,7) + xnoi=xnoi+s2(ios,i)**2 + enddo + xsnr=xsig/xnoi-1.0 + if( xsnr .lt. 0.0 ) xsnr=0.005 + xsnr=10.0*log10(xsnr)-26.3 endif 900 continue - return end subroutine ft8b + +subroutine ft8_downsample(dd,f0,c1) + +! Downconvert to complex data sampled at 187.5 Hz, 32 samples/symbol + + parameter (NMAX=15*12000,NFFT2=2812) + complex c1(0:NFFT2-1) + complex cx(0:NMAX/2-1) + real dd(NMAX),x(NMAX) + equivalence (x,cx) + + df=12000.0/NMAX + x=dd + call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain + baud=12000.0/(32.0*64.0) + i0=nint(f0/df) + ft=f0+8.0*baud + it=nint(ft/df) + fb=f0-1.0*baud + ib=nint(fb/df) + k=0 + c1=cmplx(0.0,0.0) + do i=ib,it + c1(k)=cx(i) + k=k+1 + enddo + c1=cshift(c1,i0-ib) + call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain + c1=c1/(NMAX*NFFT2) + return +end subroutine ft8_downsample diff --git a/lib/fsk4hf/ft8sim.f90 b/lib/fsk4hf/ft8sim.f90 index 393b22825..29769724f 100644 --- a/lib/fsk4hf/ft8sim.f90 +++ b/lib/fsk4hf/ft8sim.f90 @@ -58,7 +58,7 @@ program ft8sim f0=(isig+2)*100.0 phi=0.0 k=-1 + nint(xdt/dt) - do j=1,NN !Generate complex waveform + do j=1,NN !Generate complex waveform dphi=twopi*(f0+itone(j)*baud)*dt if(k.eq.0) phi=-dphi do i=1,NSPS @@ -69,12 +69,12 @@ program ft8sim if(k.ge.0 .and. k.lt.NMAX) c0(k)=cmplx(cos(xphi),sin(xphi)) enddo enddo - if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NZ,fs,delay,fspread) + if(fspread.ne.0.0 .or. delay.ne.0.0) call watterson(c,NMAX,fs,delay,fspread) c=c+c0 enddo c=c*sig if(snrdb.lt.90) then - do i=0,NZ-1 !Add gaussian noise at specified SNR + do i=0,NMAX-1 !Add gaussian noise at specified SNR xnoise=gran() ynoise=gran() c(i)=c(i) + cmplx(xnoise,ynoise) @@ -85,7 +85,6 @@ program ft8sim rms=100.0 if(snrdb.ge.90.0) iwave(1:NMAX)=nint(fac*real(c)) if(snrdb.lt.90.0) iwave(1:NMAX)=nint(rms*real(c)) - iwave(NZ+1:)=0 h=default_header(12000,NMAX) write(fname,1102) ifile diff --git a/lib/fsk4hf/sync8.f90 b/lib/fsk4hf/sync8.f90 index 671e83648..eb9a1f2a2 100644 --- a/lib/fsk4hf/sync8.f90 +++ b/lib/fsk4hf/sync8.f90 @@ -80,7 +80,7 @@ subroutine sync8(iwave,nfa,nfb,nfqso,s,candidate,ncand) candidate0=0. k=0 - syncmin=4.0 + syncmin=1.5 do i=1,100 n=ia + indx(iz+1-i) - 1 if(red(n).lt.syncmin) exit diff --git a/lib/ft8_decode.f90 b/lib/ft8_decode.f90 index 19c098ad4..4ba59f2e0 100644 --- a/lib/ft8_decode.f90 +++ b/lib/ft8_decode.f90 @@ -33,10 +33,11 @@ contains real ss(1,1) !### dummy, to be removed ### real s(NH1,NHSYM) real candidate(3,100) + real dd(15*12000) logical, intent(in) :: newdat, nagain integer*2 iwave(15*12000) character datetime*13,message*22 - + this%callback => callback write(datetime,1001) nutc !### TEMPORARY ### 1001 format("000000_",i6.6) @@ -45,7 +46,8 @@ contains call sync8(iwave,nfa,nfb,nfqso,s,candidate,ncand) call timer('sync8 ',1) - syncmin=4.0 + dd=iwave + syncmin=2.0 ! rewind 51 do icand=1,ncand sync=candidate(3,icand) @@ -54,13 +56,15 @@ contains xdt=candidate(2,icand) nsnr=min(99,nint(10.0*log10(sync) - 25.5)) !### empirical ### call timer('ft8b ',0) - call ft8b(s,nfqso,f1,xdt,nharderrors,dmin,nbadcrc,message) + call ft8b(dd,s,nfqso,f1,xdt,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, & 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') +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)