305 lines
		
	
	
		
			7.9 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
		
		
			
		
	
	
			305 lines
		
	
	
		
			7.9 KiB
		
	
	
	
		
			Fortran
		
	
	
	
	
	
|   | subroutine syncmsk(cdat,npts,jpk,ipk,idf,rmax,snr,metric,decoded)
 | ||
|  | 
 | ||
|  | ! Attempt synchronization, and if successful decode using Viterbi algorithm.
 | ||
|  | 
 | ||
|  |   use iso_c_binding, only: c_loc,c_size_t
 | ||
|  |   use packjt
 | ||
|  |   use hashing
 | ||
|  |   use timer_module, only: timer
 | ||
|  | 
 | ||
|  |   parameter (NSPM=1404,NSAVE=2000)
 | ||
|  |   complex cdat(npts)                    !Analytic signal
 | ||
|  |   complex cb(66)                        !Complex waveform for Barker-11 code
 | ||
|  |   complex cd(0:11,0:3)
 | ||
|  |   complex c(0:NSPM-1)                   !Complex data for one message length
 | ||
|  |   complex c2(0:NSPM-1)
 | ||
|  |   complex cb3(1:NSPM,3)
 | ||
|  |   real r(12000)
 | ||
|  |   real rdat(12000)
 | ||
|  |   real ss1(12000)
 | ||
|  |   real symbol(234)
 | ||
|  |   real rdata(198)
 | ||
|  |   real rd2(198)
 | ||
|  |   real rsave(NSAVE)
 | ||
|  |   real xp(29)
 | ||
|  |   complex z,z0,z1,z2,z3,cfac
 | ||
|  |   integer*1 e1(198)
 | ||
|  |   integer*1, target :: d8(13)
 | ||
|  |   integer*1 i1hash(4)
 | ||
|  |   integer*1 i1
 | ||
|  |   integer*4 i4Msg6BitWords(12)            !72-bit message as 6-bit words
 | ||
|  |   integer mettab(0:255,0:1)               !Metric table for BPSK modulation
 | ||
|  |   integer ipksave(NSAVE)
 | ||
|  |   integer jpksave(NSAVE)
 | ||
|  |   integer indx(NSAVE)
 | ||
|  |   integer b11(11)                      !Barker-11 code
 | ||
|  |   character*22 decoded
 | ||
|  |   character*72 c72
 | ||
|  |   logical first
 | ||
|  |   equivalence (i1,i4)
 | ||
|  |   equivalence (ihash,i1hash)
 | ||
|  |   data xp/0.500000, 0.401241, 0.309897, 0.231832, 0.168095,    &
 | ||
|  |           0.119704, 0.083523, 0.057387, 0.039215, 0.026890,    &
 | ||
|  |           0.018084, 0.012184, 0.008196, 0.005475, 0.003808,    &
 | ||
|  |           0.002481, 0.001710, 0.001052, 0.000789, 0.000469,    &
 | ||
|  |           0.000329, 0.000225, 0.000187, 0.000086, 0.000063,    &
 | ||
|  |           0.000017, 0.000091, 0.000032, 0.000045/
 | ||
|  |   data first/.true./
 | ||
|  |   data b11/1,1,1,0,0,0,1,0,0,1,0/
 | ||
|  |   save first,cb,cd,twopi,dt,f0,f1,mettab
 | ||
|  | 
 | ||
|  |   phi=0.
 | ||
|  |   if(first) then
 | ||
|  | ! Get the metric table
 | ||
|  |      bias=0.0
 | ||
|  |      scale=20.0
 | ||
|  |      xln2=log(2.0)
 | ||
|  |      mettab=0
 | ||
|  |      do i=128,156
 | ||
|  |         x0=log(max(0.001,2.0*xp(i-127)))/xln2
 | ||
|  |         x1=log(max(0.001,2.0*(1.0-xp(i-127))))/xln2
 | ||
|  |         mettab(i,0)=nint(scale*(x0-bias))
 | ||
|  |         mettab(i,1)=nint(scale*(x1-bias))
 | ||
|  |         mettab(256-i,0)=mettab(i,1)
 | ||
|  |         mettab(256-i,1)=mettab(i,0)
 | ||
|  |      enddo
 | ||
|  |      do i=157,255
 | ||
|  |         mettab(i,0)=mettab(156,0)
 | ||
|  |         mettab(i,1)=mettab(156,1)
 | ||
|  |         mettab(256-i,0)=mettab(i,1)
 | ||
|  |         mettab(256-i,1)=mettab(i,0)
 | ||
|  |      enddo
 | ||
|  |      j=0
 | ||
|  |      twopi=8.0*atan(1.0)
 | ||
|  |      dt=1.0/12000.0
 | ||
|  |      f0=1000.0
 | ||
|  |      f1=2000.0
 | ||
|  |      dphi=0
 | ||
|  |      do i=1,11
 | ||
|  |         if(b11(i).eq.0) dphi=twopi*f0*dt
 | ||
|  |         if(b11(i).eq.1) dphi=twopi*f1*dt
 | ||
|  |         do n=1,6
 | ||
|  |            j=j+1
 | ||
|  |            phi=phi+dphi
 | ||
|  |            cb(j)=cmplx(cos(phi),sin(phi))
 | ||
|  |         enddo
 | ||
|  |      enddo
 | ||
|  |      cb3=0.
 | ||
|  |      cb3(1:66,1)=cb
 | ||
|  |      cb3(283:348,1)=cb
 | ||
|  |      cb3(769:834,1)=cb
 | ||
|  | 
 | ||
|  |      cb3(1:66,2)=cb
 | ||
|  |      cb3(487:552,2)=cb
 | ||
|  |      cb3(1123:1188,2)=cb
 | ||
|  | 
 | ||
|  |      cb3(1:66,3)=cb
 | ||
|  |      cb3(637:702,3)=cb
 | ||
|  |      cb3(919:984,3)=cb
 | ||
|  | 
 | ||
|  |      phi=0.
 | ||
|  |      do n=0,3
 | ||
|  |         k=-1
 | ||
|  |         dphi=twopi*f0*dt
 | ||
|  |         if(n.ge.2) dphi=twopi*f1*dt
 | ||
|  |         do i=0,5
 | ||
|  |            k=k+1
 | ||
|  |            phi=phi+dphi
 | ||
|  |            if(phi.gt.twopi) phi=phi-twopi
 | ||
|  |            cd(k,n)=cmplx(cos(phi),sin(phi))
 | ||
|  |         enddo
 | ||
|  | 
 | ||
|  |         dphi=twopi*f0*dt
 | ||
|  |         if(mod(n,2).eq.1) dphi=twopi*f1*dt
 | ||
|  |         do i=6,11
 | ||
|  |            k=k+1
 | ||
|  |            phi=phi+dphi
 | ||
|  |            if(phi.gt.twopi) phi=phi-twopi
 | ||
|  |            cd(k,n)=cmplx(cos(phi),sin(phi))
 | ||
|  |         enddo
 | ||
|  |      enddo
 | ||
|  | 
 | ||
|  |      first=.false.
 | ||
|  |   endif
 | ||
|  | 
 | ||
|  |   nfft=NSPM
 | ||
|  |   jz=npts-nfft
 | ||
|  |   decoded="                      "
 | ||
|  |   ipk=0
 | ||
|  |   jpk=0
 | ||
|  |   metric=-9999
 | ||
|  |   r=0.
 | ||
|  | 
 | ||
|  |   call timer('sync1   ',0)
 | ||
|  |   do j=1,jz                               !Find the Barker-11 sync vectors
 | ||
|  |      z=0.
 | ||
|  |      ss=0.
 | ||
|  |      do i=1,66
 | ||
|  |         ss=ss + real(cdat(j+i-1))**2 + aimag(cdat(j+i-1))**2
 | ||
|  |         z=z + cdat(j+i-1)*conjg(cb(i))    !Signal matching Barker 11
 | ||
|  |      enddo
 | ||
|  |      ss=sqrt(ss/66.0)*66.0
 | ||
|  |      r(j)=abs(z)/(0.908*ss)               !Goodness-of-fit to Barker 11
 | ||
|  |      ss1(j)=ss
 | ||
|  |   enddo
 | ||
|  |   call timer('sync1   ',1)
 | ||
|  | 
 | ||
|  |   call timer('sync2   ',0)
 | ||
|  |   jz=npts-nfft
 | ||
|  |   rmax=0.
 | ||
|  | !  n1=35, n2=69, n3=94
 | ||
|  |   k=0
 | ||
|  |   do j=1,jz                               !Find best full-message sync
 | ||
|  |      if(ss1(j).lt.85.0) cycle
 | ||
|  |      r1=r(j) + r(j+282) + r(j+768)        ! 6*(12+n1) 6*(24+n1+n2)
 | ||
|  |      r2=r(j) + r(j+486) + r(j+1122)       ! 6*(12+n2) 6*(24+n2+n3)
 | ||
|  |      r3=r(j) + r(j+636) + r(j+918)        ! 6*(12+n3) 6*(24+n3+n1)
 | ||
|  |      if(r1.gt.rmax) then
 | ||
|  |         rmax=r1
 | ||
|  |         jpk=j
 | ||
|  |         ipk=1
 | ||
|  |      endif
 | ||
|  |      if(r2.gt.rmax) then
 | ||
|  |         rmax=r2
 | ||
|  |         jpk=j
 | ||
|  |         ipk=2
 | ||
|  |      endif
 | ||
|  |      if(r3.gt.rmax) then
 | ||
|  |         rmax=r3
 | ||
|  |         jpk=j
 | ||
|  |         ipk=3
 | ||
|  |      endif
 | ||
|  |      rrmax=max(r1,r2,r3)
 | ||
|  |      if(rrmax.gt.1.9) then
 | ||
|  |         k=min(k+1,NSAVE)
 | ||
|  |         if(r1.eq.rrmax) ipksave(k)=1
 | ||
|  |         if(r2.eq.rrmax) ipksave(k)=2
 | ||
|  |         if(r3.eq.rrmax) ipksave(k)=3
 | ||
|  |         jpksave(k)=j
 | ||
|  |         rsave(k)=rrmax
 | ||
|  |      endif
 | ||
|  |   enddo
 | ||
|  |   call timer('sync2   ',1)
 | ||
|  |   kmax=k
 | ||
|  | 
 | ||
|  |   call indexx(rsave,kmax,indx)
 | ||
|  | 
 | ||
|  |   call timer('sync3   ',0)
 | ||
|  |   do kk=1,kmax
 | ||
|  |      k=indx(kmax+1-kk)
 | ||
|  |      ipk=ipksave(k)
 | ||
|  |      jpk=jpksave(k)
 | ||
|  |      rmax=rsave(k)
 | ||
|  | 
 | ||
|  |      c=conjg(cb3(1:NSPM,ipk))*cdat(jpk:jpk+nfft-1)
 | ||
|  |      smax=0.
 | ||
|  |      dfx=0.
 | ||
|  |      idfbest=0
 | ||
|  |      do itry=1,25
 | ||
|  |         idf=itry/2
 | ||
|  |         if(mod(itry,2).eq.0) idf=-idf
 | ||
|  |         idf=4*idf
 | ||
|  |         twk=idf
 | ||
|  |         call tweak1(c,NSPM,-twk,c2)
 | ||
|  |         z=sum(c2)
 | ||
|  |         if(abs(z).gt.smax) then
 | ||
|  |            dfx=twk
 | ||
|  |            smax=abs(z)
 | ||
|  |            phi=atan2(aimag(z),real(z))            !Carrier phase offset
 | ||
|  |            idfbest=idf
 | ||
|  |         endif
 | ||
|  |      enddo
 | ||
|  |      idf=idfbest
 | ||
|  |      call tweak1(cdat,npts,-dfx,cdat)
 | ||
|  |      cfac=cmplx(cos(phi),-sin(phi))
 | ||
|  |      cdat=cfac*cdat
 | ||
|  | 
 | ||
|  |      sig=0.
 | ||
|  |      ref=0.
 | ||
|  |      rdat(1:npts)=cdat
 | ||
|  |      iz=11
 | ||
|  |      do k=1,234                                !Compute soft symbols
 | ||
|  |         j=jpk+6*(k-1)
 | ||
|  | 
 | ||
|  |         z0=2.0*dot_product(cdat(j:j+iz),cd(0:iz,0))
 | ||
|  |         z1=2.0*dot_product(cdat(j:j+iz),cd(0:iz,1))
 | ||
|  |         z2=2.0*dot_product(cdat(j:j+iz),cd(0:iz,2))
 | ||
|  |         z3=2.0*dot_product(cdat(j:j+iz),cd(0:iz,3))
 | ||
|  | 
 | ||
|  | !### Maybe these should be weighted by yellow() ?
 | ||
|  |         if(j+1404+iz.lt.npts) then
 | ||
|  |            z0=z0 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,0))
 | ||
|  |            z1=z1 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,1))
 | ||
|  |            z2=z2 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,2))
 | ||
|  |            z3=z3 + dot_product(cdat(j+1404:j+1404+iz),cd(0:iz,3))
 | ||
|  |         endif
 | ||
|  | 
 | ||
|  |         if(j-1404.ge.1) then
 | ||
|  |            z0=z0 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,0))
 | ||
|  |            z1=z1 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,1))
 | ||
|  |            z2=z2 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,2))
 | ||
|  |            z3=z3 + dot_product(cdat(j-1404:j-1404+iz),cd(0:iz,3))
 | ||
|  |         endif
 | ||
|  | 
 | ||
|  |         sym=max(abs(real(z2)),abs(real(z3))) - max(abs(real(z0)),abs(real(z1)))
 | ||
|  | 
 | ||
|  |         if(sym.lt.0.0) then
 | ||
|  |            phi=atan2(aimag(z0),real(z0))
 | ||
|  |            sig=sig + real(z0)**2
 | ||
|  |            ref=ref + aimag(z0)**2
 | ||
|  |         else
 | ||
|  |            phi=atan2(aimag(z1),real(z1))
 | ||
|  |            sig=sig + real(z1)**2
 | ||
|  |            ref=ref + aimag(z1)**2
 | ||
|  |         endif
 | ||
|  |         n=k
 | ||
|  |         if(ipk.eq.2) n=k+47
 | ||
|  |         if(ipk.eq.3) n=k+128
 | ||
|  |         if(n.gt.234) n=n-234
 | ||
|  |         ibit=0
 | ||
|  |         if(sym.ge.0) ibit=1
 | ||
|  |         symbol(n)=sym
 | ||
|  |      enddo
 | ||
|  |      snr=db(sig/ref-1.0)
 | ||
|  | 
 | ||
|  |      rdata(1:35)=symbol(12:46)
 | ||
|  |      rdata(36:104)=symbol(59:127)
 | ||
|  |      rdata(105:198)=symbol(140:233)
 | ||
|  | 
 | ||
|  | ! Re-order the symbols and make them i*1
 | ||
|  |      j=0
 | ||
|  |      do i=1,99
 | ||
|  |         i4=128+rdata(i)                       !### Should be nint() ??? ###
 | ||
|  |         if(i4.gt.255)  i4=255
 | ||
|  |         if(i4.lt.0) i4=0
 | ||
|  |         j=j+1
 | ||
|  |         e1(j)=i1
 | ||
|  |         rd2(j)=rdata(i)
 | ||
|  |         i4=128+rdata(i+99)
 | ||
|  |         if(i4.gt.255)  i4=255
 | ||
|  |         if(i4.lt.0) i4=0
 | ||
|  |         j=j+1
 | ||
|  |         e1(j)=i1
 | ||
|  |         rd2(j)=rdata(i+99)
 | ||
|  |      enddo
 | ||
|  | 
 | ||
|  | ! Decode the message
 | ||
|  |      nb1=87
 | ||
|  |      call vit213(e1,nb1,mettab,d8,metric)
 | ||
|  |      ihash=nhash(c_loc(d8),int(9,c_size_t),146)
 | ||
|  |      ihash=2*iand(ihash,32767)
 | ||
|  |      decoded='                      '
 | ||
|  |      if(d8(10).eq.i1hash(2) .and. d8(11).eq.i1hash(1)) then
 | ||
|  |         write(c72,1012) d8(1:9)
 | ||
|  | 1012    format(9b8.8)
 | ||
|  |         read(c72,1014) i4Msg6BitWords
 | ||
|  | 1014    format(12b6.6)
 | ||
|  |         call unpackmsg(i4Msg6BitWords,decoded,.false.,'      ') !Unpack to get msgsent
 | ||
|  |      endif
 | ||
|  |      if(decoded.ne.'                      ') exit
 | ||
|  |   enddo
 | ||
|  |   call timer('sync3   ',1)
 | ||
|  | 
 | ||
|  |   return
 | ||
|  | end subroutine syncmsk
 |