diff --git a/qmap/libqmap/CMakeLists.txt b/qmap/libqmap/CMakeLists.txt index a1752ba10..5913c6a1e 100644 --- a/qmap/libqmap/CMakeLists.txt +++ b/qmap/libqmap/CMakeLists.txt @@ -10,10 +10,10 @@ set (libq65_FSRCS decode0.f90 dot.f90 fchisq0.f90 - filbig.f90 - four2a.f90 + fftbig.f90 +# four2a.f90 ftninit.f90 - ftnquit.f90 +# ftnquit.f90 geocentric.f90 getcand2.f90 grid2deg.f90 diff --git a/qmap/libqmap/fftbig.f90 b/qmap/libqmap/fftbig.f90 new file mode 100644 index 000000000..383793378 --- /dev/null +++ b/qmap/libqmap/fftbig.f90 @@ -0,0 +1,56 @@ +subroutine fftbig(dd,nmax) + +! Filter and downsample complex data stored in array dd(2,nmax). +! Output is downsampled from 96000 Hz to 1375.125 Hz. + + use timer_module, only: timer + parameter (MAXFFT1=5376000,MAXFFT2=77175) + real*4 dd(2,nmax) !Input data + complex ca(MAXFFT1) !FFT of input + complex c4a(MAXFFT2) !Output data + real*8 df + integer*8 plan1 + logical first + include 'fftw3.f' + common/cacb/ca + equivalence (rfilt,cfilt) + data first/.true./,npatience/1/ + save + + if(nmax.lt.0) go to 900 + + nfft1=MAXFFT1 + 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 big FFT just once + call timer('FFTplan ',0) + call sfftw_plan_dft_1d(plan1,nfft1,ca,ca,FFTW_BACKWARD,nflags) + call timer('FFTplan ',1) + df=96000.d0/nfft1 + first=.false. + endif + + nz=min(nmax,nfft1) + do i=1,nz + ca(i)=cmplx(dd(1,i),dd(2,i)) + enddo + + if(nmax.lt.nfft1) then + do i=nmax+1,nfft1 + ca(i)=0. + enddo + endif + call timer('FFTbig ',0) + call sfftw_execute(plan1) + call timer('FFTbig ',1) + go to 999 + +900 call sfftw_destroy_plan(plan1) + +999 return +end subroutine fftbig diff --git a/qmap/libqmap/filbig.f90 b/qmap/libqmap/filbig.f90 deleted file mode 100644 index a78f1590c..000000000 --- a/qmap/libqmap/filbig.f90 +++ /dev/null @@ -1,117 +0,0 @@ -subroutine filbig(dd,nmax,f0,newdat,nfsample,c4a,n4) - -! Filter and downsample complex data stored in array dd(2,nmax). -! Output is downsampled from 96000 Hz to 1375.125 Hz. - - use timer_module, only: timer - parameter (MAXFFT1=5376000,MAXFFT2=77175) - real*4 dd(2,nmax) !Input data - complex ca(MAXFFT1) !FFT of input - complex c4a(MAXFFT2) !Output data - real*8 df - real halfpulse(8) !Impulse response of filter (one sided) - complex cfilt(MAXFFT2) !Filter (complex; imag = 0) - real rfilt(MAXFFT2) !Filter (real) - integer*8 plan1,plan2,plan3,plan4,plan5 - logical first - include 'fftw3.f' - common/cacb/ca - equivalence (rfilt,cfilt) - data first/.true./,npatience/1/ - data halfpulse/114.97547150,36.57879257,-20.93789101, & - 5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/ - save - - if(nmax.lt.0) go to 900 - - nfft1=MAXFFT1 - nfft2=MAXFFT2 - 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_1d(plan1,nfft1,ca,ca,FFTW_BACKWARD,nflags) - call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a,FFTW_FORWARD,nflags) - call sfftw_plan_dft_1d(plan5,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 sfftw_execute(plan5) - - base=cfilt(nfft2/2+1) - do i=1,nfft2 - rfilt(i)=real(cfilt(i))-base - enddo - - df=96000.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 ca. - - if(newdat.ne.0 .or. sum(abs(ca)).eq.0.0) then !### Test on ca should be unnecessary? - nz=min(nmax,nfft1) - do i=1,nz - ca(i)=cmplx(dd(1,i),dd(2,i)) - enddo - - if(nmax.lt.nfft1) then - do i=nmax+1,nfft1 - ca(i)=0. - enddo - endif - call timer('FFTbig ',0) - call sfftw_execute(plan1) - call timer('FFTbig ',1) -!### newdat=0 - 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 - j=i0+i-1 !and apply 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 - c4a(i)=rfilt(i)*ca(j) - enddo - -! Do the short reverse transform, to go back to time domain. - call timer('FFTsmall',0) - call sfftw_execute(plan3) - call timer('FFTsmall',1) - n4=min(nmax/64,nfft2) - go to 999 - -900 call sfftw_destroy_plan(plan1) - call sfftw_destroy_plan(plan2) - call sfftw_destroy_plan(plan3) - call sfftw_destroy_plan(plan4) - call sfftw_destroy_plan(plan5) - -999 return -end subroutine filbig diff --git a/qmap/libqmap/ftnquit.f90 b/qmap/libqmap/ftnquit.f90 index fe64fac29..e5e98d1e4 100644 --- a/qmap/libqmap/ftnquit.f90 +++ b/qmap/libqmap/ftnquit.f90 @@ -2,7 +2,7 @@ subroutine ftnquit ! Destroy the FFTW plans call four2a(a,-1,1,1,1) - call filbig(id,-1,f0,newdat,nfsample,c4a,n4) + call fftbig(id,-1) return end subroutine ftnquit diff --git a/qmap/libqmap/q65b.f90 b/qmap/libqmap/q65b.f90 index 28c4207c5..5fcd9ba7f 100644 --- a/qmap/libqmap/q65b.f90 +++ b/qmap/libqmap/q65b.f90 @@ -54,7 +54,7 @@ subroutine q65b(nutc,nqd,fcenter,nfcal,nfsample,ikhz,mousedf,ntol, & cx=fac*cx ! Here cx is frequency-domain data around the selected -! QSO frequency, taken from the full-length FFT computed in filbig(). +! QSO frequency, taken from the full-length FFT computed in fftbig(). ! Values for fsample, nfft1, nfft2, df, and the downsampled data rate ! are as follows: diff --git a/qmap/libqmap/qmapa.f90 b/qmap/libqmap/qmapa.f90 index e5322193e..32c49f39a 100644 --- a/qmap/libqmap/qmapa.f90 +++ b/qmap/libqmap/qmapa.f90 @@ -57,9 +57,9 @@ subroutine qmapa(dd,ss,savg,newdat,nutc,fcenter,ntol,nfa,nfb, & bClickDecode=(nagain.eq.1) nagain2=0 - call timer('filbig ',0) - call filbig(dd,NSMAX,f0,newdat,nfsample,cx,n5) !Do the full-length FFT - call timer('filbig ',1) + call timer('fftbig ',0) + call fftbig(dd,NSMAX) !Do the full-length FFT + call timer('fftbig ',1) do icand=1,ncand !Attempt to decode each candidate f0=cand(icand)%f