WSJT-X/lib/sh65.f90
Joe Taylor 1c0fcc6dec Remove unused code from sh65.f90. Make csync a 2-character flag.
git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@6666 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2016-05-04 18:16:09 +00:00

95 lines
2.1 KiB
Fortran

subroutine sh65(cx,n5,mode65,ntol,xdf,nspecial,snrdb)
parameter(NFFT=2048,NH=NFFT/2,MAXSTEPS=150)
complex cx(90000)
complex c(0:NFFT-1)
real s(-NH+1:NH)
real s2(-NH+1:NH,MAXSTEPS)
real ss(-NH+1:NH,8)
real sigmax(8)
integer ipk(8)
s=0.
ss=0.
jstep=NFFT/4
nblks=n5/jstep - 3
ia=-jstep+1
do iblk=1,nblks
ia=ia+jstep
ib=ia+NFFT-1
c=cx(ia:ib)
call four2a(c,nfft,1,1,1) !c2c FFT
do i=0,NFFT-1
j=i
if(j.gt.NH) j=j-NFFT
p=real(c(i))**2 + aimag(c(i))**2
s(j)=s(j) + p
s2(j,iblk)=p
enddo
n=mod(iblk-1,8) +1
ss(-NH+1:NH,n)=ss(-NH+1:NH,n) + s2(-NH+1:NH,iblk)
enddo
s=1.e-6*s
ss=1.e-6*ss
df=1378.1285/NFFT
do i=-NH+1,NH
f=i*df
write(13,1010) f,s(i),ss(i,1:8)
1010 format(10f10.3)
enddo
nfac=40*mode65
dtstep=0.25/df
! Define range of frequencies to be searched
fa=-ntol
fb=ntol
ia2=max(-NH+1,nint(fa/df))
! Upper tone is above sync tone by 4*nfac*df Hz
ib2=min(NH,nint(fb/df + 4.1*nfac))
! Find strongest line in each of the 4 phases, repeating for each drift rate.
sbest=0.
snrbest=0.
nbest=1
ipk=0
do n=1,8
sigmax(n)=0.
do i=ia2,ib2
sig=ss(i,n)
if(sig.ge.sigmax(n)) then
ipk(n)=i
sigmax(n)=sig
if(sig.ge.sbest) then
sbest=sig
nbest=n
endif
endif
enddo
enddo
n2best=nbest+4
if(n2best.gt.8) n2best=nbest-4
xdf=min(ipk(nbest),ipk(n2best))*df
nspecial=0
if(abs(xdf).le.ntol) then
idiff=abs(ipk(nbest)-ipk(n2best))
xk=float(idiff)/nfac
k=nint(xk)
iderr=nint((xk-k)*nfac)
maxerr=nint(0.008*abs(idiff) + 0.51)
if(abs(iderr).le.maxerr .and. k.ge.2 .and. k.le.4) nspecial=k
snrdb=-30.0
if(nspecial.gt.0) then
call sh65snr(ss(ia2,nbest),ib2-ia2+1,snr1)
call sh65snr(ss(ia2,n2best),ib2-ia2+1,snr2)
snr=0.5*(snr1+snr2)
snrdb=db(snr) - db(2500.0/df) - db(sqrt(nblks/4.0)) + 8.0
endif
if(snr1.lt.4.0 .or. snr2.lt.4.0 .or. snr.lt.5.0) nspecial=0
endif
return
end subroutine sh65