diff --git a/CMakeLists.txt b/CMakeLists.txt index ccbc900e8..a37f230b2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -408,10 +408,14 @@ set (wsjt_FSRCS lib/move.f90 lib/msk144d.f90 lib/msk144d2.f90 + lib/msk40decodeframe.f90 lib/msk144decodeframe.f90 + lib/msk40spd.f90 lib/msk144spd.f90 + lib/msk40sync.f90 lib/msk144sync.f90 lib/msk144_decode.f90 + lib/msk40_freq_search.f90 lib/msk144_freq_search.f90 lib/mskrtd.f90 lib/msk144sim.f90 diff --git a/lib/msk40_freq_search.f90 b/lib/msk40_freq_search.f90 new file mode 100644 index 000000000..2515fc305 --- /dev/null +++ b/lib/msk40_freq_search.f90 @@ -0,0 +1,52 @@ +subroutine msk40_freq_search(cdat,fc,if1,if2,delf,nframes,navmask,cb, & + cdat2,xmax,bestf,cs,xccs) + + parameter (NSPM=240,NZ=7*NSPM) + complex cdat(NZ) + complex cdat2(NZ) + 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) + real xcc(0:NSPM-1) + real xccs(0:NSPM-1) + integer navmask(nframes) !Tells which frames to average + + navg=sum(navmask) + n=nframes*NSPM + fac=1.0/(48.0*sqrt(float(navg))) + + do ifr=if1,if2 !Find freq that maximizes sync + ferr=ifr*delf + call tweak1(cdat,n,-(fc+ferr),cdat2) + c=0 + do i=1,nframes + ib=(i-1)*NSPM+1 + ie=ib+NSPM-1 + if( navmask(i) .eq. 1 ) c=c+cdat2(ib:ie) + 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),cb(1:42)) + enddo + + xcc=abs(cc) + xb=maxval(xcc)*fac + if(xb.gt.xmax) then + xmax=xb + bestf=ferr + cs=c + xccs=xcc + endif + enddo + +! write(71,3001) fc,delf,if1,if2,nframes,bestf,xmax +!3001 format(2f8.3,3i5,2f8.3) + + return +end subroutine msk40_freq_search diff --git a/lib/msk40decodeframe.f90 b/lib/msk40decodeframe.f90 new file mode 100644 index 000000000..ebd515942 --- /dev/null +++ b/lib/msk40decodeframe.f90 @@ -0,0 +1,140 @@ +subroutine msk40decodeframe(c,mycall,hiscall,xsnr,msgreceived,nsuccess) +! use timer_module, only: timer + + parameter (NSPM=240) + character*4 rpt(0:15) + character*6 mycall,hiscall,mycall0,hiscall0 + character*22 hashmsg,msgreceived + complex cb(42) + complex cfac,cca + complex c(NSPM) + integer*1 cw(32) + integer*1 decoded(16) + integer s8r(8),hardbits(40) + integer nhashes(0:15) + real*8 dt, fs, pi, twopi + real cbi(42),cbq(42) + real pp(12) + real softbits(40) + real llr(32) + logical first + data first/.true./ + data s8r/1,0,1,1,0,0,0,1/ + data mycall0/'dummy'/,hiscall0/'dummy'/ + data rpt/"-03 ","+00 ","+03 ","+06 ","+10 ","+13 ","+16 ", & + "R-03","R+00","R+03","R+06","R+10","R+13","R+16", & + "RRR ","73 "/ + save first,cb,fs,nhashes,pi,twopi,dt,s8r,pp,rpt,mycall0,hiscall0 + + if(first) then +! define half-sine pulse and raised-cosine edge window + pi=4d0*datan(1d0) + twopi=8d0*datan(1d0) + fs=12000.0 + dt=1.0/fs + + do i=1,12 + angle=(i-1)*pi/12.0 + pp(i)=sin(angle) + enddo + +! define the sync word waveforms + 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) + cb=cmplx(cbi,cbq) + first=.false. + endif + + if(mycall.ne.mycall0 .or. hiscall.ne.hiscall0) then + do i=0,15 + hashmsg=trim(mycall)//' '//trim(hiscall)//' '//rpt(i) + call fmtmsg(hashmsg,iz) + call hash(hashmsg,22,ihash) + nhashes(i)=iand(ihash,4095) + enddo + mycall0=mycall + hiscall0=hiscall + endif + + nsuccess=0 + msgreceived=' ' + +! Estimate carrier phase. + cca=sum(c(1:1+41)*conjg(cb)) + phase0=atan2(imag(cca),real(cca)) + +! Remove phase error - want constellation rotated so that sample points lie on I/Q axes + cfac=cmplx(cos(phase0),sin(phase0)) + c=c*conjg(cfac) + +! matched filter - + softbits(1)=sum(imag(c(1:6))*pp(7:12))+sum(imag(c(NSPM-5:NSPM))*pp(1:6)) + softbits(2)=sum(real(c(1:12))*pp) + do i=2,20 + 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 + +! sync word hard error weight is used as a discriminator for +! frames that have reasonable probability of decoding + hardbits=0 + do i=1,40 + if( softbits(i) .ge. 0.0 ) then + hardbits(i)=1 + endif + enddo + nbadsync1=(8-sum( (2*hardbits(1:8)-1)*s8r ) )/2 + nbadsync=nbadsync1 + if( nbadsync .gt. 3 ) then + return + endif + +! normalize the softsymbols before submitting to decoder + sav=sum(softbits)/40 + s2av=sum(softbits*softbits)/40 + ssig=sqrt(s2av-sav*sav) + softbits=softbits/ssig + + sigma=0.75 + if(xsnr.lt.0.0) sigma=0.75-0.0875*xsnr + llr(1:32)=softbits(9:40) + llr=2.0*llr/(sigma*sigma) + + max_iterations=5 + call bpdecode40(llr,max_iterations,decoded,niterations) + + if( niterations .ge. 0.0 ) then + call encode_msk40(decoded,cw) + nhammd=0 + cord=0.0 + do i=1,32 + if( cw(i) .ne. hardbits(i+8) ) then + nhammd=nhammd+1 + cord=cord+abs(softbits(i+8)) + endif + enddo + + imsg=0 + do i=1,16 + imsg=ishft(imsg,1)+iand(1,decoded(17-i)) + enddo + nrxrpt=iand(imsg,15) + nrxhash=(imsg-nrxrpt)/16 +!write(*,*) 'decodeframe ',nhammd,cord,nrxhash,nrxrpt,nhashes(nrxrpt) + if(nhammd.le.5 .and. cord .lt. 1.7 .and. nrxhash.eq.nhashes(nrxrpt)) then + nsuccess=1 + write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall), & + trim(hiscall),">",rpt(nrxrpt) + return + endif + endif + + return +end subroutine msk40decodeframe diff --git a/lib/msk40spd.f90 b/lib/msk40spd.f90 index a7714d3c9..7d891e343 100644 --- a/lib/msk40spd.f90 +++ b/lib/msk40spd.f90 @@ -1,9 +1,10 @@ -subroutine msk40spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret) +subroutine msk40spd(cbig,n,ntol,mycall,hiscall,nsuccess,msgreceived,fc,fret,tret) ! msk40 short-ping-decoder use timer_module, only: timer - parameter (NSPM=240, MAXSTEPS=100, NFFT=NSPM, MAXCAND=5, NPATTERNS=6) + parameter (NSPM=240, MAXSTEPS=150, NFFT=NSPM, MAXCAND=5, NPATTERNS=6) + character*6 mycall,hiscall character*22 msgreceived complex cbig(n) complex cdat(3*NSPM) !Analytic signal @@ -164,10 +165,10 @@ subroutine msk40spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret) endif cdat=cbig(ib:ie) fo=fc+ferrs(icand) + xsnr=snrs(icand) do iav=1,NPATTERNS navmask=navpatterns(1:3,iav) call msk40sync(cdat,3,ntol0,deltaf,navmask,npeaks,fo,fest,npkloc,nsyncsuccess,c) - if( nsyncsuccess .eq. 0 ) cycle do ipk=1,npeaks @@ -176,7 +177,7 @@ subroutine msk40spd(cbig,n,ntol,nsuccess,msgreceived,fc,fret,tret) if( is.eq.2) ic0=max(1,ic0-1) if( is.eq.3) ic0=min(NSPM,ic0+1) ct=cshift(c,ic0-1) - call msk40decodeframe(ct,msgreceived,ndecodesuccess) + call msk40decodeframe(ct,mycall,hiscall,xsnr,msgreceived,ndecodesuccess) if( ndecodesuccess .gt. 0 ) then tret=(nstart(icand)+NSPM/2)/fs diff --git a/lib/msk40sync.f90 b/lib/msk40sync.f90 index ab4e6872a..5a2ca1e77 100644 --- a/lib/msk40sync.f90 +++ b/lib/msk40sync.f90 @@ -10,7 +10,6 @@ subroutine msk40sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest, & complex cs(NSPM,8) complex cb(42) !Complex waveform for sync word -! integer*8 count0,count1,clkfreq integer s8r(8) integer iloc(1) integer npklocs(npeaks) @@ -25,10 +24,9 @@ subroutine msk40sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest, & real pp(12) !Half-sine pulse shape logical first data first/.true./ - data s8r/0,1,0,0,1,1,1,0/ - save first,cb,fs,pi,twopi,dt,s8,pp + data s8r/1,0,1,1,0,0,0,1/ + save first,cb,fs,pi,twopi,dt,s8r,pp -! call system_clock(count0,clkfreq) if(first) then pi=4.0*atan(1.0) twopi=8.0*atan(1.0) @@ -98,16 +96,8 @@ subroutine msk40sync(cdat,nframes,ntol,delf,navmask,npeaks,fc,fest, & xcc(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0 enddo - if( xmax .lt. 0.7 ) then - nsuccess=0 - else - nsuccess=1 - endif - -! call system_clock(count1,clkfreq) -! t=float(count1-count0)/clkfreq -! write(72,3001) nfreqs,OMP_GET_MAX_THREADS(),nthreads,t -!3001 format(3i6,f8.3) + nsuccess=0 + if( xmax .ge. 1.3 ) nsuccess=1 return end subroutine msk40sync diff --git a/lib/mskrtd.f90 b/lib/mskrtd.f90 index 43ceb8628..7b42dace4 100644 --- a/lib/mskrtd.f90 +++ b/lib/mskrtd.f90 @@ -13,6 +13,8 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,line) character*22 msgreceived !Decoded message character*22 msglast !!! temporary - used for dupechecking character*80 line !Formatted line with UTC dB T Freq Msg +!##### TEMPORARY + character*6 mycall,hiscall complex cdat(NFFT1) !Analytic signal complex c(NSPM) !Coherently averaged complex data @@ -34,6 +36,9 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,line) 1,1,1,1,1,0,0,0, & 1,1,1,1,1,1,1,0/ data xmc/2.0,4.5,2.5,3.5/ !Used to label decode with time at center of averaging mask +!###### TEMPORARY + data mycall/'K9AN'/ + data hiscall/'K1JT'/ save first,tsec0,nutc00,pnoise,nsnrlast,msglast,cdat @@ -85,6 +90,13 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,line) ! 3 frames along with 2- and 3-frame averages. np=8*NSPM call msk144spd(cdat,np,ntol,nsuccess,msgreceived,fc,fest,tdec) + +!############################################################ +!##### hardwired for testing - need to bring in Sh box status + if( nsuccess .eq. 0 .and. .false. ) then + call msk40spd(cdat,np,ntol,mycall,hiscall,nsuccess,msgreceived,fc,fest,tdec) + endif + if( nsuccess .eq. 1 ) then tdec=tsec+tdec decsym=' & ' @@ -125,6 +137,7 @@ subroutine mskrtd(id2,nutc0,tsec,ntol,nrxfreq,ndepth,line) enddo !Peak loop enddo + msgreceived=' ' ! no decode - update noise level used for calculating displayed snr.