mirror of https://github.com/saitohirga/WSJT-X.git
129 lines
3.7 KiB
Fortran
129 lines
3.7 KiB
Fortran
subroutine filbig(dd,npts,f0,newdat,c4a,n4,sq0)
|
|
|
|
! Filter and downsample the real data in array dd(npts), sampled at 12000 Hz.
|
|
! Output is complex, sampled at 1378.125 Hz.
|
|
|
|
parameter (NSZ=3413)
|
|
parameter (NFFT1=672000,NFFT2=77175)
|
|
parameter (NZ2=1000)
|
|
real*4 dd(npts) !Input data
|
|
real*4 rca(NFFT1)
|
|
complex ca(NFFT1/2+1) !FFT of input
|
|
complex c4a(NFFT2) !Output data
|
|
real*4 s(NZ2)
|
|
real*8 df
|
|
real halfpulse(8) !Impulse response of filter (one sided)
|
|
complex cfilt(NFFT2) !Filter (complex; imag = 0)
|
|
real rfilt(NFFT2) !Filter (real)
|
|
integer*8 plan1,plan2,plan3
|
|
logical first
|
|
include 'fftw3.f90'
|
|
equivalence (rfilt,cfilt),(rca,ca)
|
|
data first/.true./,npatience/1/
|
|
data halfpulse/114.97547150,36.57879257,-20.93789101, &
|
|
5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
|
|
common/refspec/dfref,ref(NSZ)
|
|
save
|
|
|
|
if(npts.lt.0) go to 900 !Clean up at end of program
|
|
|
|
if(first) then
|
|
nflags=FFTW_ESTIMATE
|
|
if(npatience.eq.1) nflags=FFTW_ESTIMATE_PATIENT
|
|
if(npatience.eq.2) nflags=FFTW_MEASURE
|
|
if(npatience.eq.3) nflags=FFTW_PATIENT
|
|
if(npatience.eq.4) nflags=FFTW_EXHAUSTIVE
|
|
! Plan the FFTs just once
|
|
call timer('FFTplans ',0)
|
|
call sfftw_plan_dft_r2c_1d(plan1,nfft1,rca,rca,nflags)
|
|
call sfftw_plan_dft_1d(plan2,nfft2,c4a,c4a,FFTW_FORWARD,nflags)
|
|
call sfftw_plan_dft_1d(plan3,nfft2,cfilt,cfilt,FFTW_BACKWARD,nflags)
|
|
call timer('FFTplans ',1)
|
|
|
|
! Convert impulse response to filter function
|
|
do i=1,nfft2
|
|
cfilt(i)=0.
|
|
enddo
|
|
fac=0.00625/nfft1
|
|
cfilt(1)=fac*halfpulse(1)
|
|
do i=2,8
|
|
cfilt(i)=fac*halfpulse(i)
|
|
cfilt(nfft2+2-i)=fac*halfpulse(i)
|
|
enddo
|
|
call timer('FFTfilt ',0)
|
|
call sfftw_execute(plan3)
|
|
call timer('FFTfilt ',1)
|
|
|
|
base=cfilt(nfft2/2+1)
|
|
do i=1,nfft2
|
|
rfilt(i)=real(cfilt(i))-base
|
|
enddo
|
|
|
|
df=12000.d0/nfft1
|
|
first=.false.
|
|
endif
|
|
|
|
! When new data comes along, we need to compute a new "big FFT"
|
|
! If we just have a new f0, continue with the existing data in ca.
|
|
|
|
if(newdat.ne.0) then
|
|
nz=min(npts,nfft1)
|
|
rca(1:nz)=dd(1:nz)
|
|
rca(nz+1:)=0.
|
|
call timer('FFTbig ',0)
|
|
call sfftw_execute(plan1)
|
|
call timer('FFTbig ',1)
|
|
|
|
do i=1,NFFT1/2 !Flatten the spectrum
|
|
j=nint(i*df/dfref)
|
|
if(j.lt.1) j=1
|
|
if(j.gt.NSZ) j=NSZ
|
|
fac=sqrt(min(30.0,1.0/ref(j)))
|
|
ca(i)=conjg(fac * ca(i))
|
|
enddo
|
|
endif
|
|
|
|
! NB: f0 is the frequency at which we want our filter centered.
|
|
! i0 is the bin number in ca closest to f0.
|
|
|
|
i0=nint(f0/df) + 1
|
|
nh=nfft2/2
|
|
do i=1,nh !Copy data into c4a and apply
|
|
j=i0+i-1 !the filter function
|
|
if(j.ge.1 .and. j.le.nfft1) then
|
|
c4a(i)=rfilt(i)*ca(j)
|
|
else
|
|
c4a(i)=0.
|
|
endif
|
|
enddo
|
|
do i=nh+1,nfft2
|
|
j=i0+i-1-nfft2
|
|
if(j.lt.1) j=j+nfft1 !nfft1 was nfft2
|
|
c4a(i)=rfilt(i)*ca(j)
|
|
enddo
|
|
|
|
nadd=nfft2/NZ2
|
|
i=0
|
|
do j=1,NZ2
|
|
s(j)=0.
|
|
do n=1,nadd
|
|
i=i+1
|
|
s(j)=s(j) + real(c4a(i))**2 + aimag(c4a(i))**2
|
|
enddo
|
|
enddo
|
|
call pctile(s,NZ2,30,sq0)
|
|
|
|
! Do the short reverse transform, to go back to time domain.
|
|
call timer('FFTsmall',0)
|
|
call sfftw_execute(plan2)
|
|
call timer('FFTsmall',1)
|
|
n4=min(npts/8,nfft2)
|
|
return
|
|
|
|
900 call sfftw_destroy_plan(plan1)
|
|
call sfftw_destroy_plan(plan2)
|
|
call sfftw_destroy_plan(plan3)
|
|
|
|
return
|
|
end subroutine filbig
|