WSJT-X/lib/detectmsk32.f90
Steven Franke aeed9e3344 First-cut at decoder for (32,16) msk32. Needs more work.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6954 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2016-07-27 14:40:38 +00:00

227 lines
6.0 KiB
Fortran

subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00)
use timer_module, only: timer
parameter (NSPM=192, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=10)
character*4 rpt(0:63)
character*6 mycall,partnercall
character*22 msg,msgsent,msgreceived,allmessages(32)
character*80 lines(100)
complex bb(6)
complex cbig(n)
complex cdat(NPTS) !Analytic signal
complex ctmp(NPTS) !Analytic signal
complex cft(512)
complex cwaveforms(192,64)
integer, dimension(1) :: iloc
integer indices(MAXSTEPS)
integer itone(144)
logical ismask(NFFT)
real detmet(-2:MAXSTEPS+3)
real detfer(MAXSTEPS)
real ferrs(MAXCAND)
real pp(12)
real rcw(12)
real snrs(MAXCAND)
real times(MAXCAND)
real tonespec(NFFT)
real*8 dt, df, fs, pi, twopi
logical first
data first/.true./
save df,first,cb,cbr,fs,nhashes,pi,twopi,dt,rcw,pp,nmatchedfilter,cwaveforms,rpt
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
do i=0,30
if( i.lt.5 ) then
write(rpt(i),'(a1,i2.2,a1)') '-',abs(i-5)
write(rpt(i+31),'(a2,i2.2,a1)') 'R-',abs(i-5)
else
write(rpt(i),'(a1,i2.2,a1)') '+',i-5
write(rpt(i+31),'(a2,i2.2,a1)') 'R+',i-5
endif
enddo
rpt(62)='RRR '
rpt(63)='73 '
dphi0=twopi*(freq-500)/12000.0
dphi1=twopi*(freq+500)/12000.0
do i=1,64
msg='<'//trim(mycall)//' '//trim(partnercall)//'> '//rpt(i-1)
call genmsk32(msg,msgsent,0,itone,itype)
! write(*,*) i,msg,msgsent,itype
nsym=32
phi=0.0
indx=1
nreps=1
do jrep=1,nreps
do isym=1,nsym
if( itone(isym) .eq. 0 ) then
dphi=dphi0
else
dphi=dphi1
endif
do j=1,6
cwaveforms(indx,i)=cmplx(cos(phi),sin(phi));
indx=indx+1
phi=mod(phi+dphi,twopi)
enddo
enddo
enddo
enddo
first=.false.
endif
! Fill the detmet, detferr arrays
nstepsize=48 ! 4ms steps
nstep=(n-NPTS)/nstepsize
detmet=0
detmax=-999.99
detfer=-999.99
do istp=1,nstep
ns=1+nstepsize*(istp-1)
ne=ns+NPTS-1
if( ne .gt. n ) exit
ctmp=cmplx(0.0,0.0)
ctmp(1:NPTS)=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(NPTS-11:NPTS)=ctmp(NPTS-11:NPTS)*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)
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)
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)
detfer(istp)=ferr
! write(*,*) istp,ilpk,ihpk,ah,al
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 i=1,nstep
!write(77,*) i,detmet(i),detfer(i)
!enddo
do ip=1,MAXCAND ! find candidates
iloc=maxloc(detmet(1:nstep))
il=iloc(1)
if( (detmet(il) .lt. 4.2) ) exit
if( abs(detfer(il)) .le. ntol ) then
ndet=ndet+1
times(ndet)=((il-1)*nstepsize+NPTS/2)*dt
ferrs(ndet)=detfer(il)
snrs(ndet)=12.0*log10(detmet(il)-1)/2-8.0
endif
detmet(max(1,il-3):min(nstep,il+3))=0.0
! detmet(il)=0.0
enddo
nmessages=0
lines=char(0)
fbest=1e6
pkbest=-1e6
imsgbest=-1
istartbest=-1
ipbest=-1
nsnrbest=-100
t0best=-1e6
do ip=1,ndet !run through the candidates and try to sync/demod/decode
imid=times(ip)*fs
if( imid .lt. NPTS/2 ) imid=NPTS/2
if( imid .gt. n-NPTS/2 ) imid=n-NPTS/2
t0=times(ip) + t00
cdat=cbig(imid-NPTS/2+1:imid+NPTS/2)
ferr=ferrs(ip)
nsnr=2*nint(snrs(ip)/2.0)
if( nsnr .lt. -4 ) nsnr=-4
if( nsnr .gt. 24 ) nsnr=24
do imsg=1,64
do istart=NSPM-NSPM/2,NPTS-NSPM
cft(1:144)=cdat(istart:istart+144-1)*conjg(cwaveforms(1:144,imsg))
cft(145:512)=0.
df=12000.0/512.0
call four2a(cft,512,1,-1,1)
iloc=maxloc(abs(cft))
ipk=iloc(1)
pk=abs(cft(ipk))
fpk=(ipk-1)*df
if( fpk.gt.12000.0 ) fpk=fpk-12000.0
if( pk .gt. pkbest .and. abs(fpk-1500.0) .le. ntol) then
ipbest=ip
pkbest=pk
fbest=fpk
imsgbest=imsg
istartbest=istart
nsnrbest=nsnr
t0best=t0
endif
enddo
enddo
enddo ! candidate loop
999 continue
msgreceived=' '
if( imsgbest .gt. 0 .and. pkbest .ge. 108.0) then
nrxrpt=iand(imsgbest-1,63)
nrxhash=(imsgbest-1-nrxrpt)/64
!write(*,*) ipbest,pkbest,fbest,imsgbest,istartbest,nsnrbest,t0best,nrxrpt,nrxhash
nmessages=1
write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall), &
trim(partnercall),">",rpt(nrxrpt)
write(lines(nmessages),1020) nutc,nsnrbest,t0best,nint(fbest),msgreceived
1020 format(i6.6,i4,f5.1,i5,' & ',a22)
endif
return
end subroutine detectmsk32