2019-02-07 19:48:38 -05:00
|
|
|
subroutine ft4_downsample(iwave,f0,c)
|
|
|
|
|
|
|
|
|
|
! Input: i*2 data in iwave() at sample rate 12000 Hz
|
|
|
|
|
! Output: Complex data in c(), sampled at 1200 Hz
|
|
|
|
|
|
|
|
|
|
include 'ft4_params.f90'
|
|
|
|
|
parameter (NFFT2=NMAX/16)
|
|
|
|
|
integer*2 iwave(NMAX)
|
|
|
|
|
complex c(0:NMAX/NDOWN-1)
|
|
|
|
|
complex c1(0:NFFT2-1)
|
|
|
|
|
complex cx(0:NMAX/2)
|
|
|
|
|
real x(NMAX), window(0:NFFT2-1)
|
|
|
|
|
equivalence (x,cx)
|
|
|
|
|
logical first
|
|
|
|
|
data first/.true./
|
|
|
|
|
save first,window
|
|
|
|
|
|
|
|
|
|
df=12000.0/NMAX
|
|
|
|
|
baud=12000.0/NSPS
|
|
|
|
|
if(first) then
|
|
|
|
|
bw_transition = 0.5*baud
|
|
|
|
|
bw_flat = 4*baud
|
|
|
|
|
iwt = bw_transition / df
|
|
|
|
|
iwf = bw_flat / df
|
|
|
|
|
pi=4.0*atan(1.0)
|
|
|
|
|
window(0:iwt-1) = 0.5*(1+cos(pi*(/(i,i=iwt-1,0,-1)/)/iwt))
|
|
|
|
|
window(iwt:iwt+iwf-1)=1.0
|
2019-02-14 14:56:02 -05:00
|
|
|
window(iwt+iwf:2*iwt+iwf-1) = 0.5*(1+cos(pi*(/(i,i=0,iwt-1)/)/iwt))
|
2019-02-07 19:48:38 -05:00
|
|
|
window(2*iwt+iwf:)=0.0
|
|
|
|
|
iws = baud / df
|
|
|
|
|
window=cshift(window,iws)
|
|
|
|
|
first=.false.
|
|
|
|
|
endif
|
|
|
|
|
|
|
|
|
|
x=iwave
|
|
|
|
|
call four2a(x,NMAX,1,-1,0) !r2c FFT to freq domain
|
|
|
|
|
i0=nint(f0/df)
|
|
|
|
|
c1=0.
|
|
|
|
|
c1(0)=cx(i0)
|
|
|
|
|
do i=1,NFFT2/2
|
|
|
|
|
if(i0+i.le.NMAX/2) c1(i)=cx(i0+i)
|
|
|
|
|
if(i0-i.ge.0) c1(NFFT2-i)=cx(i0-i)
|
|
|
|
|
enddo
|
|
|
|
|
c1=c1*window/NFFT2
|
|
|
|
|
call four2a(c1,NFFT2,1,1,1) !c2c FFT back to time domain
|
|
|
|
|
c=c1(0:NMAX/NDOWN-1)
|
|
|
|
|
|
|
|
|
|
return
|
|
|
|
|
end subroutine ft4_downsample
|
