WIP: cleanup of things related to use of FFTW. More needed!

qmap/libqmap/CMakeLists.txt

@@ -10,10 +10,10 @@ set (libq65_FSRCS
   decode0.f90
   dot.f90
   fchisq0.f90
-  filbig.f90
+  fftbig.f90
-  four2a.f90
+#  four2a.f90
   ftninit.f90
-  ftnquit.f90
+#  ftnquit.f90
   geocentric.f90
   getcand2.f90
   grid2deg.f90

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

qmap/libqmap/filbig.f90

@@ -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

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

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:
 

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 timer('fftbig  ',0)
-  call filbig(dd,NSMAX,f0,newdat,nfsample,cx,n5) !Do the full-length FFT
+  call fftbig(dd,NSMAX) !Do the full-length FFT
-  call timer('filbig  ',1)
+  call timer('fftbig  ',1)
 
   do icand=1,ncand                        !Attempt to decode each candidate
      f0=cand(icand)%f