WSJT-X/lib/syncmsk.f90

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