      subroutine rarfilt(npts,nfft,nar,rdata,work,ar,istat)
c
c  birrp
c
c  computes nar term auotoregressive filter (leading 1 suppressed)
c  from robust acvs estimated from median of section averaged power spectrum using
c  tbw 4 data window
c
c  inputs:
c
c  npts  - number of data points in data
c  nfft  - length of section ffts, must be .le. npts
c  nar   - number of filter terms less leading 1, must be .le. nfft/2
c  rdata  - real vector of data
c  work  - real work vector of length at least 18*nfft+14*npts+7*npts/nfft+32
c
c  outputs:
c
c  ar    - vector of nar ar filter terms (leading 1 suppressed)
c  istat - error flag, 0 for normal return, 1 or 2 for error in call to rtpss, 3 or more
c          for ill-conditioned toeplitz matrix 
c
      double precision theta
      dimension rdata(*),ar(*)
      double complex work(*)
      call rtpss(nfft,nfft,1,dble(4./float(nfft)),work,
     $           theta,istat,work(2*nfft+1))
      i1=1
      i2=2*nfft+i1
      i3=8*nfft+30+i2
      i4=7*npts/nfft+2+i3
      i5=4*nfft+i4
      call rarfilt1(npts,nfft,nar,rdata,work(i1),work(i2),
     $              work(i3),work(i4),work(i5),ar,istat)
      return
      end
      subroutine rarfilt1(npts,nfft,nar,rdata,dpss,wsave,
     $                    p,fft,ffts,ar,istat)
      implicit double precision (a-h,o-z)
      parameter (narm=100)
      real rdata,ar
      dimension rdata(*),ar(*)
      double complex dpss(*), wsave(*), p(*)
      double complex fft,ffts,t0,t1
      dimension fft(*),ffts(nfft,*)
      dimension r(narm+1),a(narm)
      noff=0.29*nfft+1
      nsect=float(npts-nfft-1)/float(noff)+1
      if(npts-((nsect-1)*noff+nfft+1).ge.nfft)nsect=nsect+1
      k=1
      nn=nsect/2
      do 5 i=1,nfft
    5   wsave(i)=0.d0
      call dcffti(nfft,wsave)
      do 30 n=1,2*nn,2
        do 10 i=1,nfft
          fft(i)=dpss(i)*dcmplx(rdata(k),rdata(k+noff))
   10     k=k+1
        call dcfftf(nfft,fft,wsave)
        ffts(1,n)=dreal(fft(1))
        ffts(1,n+1)=dimag(fft(1))
        do 20 i=1,nfft/2-1
          t0=fft(i+1)
          t1=fft(nfft-i+1)
          ffts(i+1,n)=0.5*(t0+conjg(t1))
          ffts(i+1,n+1)=-dcmplx(0.,0.5)*(t0-conjg(t1))
   20     continue
        ffts(nfft/2+1,n)=dreal(fft(nfft/2+1))
        ffts(nfft/2+1,n+1)=dimag(fft(nfft/2+1))
   30   k=k-nfft+2*noff
      if(mod(nsect,2).ne.0)then
        do 40 i=1,nfft
          fft(i)=dpss(i)*rdata(k)
   40     k=k+1
        call dcfftf(nfft,fft,wsave)
        do 60 i=1,nfft/2+1
          ffts(i,nsect)=fft(i)
   60     continue
       endif
       do 80 i=1,nfft/2+1
         do 70 j=1,nsect
  70       p(j)=abs(ffts(i,j))**2/nfft
         call dsort(p,p,nsect,1)
         if(mod(nsect,2).eq.0)then
           fft(i)=(p(nsect/2)+p(nsect/2+1))/2.
         else
           fft(i)=p(nsect/2+1)
         endif
  80     continue
      do 90 i=1,nfft/2-1
  90    fft(nfft-i+1)=fft(i+1)
      call dcfftb(nfft,fft,wsave)
      do 110 i=1,nar+1
        r(i)=dreal(fft(i))/nar
 110    continue
      call levinson(nar,r(1),r(2),pp,a,istat)
      if(istat.eq.1)istat=istat+3
      do 120 i=1,nar
 120    ar(i)=a(i)
      return
      end
      subroutine levinson(m,to,t,p,a,istat)
c  solves the set of linear simultaneous equations t*a=p
c  by the levinson algorithm. t is m+1 by m+1 toeplitz, a is a column
c  vector with elements 1,a(1)...a(m), and p is a column vector with p
c  as the top element and zeroes elsewhere.
c
c  input parameters
c
c  m  -- autoregressive order
c  to -- double precision scalar corresponding to matrix element t(0)
c  t  -- double precision vector of acvs lags and left column of toeplitz matrix
c
c  output parameters:
c
c  p  -- double precision scalar of top element of rhs (or ar driving noise)
c  a  -- double precision vector of solution for nar ar coefficients
c  istat--status indicator, 0 for normal exit and 1 for ill-conditioned problem
c
c  converted to real form from S.L. Marple, Digital Spectral Analysis with
c                              Applications,Prentice-Hall, 1987, pp. 104-105
      double precision to,t,p,a,temp,sav
      dimension t(*),a(*)
      p=to
      istat=0
      if(m.eq.0)return
      k=0
  100 k=k+1
      sav=t(k)
      if(k.eq.1)goto 20
      do 10 j=1,k-1
   10   sav=sav+a(j)*t(k-j)
   20 temp=-sav/p
      p=p*(1.d0-temp)*(1.d0+temp)
      if(p.gt.0.)goto 30
      istat=1
      return
   30 a(k)=temp
      if(k.eq.1)goto 50
      khalf=k/2
      do 40 j=1,khalf
        kj=k-j
        sav=a(j)
        a(j)=sav+temp*a(kj)
        if(j.eq.kj)goto 40
        a(kj)=a(kj)+temp*sav
   40   continue
   50 if(k.lt.m)goto 100
      return
      end
      subroutine prewhiten(npts,nar,rdata,ar,pdata)
c
c  birrp
c
c  prewhitens time series in data using ar filter in ar
c  returns npts-nar residuals starting at nar+1 in data
c
c  inputs:
c
c  npts  - number of data points in data, changed to npts-nar on return
c  nar   - number of filter terms less leading 1
c  rdata  - real vector of data
c  ar    - vector of nar ar filter terms (leading 1 is suppressed)
c
c  outputs:
c
c  pdata - real vector of prewhitened data of length npts-nar
c
      dimension rdata(*),ar(*),pdata(*)
      double precision sum
      do 140 i=nar+1,npts
        sum=dble(rdata(i))
        do 130 j=2,nar+1
 130     sum=sum+ar(j-1)*rdata(i-j+1)
 140    pdata(i-nar)=sum
      npts=npts-nar
      return
      end
      subroutine sft(x,n,om,ct,st)
c  calculates fourier transform of real sequence x(i),i=1,...n
c  at angular frequency om normalized so that nyquist=pi. the sine
c  transform is returned in st and the cosine transform in ct.
c  algorithm is that of goertzal with modifications by
c  gentleman, comp.j. 1969
c  transform is not normalized
c  to normalize one-sided ft, divide by sqrt(data length)
c  for positive om, the ft is defined as ct+(0.,1.)st or like slatec
c  cfftf
      implicit double precision (a-h,o-z)
      parameter (pi=3.141592653589793238d0,tp=2.d0*pi)
      dimension x(*)
      np1=n+1
      l=6.d0*om/tp
      s=sin(om)
      a=0.d0
      c=0.d0
      d=0.d0
      e=0.d0
      if(l.eq.0)then
c  recursion for low frequencies (.lt. nyq/3)
        s2=sin(om/2.d0)
        b=-4.d0*s2*s2
        do 10 k=1,n
          c=a
          d=e
          a=x(np1-k)+b*d+c
   10     e=a+d
      elseif(l.eq.1)then
c  regular goertzal algorithm for intermediate frequencies
        b=2.d0*cos(om)
        do 20 k=1,n
          a=x(np1-k)+b*e-d
          d=e
   20     e=a
      else
c  recursion for high frequencies (.gt. 2*nyq/3)
        c2=cos(om/2.d0)
        b=4.d0*c2*c2
        do 30 k=1,n
          c=a
          d=e
          a=x(np1-k)+b*d-c
   30     e=a-d
      endif
      st=-s*d
      ct=a-b*d/2.d0
      return
      end
