WSJT-X/lib/sync4.f90
2018-10-13 19:20:11 -04:00

147 lines
3.9 KiB
Fortran

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(iaa:ibb),ibb-iaa+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