subroutine opdetmsk144(cbig,n,lines,nmessages,nutc,ntol,t00) use timer_module, only: timer ! NSPM: number of samples per message frame ! NAVG: number of frames to coherently average parameter (NSPM=864, NAVG=7, NPTS=NAVG*NSPM, NSTEP=6000) character*22 msgreceived,allmessages(20) character*80 lines(100) complex cbig(n) complex cdat(NPTS) !Analytic signal complex cdat2(NPTS) complex c(NSPM) complex cr(NSPM,NAVG) complex ct(NSPM) complex cs(NSPM) complex cb(42) !Complex waveform for sync word complex cfac,cca,ccb complex cc(0:NSPM-1) complex csum integer s8(8),hardbits(144) integer, dimension(1) :: iloc integer*1 decoded(80) integer ipeaks(10) real cbi(42),cbq(42) real rcw(12) real ccm(0:NSPM-1) real ccms(0:NSPM-1) real pp(12) !Half-sine pulse shape real*8 dt, fs, pi, twopi real softbits(144) real*8 lratio(128) real llr(128) logical first data first/.true./ data s8/0,1,1,1,0,0,1,0/ save first,cb,fs,pi,twopi,dt,s8,rcw,pp,nmatchedfilter 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 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 waveforms 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) first=.false. endif nmessages=0 allmessages=char(0) lines=char(0) nshort=0 trec=NPTS/12000.0 ! Duration of the data record ! write(*,*) "number of points in opdetmsk144",n if( n .lt. NPTS .or. n .gt. 181000) return nsteps=2*n/NSTEP-1 ! Hardwired for 0.5s steps ! write(*,*) 'nsteps ',nsteps nsnr=-4 do istep=1,nsteps ib=(istep-1)*NPTS/2+1 ie=ib+NPTS-1 if( ie .gt. n ) then ie=n ib=n-NPTS+1 endif t0=t00+(ib-1)/12000.0+trec/2 cdat=cbig(ib:ie) xmax=0.0 bestf=0.0 do ifr=-ntol,ntol ! search for frequency that maximizes sync correlation ferr=ifr ! call timer('opdetavg ',0) call tweak1(cdat,NPTS,-(1500+ferr),cdat2) cr=reshape(cdat2,(/NSPM,NAVG/)) c=sum(cr,2) ! call timer('opdetavg ',1) ! call timer('opdetsync',0) cc=0 do ish=0,NSPM-1 ct=cshift(c,ish) ! cc(ish)=sum(ct(1:42)*conjg(cb))+sum(ct(56*6:56*6+41)*conjg(cb)) ! cc(ish)=sum(ct(1:42)*conjg(cb)+ ct(56*6:56*6+41)*conjg(cb)) csum=0 do j=1,42 csum=csum+(ct(j)+ct(56*6+j-1))*conjg(cb(j)) enddo cc(ish)=csum enddo ccm=abs(cc) xb=maxval(ccm) ! call timer('opdetsync',1) if( xb .gt. xmax ) then xmax=xb bestf=ferr cs=c ccms=ccm endif enddo fest=1500+bestf ! write(*,*) istep,fest,xmax c=cs ccm=ccms ! Find 2 largest peaks do ipk=1, 2 iloc=maxloc(ccm) ic2=iloc(1) ipeaks(ipk)=ic2 ccm(max(0,ic2-7):min(NSPM-1,ic2+7))=0.0 enddo do ipk=1,2 do is=1,3 ic0=ipeaks(ipk) if( is.eq.2 ) ic0=max(1,ic0-1) if( is.eq.3 ) ic0=min(NSPM,ic0+1) ct=cshift(c,ic0-1) ! Estimate final frequency error and carrier phase. cca=sum(ct(1:1+41)*conjg(cb)) ccb=sum(ct(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 rotated so that sample points lie on I/Q axes cfac=cmplx(cos(phase0),sin(phase0)) ct=ct*conjg(cfac) if( nmatchedfilter .eq. 0 ) then ! 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 else ! matched filter - softbits(1)=sum(imag(ct(1:6))*pp(7:12))+sum(imag(ct(864-5:864))*pp(1:6)) softbits(2)=sum(real(ct(1:12))*pp) do i=2,72 softbits(2*i-1)=sum(imag(ct(1+(i-1)*12-6:1+(i-1)*12+5))*pp) softbits(2*i)=sum(real(ct(7+(i-1)*12-6:7+(i-1)*12+5))*pp) enddo endif ! sync word hard error weight is a good discriminator for ! frames that have reasonable probability of decoding 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 !write(*,*) 'peak index ',ipk,' pk loc: ',ic0,' nbadsync ',nbadsync 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.75 lratio(1:48)=softbits(9:9+47) lratio(49:128)=softbits(65:65+80-1) llr=2.0*lratio/(sigma*sigma) lratio=exp(2.0*lratio/(sigma*sigma)) max_iterations=10 max_dither=1 call timer('bpdec144 ',0) call bpdecode144(llr,max_iterations,decoded,niterations) call timer('bpdec144 ',1) if( niterations .ge. 0.0 ) then call extractmessage144(decoded,msgreceived,nhashflag) if( nhashflag .gt. 0 ) then ! CRCs match, so print it 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 1020 format(i6.6,i4,f5.1,i5,' ^ ',a22) endif goto 999 else msgreceived=' ' ndither=-99 ! -99 is bad hash flag endif endif enddo ! slicer dither enddo ! peak loop enddo msgreceived=' ' ndither=-98 999 continue return end subroutine opdetmsk144