WSJT-X/libm65/mskdf.f90

76 lines
1.5 KiB
Fortran

subroutine mskdf(cdat,npts,t2,nfft1,f0,nfreeze,mousedf,ntol,dfx,snrsq2)
! Determine DF for a JTMS signal. Also find ferr, a measure of
! frequency differerence between 1st and 2nd harmonic.
! (Should be 0.000)
parameter (NZ=128*1024)
complex cdat(npts)
integer ntol
real sq(NZ)
real ccf(-2600:2600) !Correct limits?
real tmp(NZ)
complex c(NZ)
data nsps/8/
save c
df1=48000.0/nfft1
nh=nfft1/2
fac=1.0/(nfft1**2)
do i=1,npts
c(i)=fac*cdat(i)**2
enddo
c(npts+1:nfft1)=0.
call four2a(c,nfft1,1,-1,1)
! In the "doubled-frequencies" spectrum of squared cdat:
fa=2.0*(f0-400)
fb=2.0*(f0+400)
j0=nint(2.0*f0/df1)
ja=nint(fa/df1)
jb=nint(fb/df1)
jd=nfft1/nsps
do j=1,nh+1
sq(j)=real(c(j))**2 + aimag(c(j))**2
! if(j*df1.lt.6000.0) write(54,3009) j*df1,sq(j),db(sq(j))
!3009 format(3f12.3)
enddo
ccf=0.
do j=ja,jb
ccf(j-j0-1)=sq(j)+sq(j+jd)
enddo
call pctile(ccf(ja-j0-1),tmp,jb-ja+1,50,base)
ccf=ccf/base
if(NFreeze.gt.0) then
fa=2.0*(f0+MouseDF-ntol)
fb=2.0*(f0+MouseDF+ntol)
endif
ja=nint(fa/df1)
jb=nint(fb/df1)
! rewind 51
smax=0.
do j=ja,jb
k=j-j0-1
if(ccf(k).gt.smax) then
smax=ccf(k)
jpk=j
endif
f=0.5*k*df1
! write(51,3002) f,ccf(k)
!3002 format(2f12.3)
enddo
! call flush(51)
fpk=(jpk-1)*df1
dfx=0.5*fpk-f0
snrsq2=smax
return
end subroutine mskdf