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.

  use, intrinsic :: iso_c_binding
  use FFTW3

  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
  type(C_PTR) :: plan1,plan2,plan3           !Pointers to FFTW plans
  logical first
!  include 'fftw3.f90'
  equivalence (rfilt,cfilt),(rca,ca)
  data first/.true./
  data halfpulse/114.97547150,36.57879257,-20.93789101,              &
       5.89886379,1.59355187,-2.49138308,0.60910773,-0.04248129/
  common/refspec/dfref,ref(NSZ)
  common/patience/npatience,nthreads
  save first,plan1,plan2,plan3

  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
     plan1=fftwf_plan_dft_r2c_1d(nfft1,rca,ca,nflags)
     plan2=fftwf_plan_dft_1d(nfft2,c4a,c4a,-1,nflags)
     plan3=fftwf_plan_dft_1d(nfft2,cfilt,cfilt,+1,nflags)

! 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 fftwf_execute_dft(plan3,cfilt,cfilt)

     base=real(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
     call timer('FFTbig  ',0)
     nz=min(npts,nfft1)
     rca(1:nz)=dd(1:nz)
     rca(nz+1:)=0.
     call fftwf_execute_dft_r2c(plan1,rca,ca)
     call timer('FFTbig  ',1)

     call timer('flatten ',0)
     ib=0
     do j=1,NSZ
        ia=ib+1
        ib=nint(j*dfref/df)
        fac=sqrt(min(30.0,1.0/ref(j)))
        ca(ia:ib)=fac*conjg(ca(ia:ib))
     enddo
     call timer('flatten ',1)
  endif

! NB: f0 is the frequency at which we want our filter centered.
!     i0 is the bin number in ca closest to f0.
  call timer('loops   ',0)
  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/2+1) 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
     if(j.ge.1) then
        c4a(i)=rfilt(i)*ca(j)
     else
        c4a(i)=rfilt(i)*conjg(ca(2-j))
     endif
  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)
  call timer('loops   ',1)

! Do the short reverse transform, to go back to time domain.
  call timer('FFTsmall',0)
  call fftwf_execute_dft(plan2,c4a,c4a)
  call timer('FFTsmall',1)
  n4=min(npts/8,nfft2)
  return

900 call fftwf_destroy_plan(plan1)
  call fftwf_destroy_plan(plan2)
  call fftwf_destroy_plan(plan3)
  
  return
end subroutine filbig