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,0,1,1,1,1,0,0, & 1,1,1,1,1,0,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=3 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(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