WSJT-X/lib/msk144sd.f90

198 lines
5.1 KiB
Fortran

program msk144sd
!
! A simple decoder for slow msk144.
! Can be used as a (slow) brute-force multi-decoder by looping
! over a set of carrier frequencies.
!
use options
use timer_module, only: timer
use timer_impl, only: init_timer
use readwav
parameter (NRECENT=10)
parameter (NSPM=864)
parameter (NPATTERNS=4)
character ch
character*80 line
character*500 infile
character*12 mycall,hiscall
character*6 mygrid
character(len=500) optarg
character*22 msgreceived
character*12 recent_calls(NRECENT)
complex cdat(30*375)
complex c(NSPM)
complex ct(NSPM)
real softbits(144)
real xmc(NPATTERNS)
logical :: display_help=.false.
type(wav_header) :: wav
integer iavmask(8)
integer iavpatterns(8,NPATTERNS)
integer npkloc(10)
integer*2 id2(30*12000)
integer*2 ichunk(7*1024)
data iavpatterns/ &
1,1,1,1,0,0,0,0, &
0,1,1,1,1,0,0,0, &
0,0,1,1,1,1,0,0, &
1,1,1,1,1,1,0,0/
data xmc/2.0,4.5,2.5,3.0/
type (option) :: long_options(2) = [ &
option ('frequency',.true.,'f','rxfreq',''), &
option ('help',.false.,'h','Display this help message','') &
]
t0=0.0
ntol=100
nrxfreq=1500
do
call getopt('f:h',long_options,ch,optarg,narglen,nstat,noffset,nremain,.true.)
if( nstat .ne. 0 ) then
exit
end if
select case (ch)
case ('f')
read (optarg(:narglen), *) nrxfreq
case ('h')
display_help = .true.
end select
end do
if(display_help .or. nstat.lt.0 .or. nremain.lt.1) then
print *, ''
print *, 'Usage: msk144sd [OPTIONS] file1 [file2 ...]'
print *, ''
print *, ' decode pre-recorded .WAV file(s)'
print *, ''
print *, 'OPTIONS:'
do i = 1, size (long_options)
call long_options(i) % print (6)
end do
go to 999
endif
call init_timer ('timer.out')
call timer('msk144 ',0)
ndecoded=0
do ifile=noffset+1,noffset+nremain
call get_command_argument(ifile,optarg,narglen)
infile=optarg(:narglen)
call timer('read ',0)
call wav%read (infile)
i1=index(infile,'.wav')
if( i1 .eq. 0 ) i1=index(infile,'.WAV')
read(infile(i1-6:i1-1),*,err=998) nutc
inquire(FILE=infile,SIZE=isize)
npts=min((isize-216)/2,360000)
read(unit=wav%lun) id2(1:npts)
close(unit=wav%lun)
call timer('read ',1)
! do if=1,89 ! brute force multi-decoder
fo=nrxfreq
! fo=(if-1)*25.0+300.0
call msksddc(id2,npts,fo,cdat)
np=npts/32
ntol=200 ! actual ntol is ntol/32=6.25 Hz. Detection window is 12.5 Hz wide
fc=1500.0
call msk144spd(cdat,np,ntol,ndecodesuccess,msgreceived,fc,fest,tdec,navg,ct, &
softbits,recent_calls,nrecent)
nsnr=0 ! need an snr estimate
if( ndecodesuccess .eq. 1 ) then
fest=fo+fest-fc ! fudging because spd thinks input signal is at 1500 Hz
goto 900
endif
! If short ping decoder doesn't find a decode
npat=NPATTERNS
do iavg=1,npat
iavmask=iavpatterns(1:8,iavg)
navg=sum(iavmask)
deltaf=4.0/real(navg) ! search increment for frequency sync
npeaks=4
ntol=200
fc=1500.0
call msk144sync(cdat(1:6*NSPM),6,ntol,deltaf,iavmask,npeaks,fc, &
fest,npkloc,nsyncsuccess,xmax,c)
if( nsyncsuccess .eq. 0 ) cycle
do ipk=1,npeaks
do is=1,3
ic0=npkloc(ipk)
if(is.eq.2) ic0=max(1,ic0-1)
if(is.eq.3) ic0=min(NSPM,ic0+1)
ct=cshift(c,ic0-1)
call msk144decodeframe(ct,softbits,msgreceived,ndecodesuccess, &
recent_calls,nrecent)
if(ndecodesuccess .gt. 0) then
tdec=tsec+xmc(iavg)*tframe
fest=fo+(fest-fc)/32.0
goto 900
endif
enddo !Slicer dither
enddo !Peak loop
enddo
! enddo
900 continue
if( ndecodesuccess .gt. 0 ) then
write(*,1020) nutc,nsnr,tdec,nint(fest),' % ',msgreceived,navg
1020 format(i6.6,i4,f5.1,i5,a3,a22,i4)
endif
enddo
call timer('msk144 ',1)
call timer('msk144 ',101)
go to 999
998 print*,'Cannot read from file:'
print*,infile
999 continue
end program msk144sd
subroutine msksddc(id2,npts,fc,cdat)
! The msk144 detector/demodulator/decoder will decode signals
! with carrier frequency, fc, in the range fN/4 +/- 0.03333*fN.
!
! For slow MSK144 with nslow=32:
! fs=12000/32=375 Hz, fN=187.5 Hz
!
! This routine accepts input samples with fs=12000 Hz. It
! downconverts and decimates by 32 to center a signal with input carrier
! frequency fc at new carrier frequency 1500/32=46.875 Hz.
! The analytic signal is returned.
parameter (NFFT1=30*12000,NFFT2=30*375)
integer*2 id2(npts)
complex cx(0:NFFT1)
complex cdat(30*375)
dt=1.0/12000.0
df=1.0/(NFFT1*dt)
icenter=int(fc/df+0.5)
i46p875=int(46.875/df+0.5)
ishift=icenter-i46p875
cx=cmplx(0.0,0.0)
cx(1:npts)=id2
call four2a(cx,NFFT1,1,-1,1)
cx=cshift(cx,ishift)
cx(1)=0.5*cx(1)
cx(2*i46p875+1:)=cmplx(0.0,0.0)
call four2a(cx,NFFT2,1,1,1)
cdat(1:npts/32)=cx(0:npts/32-1)/NFFT1
return
end subroutine msksddc