WSJT-X/lib/detectmsk32.f90

subroutine detectmsk32(cbig,n,mycall,partnercall,lines,nmessages,nutc,ntol,t00)
  use timer_module, only: timer

  parameter (NSPM=192, NPTS=3*NSPM, MAXSTEPS=7500, NFFT=3*NSPM, MAXCAND=10)
  character*4 rpt(0:63)
  character*6 mycall,partnercall
  character*22 msg,msgsent,msgreceived,allmessages(32)
  character*80 lines(100)
  complex bb(6)
  complex cbig(n)
  complex cdat(NPTS)                    !Analytic signal
  complex ctmp(NPTS)                    !Analytic signal
  complex cft(512)
  complex cwaveforms(192,64)
  integer, dimension(1) :: iloc
  integer indices(MAXSTEPS)
  integer itone(144)
  logical ismask(NFFT)
  real detmet(-2:MAXSTEPS+3)
  real detfer(MAXSTEPS)
  real ferrs(MAXCAND)
  real pp(12)
  real rcw(12)
  real snrs(MAXCAND)
  real times(MAXCAND)
  real tonespec(NFFT)
  real*8 dt, df, fs, pi, twopi
  logical first
  data first/.true./

  save df,first,cb,cbr,fs,nhashes,pi,twopi,dt,rcw,pp,nmatchedfilter,cwaveforms,rpt

  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
     df=fs/NFFT

     do i=1,12
       angle=(i-1)*pi/12.0
       pp(i)=sin(angle)
       rcw(i)=(1-cos(angle))/2
     enddo

     do i=0,30
       if( i.lt.5 ) then
         write(rpt(i),'(a1,i2.2,a1)') '-',abs(i-5)
         write(rpt(i+31),'(a2,i2.2,a1)') 'R-',abs(i-5)
       else
         write(rpt(i),'(a1,i2.2,a1)') '+',i-5
         write(rpt(i+31),'(a2,i2.2,a1)') 'R+',i-5
       endif
     enddo
     rpt(62)='RRR '
     rpt(63)='73  '

     dphi0=twopi*(freq-500)/12000.0
     dphi1=twopi*(freq+500)/12000.0
     do i=1,64
       msg='<'//trim(mycall)//' '//trim(partnercall)//'> '//rpt(i-1)
       call genmsk32(msg,msgsent,0,itone,itype)
!     write(*,*) i,msg,msgsent,itype
       nsym=32
       phi=0.0
       indx=1
       nreps=1
       do jrep=1,nreps
         do isym=1,nsym
           if( itone(isym) .eq. 0 ) then
             dphi=dphi0
           else
             dphi=dphi1
           endif
           do j=1,6
             cwaveforms(indx,i)=cmplx(cos(phi),sin(phi));
             indx=indx+1
             phi=mod(phi+dphi,twopi)
           enddo
         enddo
       enddo
     enddo

     first=.false.
  endif

! Fill the detmet, detferr arrays
  nstepsize=48  ! 4ms steps
  nstep=(n-NPTS)/nstepsize  
  detmet=0
  detmax=-999.99
  detfer=-999.99
  do istp=1,nstep
    ns=1+nstepsize*(istp-1)
    ne=ns+NPTS-1
    if( ne .gt. n ) exit
    ctmp=cmplx(0.0,0.0)
    ctmp(1:NPTS)=cbig(ns:ne)

! Coarse carrier frequency sync - seek tones at 2000 Hz and 4000 Hz in 
! squared signal spectrum.
! search range for coarse frequency error is +/- 100 Hz

    ctmp=ctmp**2
    ctmp(1:12)=ctmp(1:12)*rcw
    ctmp(NPTS-11:NPTS)=ctmp(NPTS-11:NPTS)*rcw(12:1:-1)
    call four2a(ctmp,NFFT,1,-1,1)
    tonespec=abs(ctmp)**2

    ihlo=(4000-2*ntol)/df+1
    ihhi=(4000+2*ntol)/df+1
    ismask=.false.
    ismask(ihlo:ihhi)=.true.  ! high tone search window
    iloc=maxloc(tonespec,ismask)
    ihpk=iloc(1)
    deltah=-real( (ctmp(ihpk-1)-ctmp(ihpk+1)) / (2*ctmp(ihpk)-ctmp(ihpk-1)-ctmp(ihpk+1)) )
    ah=tonespec(ihpk)
    illo=(2000-2*ntol)/df+1
    ilhi=(2000+2*ntol)/df+1
    ismask=.false.
    ismask(illo:ilhi)=.true.   ! window for low tone
    iloc=maxloc(tonespec,ismask)
    ilpk=iloc(1)
    deltal=-real( (ctmp(ilpk-1)-ctmp(ilpk+1)) / (2*ctmp(ilpk)-ctmp(ilpk-1)-ctmp(ilpk+1)) )
    al=tonespec(ilpk)
    fdiff=(ihpk+deltah-ilpk-deltal)*df
    i2000=2000/df+1
    i4000=4000/df+1
    ferrh=(ihpk+deltah-i4000)*df/2.0
    ferrl=(ilpk+deltal-i2000)*df/2.0
    if( ah .ge. al ) then
      ferr=ferrh
    else
      ferr=ferrl
    endif
    detmet(istp)=max(ah,al)
    detfer(istp)=ferr
!    write(*,*) istp,ilpk,ihpk,ah,al
  enddo  ! end of detection-metric and frequency error estimation loop

  call indexx(detmet(1:nstep),nstep,indices) !find median of detection metric vector
  xmed=detmet(indices(nstep/4))
  detmet=detmet/xmed ! noise floor of detection metric is 1.0
  ndet=0

!do i=1,nstep
!write(77,*) i,detmet(i),detfer(i)
!enddo

  do ip=1,MAXCAND ! find candidates
    iloc=maxloc(detmet(1:nstep))
    il=iloc(1)
    if( (detmet(il) .lt. 4.2) ) exit 
    if( abs(detfer(il)) .le. ntol ) then 
      ndet=ndet+1
      times(ndet)=((il-1)*nstepsize+NPTS/2)*dt
      ferrs(ndet)=detfer(il)
      snrs(ndet)=12.0*log10(detmet(il)-1)/2-8.0
    endif
    detmet(max(1,il-3):min(nstep,il+3))=0.0
!    detmet(il)=0.0
  enddo

  nmessages=0
  lines=char(0)

  fbest=1e6
  pkbest=-1e6
  imsgbest=-1
  istartbest=-1
  ipbest=-1
  nsnrbest=-100
  t0best=-1e6
  do ip=1,ndet  !run through the candidates and try to sync/demod/decode
    imid=times(ip)*fs
    if( imid .lt. NPTS/2 ) imid=NPTS/2
    if( imid .gt. n-NPTS/2 ) imid=n-NPTS/2
    t0=times(ip) + t00
    cdat=cbig(imid-NPTS/2+1:imid+NPTS/2)
    ferr=ferrs(ip)
    nsnr=2*nint(snrs(ip)/2.0)
    if( nsnr .lt. -4 ) nsnr=-4
    if( nsnr .gt. 24 ) nsnr=24

    do imsg=1,64
      do istart=NSPM-NSPM/2,NPTS-NSPM
        cft(1:144)=cdat(istart:istart+144-1)*conjg(cwaveforms(1:144,imsg))
        cft(145:512)=0.
        df=12000.0/512.0
        call four2a(cft,512,1,-1,1)
        iloc=maxloc(abs(cft)) 
        ipk=iloc(1)
        pk=abs(cft(ipk))
        fpk=(ipk-1)*df
        if( fpk.gt.12000.0 ) fpk=fpk-12000.0
        if( pk .gt. pkbest .and. abs(fpk-1500.0) .le. ntol) then
          ipbest=ip
          pkbest=pk
          fbest=fpk
          imsgbest=imsg
          istartbest=istart
          nsnrbest=nsnr
          t0best=t0
        endif
      enddo
    enddo

  enddo  ! candidate loop

999 continue
  msgreceived=' '
  if( imsgbest .gt. 0 .and. pkbest .ge. 108.0) then
           nrxrpt=iand(imsgbest-1,63)
           nrxhash=(imsgbest-1-nrxrpt)/64
!write(*,*) ipbest,pkbest,fbest,imsgbest,istartbest,nsnrbest,t0best,nrxrpt,nrxhash
           nmessages=1
           write(msgreceived,'(a1,a,1x,a,a1,1x,a4)') "<",trim(mycall),      &
                trim(partnercall),">",rpt(nrxrpt)
           write(lines(nmessages),1020) nutc,nsnrbest,t0best,nint(fbest),msgreceived
1020       format(i6.6,i4,f5.1,i5,' & ',a22)
  endif

  return
end subroutine detectmsk32