WSJT-X/lib/sync4.f90

subroutine sync4(dat,jz,ntol,NFreeze,MouseDF,mode,mode4,minwidth,    &
     dtx,dfx,snrx,snrsync,ccfblue,ccfred1,flip,width,ps0)

! 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 psavg(NHMAX)                !Average spectrum of whole record
  real ps0(450)                    !Avg spectrum for plotting
  real s2(NHMAX,NSMAX)             !2d spectrum, stepped by half-symbols
  real ccfblue(-5:540)             !CCF with pseudorandom sequence
  real ccfred(-450:450)            !Peak of ccfblue, as function of freq
  real red(-450:450)               !Peak of ccfblue, as function of freq
  real ccfred1(-224:224)           !Peak of ccfblue, as function of freq
  real tmp(1260)
  integer ipk1(1)
  integer nch(7)
  logical savered
  equivalence (ipk1,ipk1a)
  data nch/1,2,4,9,18,36,72/
  save

!  write(*,3001) 'A',ntol,nfreeze,mousedf,mode,mode4,minwidth
!3001 format(a1,6i6)

! 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
  psavg(1:nh)=0.
  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))
     psavg(1:nh)=psavg(1:nh) + s2(1:nh,j)
  enddo

  nsmo=min(10*mode4,150)
  call flat1b(psavg,nsmo,s2,nh,nsteps,NHMAX,NSMAX)        !Flatten spectra

  if(mode4.ge.9) call smo(psavg,nh,tmp,mode4/4)

  i0=132
  do i=1,450
     ps0(i)=5.0*(psavg(i0+2*i) + psavg(i0+2*i+1) - 2.0)
  enddo

! Set freq and lag ranges
  famin=200.0 + 3*mode4*df
  fbmax=2700.0 - 3*mode4*df
  fa=famin
  fb=fbmax
  if(NFreeze.eq.1) then
     fa=max(famin,1270.46+MouseDF-ntol)
     fb=min(fbmax,1270.46+MouseDF+ntol)
  else
     fa=max(famin,1270.46+MouseDF-600)
     fb=min(fbmax,1270.46+MouseDF+600)
  endif
  ia=fa/df - 3*mode4                   !Index of lowest tone, bottom of range
  ib=fb/df - 3*mode4                   !Index of lowest tone, top of range
  i0=nint(1270.46/df)
  irange=450
  if(ia-i0.lt.-irange) ia=i0-irange
  if(ib-i0.gt.irange)  ib=i0+irange
  lag1=-5
  lag2=59
  syncbest=-1.e30
  ccfred=0.
  jmax=-1000
  jmin=1000

  do ich=minwidth,7                       !Find best width
     kz=nch(ich)/2
     savered=.false.
     do i=ia+kz,ib-kz                     !Find best frequency channel for CCF
        call xcor4(s2,i,nsteps,nsym,lag1,lag2,ich,mode4,ccfblue,ccf0,   &
             lagpk0,flip)
        j=i-i0 + 3*mode4
        if(j.ge.-372 .and. j.le.372) then
           ccfred(j)=ccf0
           jmax=max(j,jmax)
           jmin=min(j,jmin)
        endif

! Find rms of the CCF, without main peak
        call slope(ccfblue(lag1),lag2-lag1+1,lagpk0-lag1+1.0)
        sync=abs(ccfblue(lagpk0))

! Find best sync value
        if(sync.gt.syncbest) then
           ipk=i
           lagpk=lagpk0
           ichpk=ich
           syncbest=sync
           savered=.true.
        endif
     enddo
     if(savered) red=ccfred
  enddo

  ccfred=red
!  width=df*nch(ichpk)
  dfx=(ipk-i0 + 3*mode4)*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)
  snrx=-26.
  if(mode4.eq.2)  snrx=-25.
  if(mode4.eq.4)  snrx=-24.
  if(mode4.eq.9)  snrx=-23.
  if(mode4.eq.18) snrx=-22.
  if(mode4.eq.36) snrx=-21.
  if(mode4.eq.72) snrx=-20.
  snrx=snrx + snrsync

  dt=2.0/11025.0
  istart=xlag*nq
  dtx=istart*dt
  ccfred1=0.
  jmin=max(jmin,-224)
  jmax=min(jmax,224)
  do i=jmin,jmax
     ccfred1(i)=ccfred(i)
  enddo

  ipk1=maxloc(ccfred1) - 225
  ns=0
  s=0.
  iw=min(mode4,(ib-ia)/4)
  do i=jmin,jmax
     if(abs(i-ipk1a).gt.iw) then
        s=s+ccfred1(i)
        ns=ns+1
     endif
  enddo
  base=s/ns
  ccfred1=ccfred1-base
  ccf10=0.5*maxval(ccfred1)
  do i=ipk1a,jmin,-1
     if(ccfred1(i).le.ccf10) exit
  enddo
  i1=i
  do i=ipk1a,jmax
     if(ccfred1(i).le.ccf10) exit
  enddo
  width=(i-i1)*df

  return
end subroutine sync4

include 'flat1b.f90'