More work on real time decoder.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7105 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2016-09-22 00:03:44 +00:00
parent 5df0bbcd50
commit be0c4f15c1
6 changed files with 459 additions and 77 deletions

View File

@ -67,7 +67,7 @@ subroutine hspec(id2,k,nutc0,ntrperiod,ntol,bmsk144,ingain,green,s,jh,line1)
!###
if(bmsk144) then
if(k.ge.7168) then
tsec=(k-3584)/12000.0
tsec=(k-7168)/12000.0
call mskrtd(id2(k-7168+1:k),nutc0,tsec,ntol,line1)
endif
endif

View File

@ -84,9 +84,9 @@ program msk144d2
call timer('read ',1)
do i=1,npts,7*512
ichunk=id2(i:i+7*1024-1)
tsec=(i-1+7*512)/12000.0
tsec=(i-1)/12000.0
call mskrtd(ichunk,nutc,tsec,ntol,line)
if( line .ne. ' ' ) then
if( index(line,"^") .ne. 0 .or. index(line,"&") .ne. 0 ) then
write(*,*) line
endif
enddo

View File

@ -1,5 +1,5 @@
subroutine msk144decodeframe(c,msgreceived,nsuccess)
use timer_module, only: timer
! use timer_module, only: timer
parameter (NSPM=864)
character*22 msgreceived
@ -16,10 +16,9 @@ subroutine msk144decodeframe(c,msgreceived,nsuccess)
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
save df,first,cb,fs,pi,twopi,dt,s8,pp,nmatchedfilter
save df,first,cb,fs,pi,twopi,dt,s8,pp
if(first) then
nmatchedfilter=1
! define half-sine pulse and raised-cosine edge window
pi=4d0*datan(1d0)
twopi=8d0*datan(1d0)
@ -95,9 +94,9 @@ subroutine msk144decodeframe(c,msgreceived,nsuccess)
max_iterations=10
max_dither=1
call timer('bpdec144 ',0)
! call timer('bpdec144 ',0)
call bpdecode144(llr,max_iterations,decoded,niterations)
call timer('bpdec144 ',1)
! call timer('bpdec144 ',1)
if( niterations .ge. 0.0 ) then
call extractmessage144(decoded,msgreceived,nhashflag)

312
lib/msk144spd.f90 Normal file
View File

@ -0,0 +1,312 @@
subroutine msk144spd(cbig,n,ntol,nsuccess,msgreceived,fret,snrret,tret)
! msk144 short-ping-decoder
use timer_module, only: timer
parameter (NSPM=864, NPTS=3*NSPM, MAXSTEPS=1700, NFFT=NSPM, MAXCAND=5)
character*22 msgreceived
complex cbig(n)
complex cdat(NPTS) !Analytic signal
complex cdat2(NPTS)
complex c(NSPM)
complex ctmp(NFFT)
complex cb(42) !Complex waveform for sync word
complex cbr(42) !Complex waveform for reversed sync word
complex cfac,cca,ccb
complex cc(NPTS)
complex ccr(NPTS)
complex cc1(NPTS)
complex cc2(NPTS)
complex ccr1(NPTS)
complex ccr2(NPTS)
complex bb(6)
integer s8(8),s8r(8)
integer, dimension(1) :: iloc
integer indices(MAXSTEPS)
integer ipeaks(10)
logical ismask(NFFT)
real cbi(42),cbq(42)
real detmet(-2:MAXSTEPS+3)
real detmet2(-2:MAXSTEPS+3)
real detfer(MAXSTEPS)
real rcw(12)
real dd(NPTS)
real ferrs(MAXCAND)
real pp(12) !Half-sine pulse shape
real snrs(MAXCAND)
real times(MAXCAND)
real tonespec(NFFT)
real*8 dt, df, fs, pi, twopi
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
data s8r/1,0,1,1,0,0,0,1/
save df,first,cb,fs,pi,twopi,dt,s8,rcw,pp,nmatchedfilter
if(first) then
nmatchedfilter=1
! define half-sine pulse and raised-cosine edge window
pi=4d0*datan(1d0)
twopi=8d0*datan(1d0)
fs=12000.0
dt=1.0/fs
df=fs/NFFT
do i=1,12
angle=(i-1)*pi/12.0
pp(i)=sin(angle)
rcw(i)=(1-cos(angle))/2
enddo
! define the sync word waveforms
s8=2*s8-1
cbq(1:6)=pp(7:12)*s8(1)
cbq(7:18)=pp*s8(3)
cbq(19:30)=pp*s8(5)
cbq(31:42)=pp*s8(7)
cbi(1:12)=pp*s8(2)
cbi(13:24)=pp*s8(4)
cbi(25:36)=pp*s8(6)
cbi(37:42)=pp(1:6)*s8(8)
cb=cmplx(cbi,cbq)
s8r=2*s8r-1
cbq(1:6)=pp(7:12)*s8r(1)
cbq(7:18)=pp*s8r(3)
cbq(19:30)=pp*s8r(5)
cbq(31:42)=pp*s8r(7)
cbi(1:12)=pp*s8r(2)
cbi(13:24)=pp*s8r(4)
cbi(25:36)=pp*s8r(6)
cbi(37:42)=pp(1:6)*s8r(8)
cbr=cmplx(cbi,cbq)
first=.false.
endif
! fill the detmet, detferr arrays
nstep=(n-NPTS)/216 ! 72ms/4=18ms steps
detmet=0
detmet2=0
detfer=-999.99
do istp=1,nstep
ns=1+216*(istp-1)
ne=ns+NSPM-1
if( ne .gt. n ) exit
ctmp=cmplx(0.0,0.0)
ctmp(1:NSPM)=cbig(ns:ne)
! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in
! squared signal spectrum.
! search range for coarse frequency error is +/- 100 Hz
ctmp=ctmp**2
ctmp(1:12)=ctmp(1:12)*rcw
ctmp(NSPM-11:NSPM)=ctmp(NSPM-11:NSPM)*rcw(12:1:-1)
call four2a(ctmp,NFFT,1,-1,1)
tonespec=abs(ctmp)**2
ihlo=(4000-2*ntol)/df+1
ihhi=(4000+2*ntol)/df+1
ismask=.false.
ismask(ihlo:ihhi)=.true. ! high tone search window
iloc=maxloc(tonespec,ismask)
ihpk=iloc(1)
deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
ah=tonespec(ihpk)
ahavp=(sum(tonespec,ismask)-ah)/count(ismask)
trath=ah/(ahavp+0.01)
illo=(2000-2*ntol)/df+1
ilhi=(2000+2*ntol)/df+1
ismask=.false.
ismask(illo:ilhi)=.true. ! window for low tone
iloc=maxloc(tonespec,ismask)
ilpk=iloc(1)
deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
al=tonespec(ilpk)
alavp=(sum(tonespec,ismask)-al)/count(ismask)
tratl=al/(alavp+0.01)
fdiff=(ihpk+deltah-ilpk-deltal)*df
i2000=2000/df+1
i4000=4000/df+1
ferrh=(ihpk+deltah-i4000)*df/2.0
ferrl=(ilpk+deltal-i2000)*df/2.0
if( ah .ge. al ) then
ferr=ferrh
else
ferr=ferrl
endif
detmet(istp)=max(ah,al)
detmet2(istp)=max(trath,tratl)
detfer(istp)=ferr
enddo ! end of detection-metric and frequency error estimation loop
call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
xmed=detmet(indices(nstep/4))
detmet=detmet/xmed ! noise floor of detection metric is 1.0
ndet=0
do ip=1,MAXCAND ! Find candidates
iloc=maxloc(detmet(1:nstep))
il=iloc(1)
! if( (detmet(il) .lt. 4.0) ) exit
if( (detmet(il) .lt. 3.0) ) exit
if( abs(detfer(il)) .le. ntol ) then
ndet=ndet+1
times(ndet)=((il-1)*216+NSPM/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet(il))/2-9.0
endif
! detmet(max(1,il-1):min(nstep,il+1))=0.0
detmet(il)=0.0
enddo
if( ndet .lt. 3 ) then
do ip=1,MAXCAND-ndet ! Find candidates
iloc=maxloc(detmet2(1:nstep))
il=iloc(1)
if( (detmet2(il) .lt. 12.0) ) exit
if( abs(detfer(il)) .le. ntol ) then
ndet=ndet+1
times(ndet)=((il-1)*216+NSPM/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet2(il))/2-9.0
endif
! detmet2(max(1,il-1):min(nstep,il+1))=0.0
detmet2(il)=0.0
enddo
endif
nsuccess=0
msgreceived=' '
do ip=1,ndet ! Try to sync/demod/decode each candidate.
imid=times(ip)*fs
if( imid .lt. NPTS/2 ) imid=NPTS/2
if( imid .gt. n-NPTS/2 ) imid=n-NPTS/2
cdat=cbig(imid-NPTS/2+1:imid+NPTS/2)
ferr=ferrs(ip)
! remove coarse freq error - should now be within a few Hz
call tweak1(cdat,NPTS,-(1500+ferr),cdat)
! attempt frame synchronization
! correlate with sync word waveforms
cc=0
ccr=0
cc1=0
cc2=0
ccr1=0
ccr2=0
do i=1,NPTS-(56*6+41)
cc1(i)=sum(cdat(i:i+41)*conjg(cb))
cc2(i)=sum(cdat(i+56*6:i+56*6+41)*conjg(cb))
enddo
cc=cc1+cc2
dd=abs(cc1)*abs(cc2)
cmax=maxval(abs(cc))
! Find 6 largest peaks
do ipk=1, 6
iloc=maxloc(abs(cc))
ic1=iloc(1)
iloc=maxloc(dd)
ic2=iloc(1)
ipeaks(ipk)=ic2
dd(max(1,ic2-7):min(NPTS-56*6-41,ic2+7))=0.0
enddo
do ipk=1,4
! we want ic to be the index of the first sample of the frame
ic0=ipeaks(ipk)
! fine adjustment of sync index
do i=1,6
if( ic0+11+NSPM .le. NPTS ) then
bb(i) = sum( ( cdat(ic0+i-1+6:ic0+i-1+6+NSPM:6) * conjg( cdat(ic0+i-1:ic0+i-1+NSPM:6) ) )**2 )
else
bb(i) = sum( ( cdat(ic0+i-1+6:NPTS:6) * conjg( cdat(ic0+i-1:NPTS-6:6) ) )**2 )
endif
enddo
iloc=maxloc(abs(bb))
ibb=iloc(1)
bba=abs(bb(ibb))
bbp=atan2(-imag(bb(ibb)),-real(bb(ibb)))/(2*twopi*6*dt)
if( ibb .le. 3 ) ibb=ibb-1
if( ibb .gt. 3 ) ibb=ibb-7
do id=1,3 ! Slicer dither.
if( id .eq. 1 ) is=0
if( id .eq. 2 ) is=-1
if( id .eq. 3 ) is=1
! Adjust frame index to place peak of bb at desired lag
ic=ic0+ibb+is
if( ic .lt. 1 ) ic=ic+864
! Estimate fine frequency error.
! Should a larger separation be used when frames are averaged?
cca=sum(cdat(ic:ic+41)*conjg(cb))
if( ic+56*6+41 .le. NPTS ) then
ccb=sum(cdat(ic+56*6:ic+56*6+41)*conjg(cb))
cfac=ccb*conjg(cca)
ferr2=atan2(imag(cfac),real(cfac))/(twopi*56*6*dt)
else
ccb=sum(cdat(ic-88*6:ic-88*6+41)*conjg(cb))
cfac=cca*conjg(ccb)
ferr2=atan2(imag(cfac),real(cfac))/(twopi*88*6*dt)
endif
! Final estimate of the carrier frequency - returned to the calling program
fest=1500+ferr+ferr2
do idf=0,4 ! frequency jitter
if( idf .eq. 0 ) then
deltaf=0.0
elseif( mod(idf,2) .eq. 0 ) then
deltaf=idf
else
deltaf=-(idf+1)
endif
! Remove fine frequency error
call tweak1(cdat,NPTS,-(ferr2+deltaf),cdat2)
! place the beginning of frame at index NSPM+1
cdat2=cshift(cdat2,ic-(NSPM+1))
do iav=1,7 ! Hopefully we can eliminate some of these after looking at more examples
if( iav .eq. 1 ) then
c=cdat2(NSPM+1:2*NSPM)
elseif( iav .eq. 2 ) then
c=cdat2(NSPM-431:NSPM+432)
c=cshift(c,-432)
elseif( iav .eq. 3 ) then
c=cdat2(2*NSPM-431:2*NSPM+432)
c=cshift(c,-432)
elseif( iav .eq. 4 ) then
c=cdat2(1:NSPM)
elseif( iav .eq. 5 ) then
c=cdat2(2*NSPM+1:NPTS)
elseif( iav .eq. 6 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)
elseif( iav .eq. 7 ) then
c=cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
elseif( iav .eq. 8 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
endif
call msk144decodeframe(c,msgreceived,nsuccess)
if( nsuccess .eq. 1 ) then
fret=1500+ferrs(ip)
snrret=snrs(ip)
tret=times(ip)
return
endif
enddo ! frame averaging loop
enddo ! frequency dithering loop
enddo ! sample-time dither loop
enddo ! peak loop
enddo ! candidate loop
return
end subroutine msk144spd

102
lib/msk144sync.f90 Normal file
View File

@ -0,0 +1,102 @@
subroutine msk144sync(cdat,n,ntol,ndf,iavmask,npeaks,fest,snr,ipklocs,pkamps,c)
parameter (NSPM=864)
complex cdat(n)
complex cdat2(n)
complex c(NSPM) !Coherently averaged complex data
complex ct2(2*NSPM)
complex cs(NSPM)
complex cb(42) !Complex waveform for sync word
complex cc(0:NSPM-1)
integer s8(8)
integer iloc(1)
integer ipklocs(npeaks)
integer iavmask(8) ! defines which frames to average
real cbi(42),cbq(42)
real pkamps(npeaks)
real xcc(0:NSPM-1)
real xccs(0:NSPM-1)
real pp(12) !Half-sine pulse shape
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
save first,cb,fs,pi,twopi,dt,s8,pp
! call system_clock(count0,clkfreq)
if(first) then
pi=4.0*atan(1.0)
twopi=8.0*atan(1.0)
fs=12000.0
dt=1.0/fs
do i=1,12 !Define half-sine pulse
angle=(i-1)*pi/12.0
pp(i)=sin(angle)
enddo
! Define the sync word waveforms
s8=2*s8-1
cbq(1:6)=pp(7:12)*s8(1)
cbq(7:18)=pp*s8(3)
cbq(19:30)=pp*s8(5)
cbq(31:42)=pp*s8(7)
cbi(1:12)=pp*s8(2)
cbi(13:24)=pp*s8(4)
cbi(25:36)=pp*s8(6)
cbi(37:42)=pp(1:6)*s8(8)
cb=cmplx(cbi,cbq)
first=.false.
endif
navg=sum(iavmask)
xmax=0.0
bestf=0.0
do ifr=-ntol,ntol,ndf !Find freq that maximizes sync
ferr=ifr
call tweak1(cdat,n,-(1500+ferr),cdat2)
c=0
do i=1,8
ib=(i-1)*NSPM+1
ie=ib+NSPM-1
if( iavmask(i) .eq. 1 ) then
c(1:NSPM)=c(1:NSPM)+cdat2(ib:ie)
endif
enddo
cc=0
ct2(1:NSPM)=c
ct2(NSPM+1:2*NSPM)=c
do ish=0,NSPM-1
cc(ish)=dot_product(ct2(1+ish:42+ish)+ct2(336+ish:377+ish),cb(1:42))
enddo
xcc=abs(cc)
xb=maxval(xcc)/(48.0*sqrt(float(navg)))
if(xb.gt.xmax) then
xmax=xb
bestf=ferr
cs=c
xccs=xcc
endif
enddo
fest=1500+bestf
c=cs
xcc=xccs
! Find npeaks largest peaks
do ipk=1,npeaks
iloc=maxloc(xcc)
ic2=iloc(1)
ipklocs(ipk)=ic2
pkamps(ipk)=xcc(ic2-1)
xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
enddo
!write(*,*) xmax,bestf,fest,pkamps(1)/(48.0*sqrt(float(navg))),ipklocs(1),ipklocs(2)
snr=-6.0
return
end subroutine msk144sync

View File

@ -10,34 +10,41 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,line)
parameter (NAVGMAX=7) !Coherently average up to 7 frames
parameter (NPTSMAX=7*NSPM) !Max points analyzed at once
integer*2 id2(NZ) !Raw 16-bit data
character*3 decsym !"&" for mskspd or "^" for long averages
character*22 msgreceived !Decoded message
character*80 line !Formatted line with UTC dB T Freq Msg
complex cdat(NFFT1) !Analytic signal
complex cdat2(NFFT1) !Signal shifted to baseband
complex c(NSPM) !Coherently averaged complex data
complex ct(NSPM)
complex ct2(2*NSPM)
complex cs(NSPM)
complex cb(42) !Complex waveform for sync word
complex cc(0:NSPM-1)
! integer*8 count0,count1,count2,count3,clkfreq
integer*2 id2(NZ) !Raw 16-bit data
integer iavmask(8)
integer iavpatterns(8,6)
integer s8(8)
integer iloc(1)
integer ipeaks(10)
integer nav(6)
real cbi(42),cbq(42)
real d(NFFT1)
real xcc(0:NSPM-1)
real xccs(0:NSPM-1)
real pkamps(10)
real pp(12) !Half-sine pulse shape
real xmc(6)
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
data nav/1,2,3,5,7,9/
data iavpatterns/ &
1,1,1,0,0,0,0,0, &
0,1,1,1,0,0,0,0, &
0,0,1,1,1,0,0,0, &
1,1,1,1,1,0,0,0, &
0,0,1,1,1,1,1,0, &
1,1,1,1,1,1,1,0/
data xmc/1.5,2.5,3.5,2.5,4.5,3.5/ !Used to label decode with time at center of averaging mask
save first,cb,fs,pi,twopi,dt,s8,pp,t03,t12,nutc00
! call system_clock(count0,clkfreq)
@ -80,58 +87,28 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,line)
d(1:NZ)=fac*d(1:NZ)
d(NZ+1:NFFT1)=0.
call analytic(d,NZ,NFFT1,cdat) !Convert to analytic signal and filter
np=7*NSPM
call msk144spd(cdat,np,ntol,nsuccess,msgreceived,fest,snr,tdec)
if( nsuccess .eq. 1 ) then
tdec=tsec+tdec
decsym=' & '
goto 999
endif
tframe=float(NSPM)/12000.0
nmessages=0
line=char(0)
nshort=0
npts=7168
nsnr=-4 !### Temporary ###
do iavg=1,5
navg=nav(iavg)
ndf=nint(7.0/navg) + 1
xmax=0.0
bestf=0.0
! call system_clock(count1,clkfreq)
do ifr=-ntol,ntol,ndf !Find freq that maximizes sync
ferr=ifr
call tweak1(cdat,NPTS,-(1500+ferr),cdat2)
c=0
do i=1,navg
ib=(i-1)*NSPM+1
ie=ib+NSPM-1
c(1:NSPM)=c(1:NSPM)+cdat2(ib:ie)
enddo
do iavg=1,6
iavmask=iavpatterns(1:8,iavg)
navg=sum(iavmask)
! ndf=nint(7.0/navg) + 1
ndf=nint(7.0/navg)
cc=0
ct2(1:NSPM)=c
ct2(NSPM+1:2*NSPM)=c
do ish=0,NSPM-1
cc(ish)=dot_product(ct2(1+ish:42+ish)+ct2(336+ish:377+ish),cb(1:42))
enddo
xcc=abs(cc)
xb=maxval(xcc)/(48.0*sqrt(float(navg)))
if(xb.gt.xmax) then
xmax=xb
bestf=ferr
cs=c
xccs=xcc
endif
enddo
! call system_clock(count2,clkfreq)
fest=1500+bestf
c=cs
xcc=xccs
! Find 2 largest peaks
do ipk=1,2
iloc=maxloc(xcc)
ic2=iloc(1)
ipeaks(ipk)=ic2
xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0
enddo
npeaks=2
call msk144sync(cdat(1:8*NSPM),8*864,ntol,ndf,iavmask,npeaks,fest,snr,ipeaks,pkamps,c)
do ipk=1,2
do is=1,3
@ -143,8 +120,8 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,line)
call msk144decodeframe(ct,msgreceived,nsuccess)
if(nsuccess .gt. 0) then
write(line,1020) nutc0,nsnr,tsec,nint(fest),msgreceived,char(0)
1020 format(i6.6,i4,f5.1,i5,' ^ ',a22,a1)
tdec=tsec+xmc(iavg)*tframe
decsym=' ^ '
goto 999
endif
enddo !Slicer dither
@ -152,18 +129,10 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,line)
enddo
msgreceived=' '
ndither=-98
return
999 continue
! call system_clock(count3,clkfreq)
! t12=t12 + float(count2-count1)/clkfreq
! t03=t03 + float(count3-count0)/clkfreq
! if(navg.gt.7) navg=0
! write(*,3002) nutc0,tsec,t12,t03,xmax,nint(bestf),navg, &
! nbadsync,niterations,ipk,is,msgreceived(1:19)
! write(62,3002) nutc0,tsec,t12,t03,xmax,nint(bestf),navg, &
! nbadsync,niterations,ipk,is,msgreceived(1:19)
!3002 format(i6,f6.2,2f7.2,f6.2,i5,5i3,1x,a19)
nsnr=nint(snr)
write(line,1020) nutc0,nsnr,tdec,nint(fest),decsym,msgreceived,char(0)
1020 format(i6.6,i4,f5.1,i5,a3,a22,a1)
return
end subroutine mskrtd