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
This commit is contained in:
Steven Franke 2017-07-01 13:08:30 +00:00
parent 2687a308d6
commit 304cbaa672
5 changed files with 175 additions and 33 deletions

View File

@ -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)

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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)