new combined detector/demodulator/decoder routine for msk144

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6768 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Steven Franke 2016-06-12 02:42:51 +00:00
parent f585500505
commit d45dc65a87

373
lib/detectmsk144.f90 Normal file
View File

@ -0,0 +1,373 @@
subroutine detectmsk144(cbig,n,pchk_file,lines,nmessages,nutc)
!nutc and t0 are for debug output
use iso_c_binding, only: c_loc,c_size_t
use packjt
use hashing
use timer_module, only: timer
parameter (NSPM=864, NPTS=3*NSPM)
character*22 msgreceived,allmessages(20)
character*80 lines(100)
character*512 pchk_file,gen_file
complex cbig(n)
complex cdat(NPTS) !Analytic signal
complex cdat2(NPTS)
complex c(NSPM)
complex ctmp(6000)
complex cb(42) !Complex waveform for sync word
complex cfac,cca,ccb
complex cc(NPTS)
complex cc1(NPTS)
complex cc2(NPTS)
complex bb(6)
integer s8(8),hardbits(144)
integer*1, target:: i1Dec8BitBytes(10)
integer, dimension(1) :: iloc
integer*4 i4Dec6BitWords(12)
integer*1 decoded(80)
integer*1 i1hashdec
integer indices(700)
integer ipeaks(10)
logical ismask(6000)
real cbi(42),cbq(42)
real detmet(700)
real detfer(700)
real tonespec(6000)
real rcw(12)
real dd(NPTS)
real ferrs(20)
real pp(12) !Half-sine pulse shape
real snrs(20)
real times(20)
real*8 dt, df, fs, pi, twopi
real softbits(144)
real*8 unscrambledsoftbits(128)
real lratio(128)
real ttt(NPTS)
logical first
data first/.true./
data s8/0,1,1,1,0,0,1,0/
save df,first,cb,fs,nfft,pi,twopi,dt,s8,rcw,pp
if(first) then
allmessages=' '
lines=' '
i=index(pchk_file,".pchk")
gen_file=pchk_file(1:i-1)//".gen"
call init_ldpc(trim(pchk_file)//char(0),trim(gen_file)//char(0))
! define half-sine pulse and raised-cosine edge window
pi=4d0*datan(1d0)
twopi=8d0*datan(1d0)
fs=12000.0
dt=1.0/fs
nfft=6000 !using a zero-padded fft to get 2 Hz bins
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 waveform
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)
ttt=(/ (i, i=1,NPTS) /)
ttt=(ttt-NPTS/2)/(NPTS/2)
first=.false.
endif
nstep=694 ! fill the detmet, detsnr, detferr arrays
do istp=1,nstep
ns=1+256*(istp-1)
ne=ns+NPTS-1
tt=(ns+ne)/2.0/12000.0
cdat=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=cmplx(0.0,0.0)
ctmp(1:NPTS)=cdat**2
ctmp(1:12)=ctmp(1:12)*rcw
ctmp(NPTS-11:NPTS)=ctmp(NPTS-11:NPTS)*rcw(12:1:-1)
call four2a(ctmp,nfft,1,-1,1)
tonespec=abs(ctmp)**2
ismask=.false.
ismask(1901:2101)=.true. ! high tone search window
iloc=maxloc(tonespec,ismask)
ihpk=iloc(1)
ah=tonespec(ihpk)
ismask=.false.
ismask(901:1101)=.true. ! window for low tone
iloc=maxloc(tonespec,ismask)
ilpk=iloc(1)
al=tonespec(ilpk)
fdiff=(ihpk-ilpk)*df
ferrh=(ihpk-2001)*df/2.0
ferrl=(ilpk-1001)*df/2.0
if( abs(fdiff-2000) .le. 16.0 ) then
if( ah .ge. al ) then
ferr=ferrh
else
ferr=ferrl
endif
else
ferr=-999.99
endif
detmet(istp)=ah+al
detfer(istp)=ferr
enddo ! end of detection-metric, snr, and frequency error estimation loop
call indexx(detmet,nstep,indices) !find median of detection metric vector
xmed=detmet(indices(nstep/2))
detmet=detmet/xmed ! noise floor of detection metric is 1.0
ndet=0
do ip=1,20 ! use something like the "clean" algorithm to find candidates
iloc=maxloc(detmet)
il=iloc(1)
if( (detmet(il) .lt. 2.0) .or. (abs(detfer(il)) .gt. 100.0) ) cycle
ndet=ndet+1
times(ndet)=((il-1)*256+NPTS/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=10.0*log10(detmet(il))/2-5.0 !/2 because detmet is a 4th order moment
detmet(il-3:il+3)=0.0
enddo
nmessages=0
do ip=1,ndet
imid=times(ip)*fs
t0=times(ip)
cdat=cbig(imid-NPTS/2+1:imid+NPTS/2)
ferr=ferrs(ip)
nsnr=snrs(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
cc1=0
cc2=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)
! Find 5 largest peaks
do ipk=1,5
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,5
! we want ic to be the index of the first sample of the frame
ic0=ipeaks(ipk)
! fine adjustment of sync index
! bb lag used to place the sampling index at the center of the eye
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))
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
! Remove fine frequency error
call tweak1(cdat,NPTS,-ferr2,cdat2)
! place the beginning of frame at index NSPM+1
cdat2=cshift(cdat2,ic-(NSPM+1))
do iav=1,7 ! try each of 7 averaging patterns, hope that one works
if( iav .eq. 1 ) then
c(1:NSPM)=cdat2(NSPM+1:2*NSPM) !avg 1 frame to the right of ic
elseif( iav .eq. 2 ) then
c=cdat2(NSPM-431:NSPM+432) !1 frame centered on ic
c=cshift(c,-432)
elseif( iav .eq. 3 ) then !1 frame to the left of ic
c=cdat2(1:NSPM)
elseif( iav .eq. 4 ) then
c=cdat2(NSPM+432:NSPM+432+863) !1 frame beginning 36ms to the right of ic
c=cshift(c,432)
elseif( iav .eq. 5 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)
elseif( iav .eq. 6 ) then
c=cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
elseif( iav .eq. 7 ) then
c=cdat2(1:NSPM)+cdat2(NSPM+1:2*NSPM)+cdat2(2*NSPM+1:NPTS)
endif
! Estimate final frequency error and carrier phase.
cca=sum(c(1:1+41)*conjg(cb))
ccb=sum(c(1+56*6:1+56*6+41)*conjg(cb))
cfac=ccb*conjg(cca)
ffin=atan2(imag(cfac),real(cfac))/(twopi*56*6*dt)
phase0=atan2(imag(cca+ccb),real(cca+ccb))
! Remove phase error - want constellation such that sample points lie on re,im axes
cfac=cmplx(cos(phase0),sin(phase0))
c=c*conjg(cfac)
! sample to get softsamples
! do i=1,72
! softbits(2*i-1)=imag(c(1+(i-1)*12))
! softbits(2*i)=real(c(7+(i-1)*12))
! enddo
! matched filter - might be possible to improve this
softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(864-5:864))*pp(1:6))
softbits(2)=sum(real(c(1:12))*pp)
do i=2,72
softbits(2*i-1)=sum(imag(c(1+(i-1)*12-6:1+(i-1)*12+5))*pp)
softbits(2*i)=sum(real(c(7+(i-1)*12-6:7+(i-1)*12+5))*pp)
enddo
hardbits=0
do i=1,144
if( softbits(i) .ge. 0.0 ) then
hardbits(i)=1
endif
enddo
nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8 ) )/2
nbadsync2=(8-sum( (2*hardbits(1+56:8+56)-1)*s8 ) )/2
nbadsync=nbadsync1+nbadsync2
if( nbadsync .gt. 4 ) cycle
! normalize the softsymbols before submitting to decoder
sav=sum(softbits)/144
s2av=sum(softbits*softbits)/144
ssig=sqrt(s2av-sav*sav)
softbits=softbits/ssig
sigma=0.65
lratio(1:48)=softbits(9:9+47)
lratio(49:128)=softbits(65:65+80-1)
lratio=exp(2.0*lratio/(sigma*sigma))
unscrambledsoftbits(1:127:2)=lratio(1:64)
unscrambledsoftbits(2:128:2)=lratio(65:128)
max_iterations=20
max_dither=50
call ldpc_decode(unscrambledsoftbits, decoded, &
max_iterations, niterations, max_dither, ndither)
if( niterations .ge. 0.0 ) then
goto 778
endif
enddo ! frame averaging loop
enddo ! sample-time dither loop
enddo ! peak loop - could be made more efficient
msgreceived=' '
phase0=-98
i1hashdec=0
goto 999
778 continue
! The decoder found a codeword - compare decoded hash with calculated
! Collapse 80 decoded bits to 10 bytes. Bytes 1-9 are the message, byte 10 is the hash
do ibyte=1,10
itmp=0
do ibit=1,8
itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*8+ibit))
enddo
i1Dec8BitBytes(ibyte)=itmp
enddo
! Calculate the hash using the first 9 bytes.
ihashdec=nhash(c_loc(i1Dec8BitBytes),int(9,c_size_t),146)
ihashdec=2*iand(ihashdec,255)
! Compare calculated hash with received byte 10 - if they agree, keep the message.
i1hashdec=ihashdec
if( i1hashdec .eq. i1Dec8BitBytes(10) ) then
! Good hash --- unpack 72-bit message
do ibyte=1,12
itmp=0
do ibit=1,6
itmp=ishft(itmp,1)+iand(1,decoded((ibyte-1)*6+ibit))
enddo
i4Dec6BitWords(ibyte)=itmp
enddo
call unpackmsg(i4Dec6BitWords,msgreceived)
else
msgreceived=' '
phase0=-99
endif
if( msgreceived .ne. ' ' ) then
ndupe=0
do im=1,nmessages
if( allmessages(im) .eq. msgreceived ) ndupe=1
enddo
if( ndupe .eq. 0 ) then
nmessages=nmessages+1
allmessages(nmessages)=msgreceived
write(lines(nmessages),1020) nutc,nsnr,t0,nint(fest),msgreceived,char(0)
1020 format(i6.6,i4,f5.1,i5,' & ',a22,a1)
endif
endif
999 continue
! write(78,1001) nutc,t0,iav,ipk,is,fest,ffin,nbadsync1,nbadsync2, &
! phase0,niterations,ndither,i1hashdec,i1Dec8BitBytes(10),msgreceived
!1001 format(i6,f8.2,i4,i4,i4,f8.2,f8.2,i4,i4,f8.2,i4,i4,i4,i4,2x,a22)
enddo
return
end subroutine detectmsk144