WSJT-X/lib/sync4.f90

subroutine sync4(dat,jz,ntol,nfqso,mode,mode4,minwidth,dtx,dfx,snrx,    &
     snrsync,flip,width)

! Synchronizes JT4 data, finding the best-fit DT and DF.  

  parameter (NFFTMAX=2520)         !Max length of FFTs
  parameter (NHMAX=NFFTMAX/2)      !Max length of power spectra
  parameter (NSMAX=525)            !Max number of half-symbol steps
  integer ntol                     !Range of DF search
  real dat(jz)
  real s2(NHMAX,NSMAX)             !2d spectrum, stepped by half-symbols
  real ccfblue(-5:540)             !CCF with pseudorandom sequence
  real ccfred(NHMAX)               !Peak of ccfblue, as function of freq
  real red(NHMAX)                  !Peak of ccfblue, as function of freq
  integer ipk1(1)
  integer nch(7)
  logical savered
  equivalence (ipk1,ipk1a)
  data nch/1,2,4,9,18,36,72/
  save

! Do FFTs of twice symbol length, stepped by half symbols.  Note that 
! we have already downsampled the data by factor of 2.
  nsym=207
  nfft=2520
  nh=nfft/2
  nq=nfft/4
  nsteps=jz/nq - 1
  df=0.5*11025.0/nfft
  ftop=nfqso + 7*mode4*df
  if(ftop.gt.11025.0/4.0) then
     print*,'*** Rx Freq is set too high for this submode ***'
     go to 900
  endif

  if(mode.eq.-999) width=0.                        !Silence compiler warning

  do j=1,nsteps                     !Compute spectrum for each step, get average
     k=(j-1)*nq + 1
     call ps4(dat(k),nfft,s2(1,j))
  enddo

! Set freq and lag ranges
  ia=(nfqso-ntol)/df              !Index of lowest tone, bottom of search range
  ib=(nfqso+ntol)/df              !Index of lowest tone, top of search range
  iamin=nint(100.0/df)
  if(ia.lt.iamin) ia=iamin
  ibmax=nint(2700.0/df) - 6*mode4
  if(ib.gt.ibmax) ib=ibmax

  lag1=-5
  lag2=59
  syncbest=-1.e30
  snrx=-26.0
  ccfred=0.
  red=0.
  i0=nint(nfqso/df)

  do ich=minwidth,7                       !Find best width
     kz=nch(ich)/2
     savered=.false.
     iaa=ia+kz
     ibb=ib-kz
     do i=iaa,ibb                       !Find best frequency channel for CCF
        call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0,   &
             lagpk0,flip)
        ccfred(i)=ccf0
        
! Find rms of the CCF, without main peak
        call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0)
        sync=abs(ccfblue(lagpk0))
!        write(*,3000) ich,i,i*df,ccf0,sync,syncbest
!3000    format(2i5,4f12.3)

! Find best sync value
        if(sync.gt.syncbest*1.03) then
           ipk=i
           lagpk=lagpk0
           ichpk=ich
           syncbest=sync
           savered=.true.
        endif
     enddo
     if(savered) red=ccfred
  enddo
  if(syncbest.lt.-1.e29) go to 900
  ccfred=red
  call pctile(ccfred(ia:ib),ib-ia+1,45,base)
  ccfred=ccfred-base
  
  dfx=ipk*df

! Peak up in time, at best whole-channel frequency
  call xcor4(s2,ipk,nsteps,nsym,lag1,lag2,ichpk,mode4,ccfblue,ccfmax,   &
       lagpk,flip)
  xlag=lagpk
  if(lagpk.gt.lag1 .and. lagpk.lt.lag2) then
     call peakup(ccfblue(lagpk-1),ccfmax,ccfblue(lagpk+1),dx2)
     xlag=lagpk+dx2
  endif

! Find rms of the CCF, without the main peak
  call slope(ccfblue(lag1),lag2-lag1+1,xlag-lag1+1.0)
  sq=0.
  nsq=0
  do lag=lag1,lag2
     if(abs(lag-xlag).gt.2.0) then
        sq=sq+ccfblue(lag)**2
        nsq=nsq+1
     endif
  enddo
  rms=sqrt(sq/nsq)
  snrsync=max(0.0,db(abs(ccfblue(lagpk)/rms - 1.0)) - 4.5)
  dt=2.0/11025.0
  istart=xlag*nq
  dtx=istart*dt

  ipk1=maxloc(ccfred)
  ccf10=0.5*maxval(ccfred)
  do i=ipk1a,ia,-1
     if(ccfred(i).le.ccf10) exit
  enddo
  i1=i
  do i=ipk1a,ib
     if(ccfred(i).le.ccf10) exit
  enddo
  nw=i-i1
  width=nw*df

  sq=0.
  ns=0
  iaa=max(ipk1a-10*nw,ia)
  ibb=min(ipk1a+10*nw,ib)
  jmax=2*mode4/3
  do i=iaa,ibb
     j=abs(i-ipk1a)
     if(j.gt.nw .and. j.lt.jmax) then
        sq=sq + ccfred(j)*ccfred(j)
        ns=ns+1
     endif
  enddo
  rms=sqrt(sq/ns)
  snrx=10.0*log10(ccfred(ipk1a)/rms) - 41.2

900  return
end subroutine sync4