subroutine filbig(dd,nmax,f0,newdat,nfsample,xpol,c4a,c4b,n4)

C  Filter and downsample complex data stored in array dd(4,nmax).  
C  Output is downsampled from 96000 Hz to 1375.125 Hz.

      parameter (MAXFFT1=5376000,MAXFFT2=77175)
      real*4  dd(4,nmax)                         !Input data
      complex ca(MAXFFT1),cb(MAXFFT1)            !FFTs of input
      complex c4a(MAXFFT2),c4b(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,xpol
      include 'fftw3.f'
      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

      nfft1=MAXFFT1
      nfft2=MAXFFT2
      if(nfsample.eq.95238) then
         nfft1=5120000
         nfft2=74088
      endif
      if(nmax.lt.0) go to 900
      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
C  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(plan2,nfft1,cb,cb,
     +        FFTW_BACKWARD,nflags)
         call sfftw_plan_dft_1d(plan3,nfft2,c4a,c4a,
     +        FFTW_FORWARD,nflags)
         call sfftw_plan_dft_1d(plan4,nfft2,c4b,c4b,
     +        FFTW_FORWARD,nflags)
         call sfftw_plan_dft_1d(plan5,nfft2,cfilt,cfilt,
     +        FFTW_BACKWARD,nflags)
         call timer('FFTplans ',1)

C  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(plan5)
         call timer('FFTfilt ',1)

         base=cfilt(nfft2/2+1)
         do i=1,nfft2
            rfilt(i)=real(cfilt(i))-base
         enddo

         df=96000.d0/nfft1
         if(nfsample.eq.95238) df=95238.1d0/nfft1
         first=.false.
      endif

C  When new data comes along, we need to compute a new "big FFT"
C  If we just have a new f0, continue with the existing ca and cb.

      if(newdat.ne.0) then
         nz=min(nmax,nfft1)
         do i=1,nz
            ca(i)=cmplx(dd(1,i),dd(2,i))
            if(xpol) cb(i)=cmplx(dd(3,i),dd(4,i))
         enddo

         if(nmax.lt.nfft1) then
            do i=nmax+1,nfft1
               ca(i)=0.
               if(xpol) cb(i)=0.
            enddo
         endif
         call timer('FFTbig  ',0)
         call sfftw_execute(plan1)
         if(xpol) call sfftw_execute(plan2)
         call timer('FFTbig  ',1)
         newdat=0
      endif

C  NB: f0 is the frequency at which we want our filter centered.
C      i0 is the bin number in ca and cb closest to f0.

      i0=nint(f0/df) + 1
      nh=nfft2/2
      do i=1,nh                                !Copy data into c4a and c4b,
         j=i0+i-1                              !and apply the filter function
         if(j.ge.1 .and. j.le.nfft1) then
            c4a(i)=rfilt(i)*ca(j)
            if(xpol) c4b(i)=rfilt(i)*cb(j)
         else
            c4a(i)=0.
            if(xpol) c4b(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)
         if(xpol) c4b(i)=rfilt(i)*cb(j)
      enddo

C  Do the short reverse transform, to go back to time domain.
      call timer('FFTsmall',0)
      call sfftw_execute(plan3)
      if(xpol) call sfftw_execute(plan4)
      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