WSJT-X/syncf1.f
J C Dutton a39f662b4d Summary: Merge in linux branch
svn+ssh://svn.berlios.de/svnroot/repos/wsjt/WSJT/branches/linux
merged into
svn+ssh://svn.berlios.de/svnroot/repos/wsjt/trunk



git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/trunk@155 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
2006-04-05 20:07:32 +00:00

118 lines
2.8 KiB
Fortran

subroutine syncf1(data,jz,jstart,f0,NFreeze,DFTolerance,smax,red)
C Does 16k FFTs of data with stepsize 15360, using only "sync on" intervals.
C Returns a refined value of f0, the sync-tone frequency.
parameter (NFFT=16384)
parameter (NH=NFFT/2)
parameter (NQ=NFFT/4)
parameter (NB3=3*512)
real data(jz) !Raw data
integer DFTolerance
real x(NFFT)
real red(512)
real s(NQ) !Ref spectrum for flattening and birdie-zapping
complex c(0:NH)
complex z
equivalence (x,c)
ps(z)=real(z)**2 + aimag(z)**2 !Power spectrum ASF
C Accumulate a high-resolution average spectrum
df=11025.0/NFFT
jstep=10*NB3
nz=(jz-jstart)/jstep -1
call zero(s,NQ)
do n=1,nz
call zero(x,NFFT)
k=(n-1)*jstep
do i=1,10
j=(i-1)*NB3 + 1
call move(data(jstart+k+j),x(j),512)
enddo
call xfft(x,NFFT)
do i=1,NQ
x(i)=ps(c(i))
enddo
call add(s,x,s,NQ)
enddo
fac=(1.0/NFFT)**2
do i=1,NQ !Normalize
s(i)=fac*s(i)
enddo
call smooth(s,NQ)
C NB: could also compute a "blue" spectrum, using the sync-off intervals.
n8=NQ/8
do i=1,n8
red(i)=0.
do k=8*i-7,8*i
red(i)=red(i)+s(k)
enddo
red(i)=10.0*red(i)/(8.0*nz)
enddo
dftol=min(DFTolerance,25)
if(nfreeze.eq.1) dftol=DFTolerance
C Find improved value for f0
smax=0.
ia=(f0-dftol)/df
ib=(f0+dftol)/df + 0.999
! if(NFreeze.eq.1) then
! ia=(f0-5.)/df
! ib=(f0+5.)/df
! endif
do i=ia,ib
if(s(i).gt.smax) then
smax=s(i)
ipk=i
endif
enddo
f0=ipk*df
C Remove line at f0 from spectrum -- if it's strong enough.
ia=(f0-150)/df
ib=(f0+150)/df
a1=0.
a2=0.
nsum=50
do i=1,nsum
a1=a1+s(ia-i)
a2=a2+s(ib+i)
enddo
a1=a1/nsum
a2=a2/nsum
smax=2.0*smax/(a1+a2)
if(smax.gt.3.0) then
b=(a2-a1)/(ib-ia)
do i=ia,ib
s(i)=a1 + (i-ia)*b
enddo
endif
C Make a smoothed version of the spectrum.
nsum=50
fac=1./(2*nsum+1)
call zero(x,nsum)
call zero(s,50)
call zero(s(NQ-nsum),nsum)
sum=0.
do i=nsum+1,NQ-nsum
sum=sum+s(i+nsum)-s(i-nsum)
x(i)=fac*sum
enddo
call zero(x(NQ-nsum),nsum+1)
C To zap birdies, compare s(i) and x(i). If s(i) is larger by more
C than some limit, replace x(i) by s(i). That will put narrow birdies
C on top of the smoothed spectrum.
call move(x,s,NQ) !Copy smoothed spectrum into s
return
end