Use a sin^2 window on FFT in freqcal. Fix the table header for fcal.f90.

git-svn-id: svn+ssh://svn.code.sf.net/p/wsjt/wsjt/branches/wsjtx@7472 ab8295b8-cf94-4d9e-aec4-7959e3be5d79
This commit is contained in:
Joe Taylor 2017-01-10 14:51:32 +00:00
parent 1009d87ab4
commit edef3bc39e
2 changed files with 16 additions and 5 deletions

View File

@ -47,9 +47,9 @@ program fcal
call fit(fd,deltaf,r,iz,a,b,sigmaa,sigmab,rms)
write(*,1002)
1002 format(' Freq DF Meas Freq Resid Call'/ &
1002 format(' Freq DF Meas Freq Resid'/ &
' (MHz) (Hz) (MHz) (Hz)'/ &
'------------------------------------------------')
'-----------------------------------------')
do i=1,iz
fm=fd(i) + 1.d-6*deltaf(i)
calfac=1.d0 + 1.d-6*deltaf(i)/fd(i)

View File

@ -3,18 +3,29 @@ subroutine freqcal(id2,k,nkhz,noffset,ntol,line)
parameter (NZ=30*12000,NFFT=55296,NH=NFFT/2)
integer*2 id2(0:NZ-1)
real x(0:NFFT-1)
real w(0:NFFT-1) !Window function
real s(0:NH)
character line*80,cflag*1
logical first
complex cx(0:NH)
equivalence (x,cx)
data n/0/,k0/9999999/
save n,k0
data n/0/,k0/9999999/,first/.true./
save n,k0,w,first
if(first) then
pi=4.0*atan(1.0)
do i=0,NFFT-1
ww=sin(i*pi/NFFT)
w(i)=ww*ww/NFFT
enddo
first=.false.
endif
if(k.lt.NFFT) go to 900
if(k.lt.k0) n=0
k0=k
x=0.001*id2(k-NFFT:k-1)
x=w*id2(k-NFFT:k-1) !Apply window
call four2a(x,NFFT,1,-1,0) !Compute spectrum, r2c
df=12000.0/NFFT
if (ntol.gt.noffset) then