c
c   Shell that makes convolution filter compatible with ANZ
c
      subroutine CheckTcvOpt( CMND,cn,on,lf,hf,order,ftype,logfile,error )
      implicit none
      
      include 'a_struct.inc'

      record /CmndLine/ CMND
                                                         
      integer*4 order,ftype,cn,on                        
      real*4    lf,hf                                              
      logical*4 error,logfile

      error= .true.
      if ( CMND.OptText(on,cn) .eq. 'LOW' ) then
        if ( .not. CMND.OptNumOK(on,cn) ) then
          write(6,100) CMND.OptText(on,cn)
          if ( logfile ) write(3,100) CMND.OptText(on,cn)
100       format(' ERROR- ',a3,' is not a valid low frequency cutoff value'/
     &           '        for the time convolution filter')
        else
          lf= CMND.OptRealVal(on,cn)
          error= .false.
        end if

      else if ( CMND.OptText(on,cn) .eq. 'HIG' ) then
        if ( .not. CMND.OptNumOK(on,cn) ) then
          write(6,110) CMND.OptText(on,cn)
          if ( logfile ) write(3,110) CMND.OptText(on,cn)
110       format(' ERROR- ',a3,' is not a valid high frequency cutoff value'/
     &           '        for the time convolution filter')
        else
          hf= CMND.OptRealVal(on,cn)
          error= .false.
        end if
 
      else if ( CMND.OptText(on,cn) .eq. 'ORD' ) then
        if ( .not. CMND.OptNumOK(on,cn) ) then
          write(6,120) CMND.OptText(on,cn)
          if ( logfile ) write(3,120) CMND.OptText(on,cn)
120       format(' ERROR- ',a3,' is not a valid filter order for the time'/
     &           '        convolution filter. Valid orders are: 4,8,12,16,20')
        else
          order= CMND.OptIntVal(on,cn)
          if ( (mod(order,4) .ne. 0) .or. (order .lt. 4) .or. 
     &         (order .gt. 20) ) then
            order= 1
            write(6,122) order
            if ( logfile ) write(3,122) order
122         format(' ERROR- ',i3,' is not a valid order for the time convolution'/
     &             '         filter: Valid orders are: 4,8,12,16,20')
          else
            order= order/4
            error= .false.                             
          end if
        end if
            
      else if ( CMND.OptText(on,cn) .eq. 'BUT' ) then
        ftype= 0
        error= .false.

      else if ( CMND.OptText(on,cn) .eq. 'CHE' ) then
        ftype= 1
        error= .false.
                                        
      end if
      return
      end                                                                 
c
c   Copyright 1992 by Dwight Meglan
c                               
c   Purpose:
c   Discrete time convolution filter capable of low pass, high pass,
c   bandpass, or notch filtering with variable order of 2,4,6,8,or10
c   and the use of either butterworth or chebychev filter characteristics
c
c   The algorithm used in conv_fltr and conv_lp is based upon a paper by
c   Robert Walraven entitled "Digital Filters" from 1979. 
c                                   
c   Value of conv_fltr upon exiting:    
c		OK : successful filtering
c		ERR_LF1 : low frequency cutoff < 0.0
c		ERR_LF2 : low frequency cutoff > sample frequency
c		ERR_HF1 : high frequency cutoff < 0.0
c		ERR_HF2 : high frequency cutoff > sample frequency
c		ERR_TYP : invalid filter type given
c		ERR_ORD : invalid filter order given
c     ERR_PNT : too many data points for temporary arrays <5000
c
      function conv_fltr( data,np,sf,lf,hf,ftype,order )
      implicit none

      include 'convolve.inc'

      real*4    data(*),sf,lf,hf
      integer*4 np,ftype,order,conv_fltr
                     
      real*4    tmp1(5000),tmp2(5000),f                        
      integer*4 i,filter
c                                       
c     First check over the filter characteristics and make sure that they
c     are logical. Crash gracefully if something is wrong.
c
      if ( lf .lt. 0.0 ) then
        conv_fltr= ERR_CFLTR_LF1
        return                                          
      else if ( lf .gt. sf/2.0 ) then
        conv_fltr= ERR_CFLTR_LF2
        return              
      else if ( hf .lt. 0.0 ) then
        conv_fltr= ERR_CFLTR_HF1
        return    
      else if ( hf .gt. sf/2.0 ) then                               
        conv_fltr= ERR_CFLTR_HF2
        return              
      else if ( (ftype .lt. 0) .and. (ftype .gt. 1) ) then
        conv_fltr= ERR_CFLTR_TYP
        return
      else if ( (order .lt. 1) .and. (order .gt. 10) ) then
        conv_fltr= ERR_CFLTR_ORD
        return
      end if
c
c     Figure out what filtering combination will be carried out
c
      if ( abs(lf) .lt. 1.0e-6 ) then
        f= 2.0*hf/sf
        filter= LP_FLTR
      else if ( abs(hf-sf/2.0) .lt. 1.0e-6 ) then
        f= 2.0*lf/sf
        filter= HP_FLTR
      else
        f= 2.0*lf/sf
        filter= BP_FLTR
        if ( hf .lt. lf ) filter= NTCH_FLTR
      end if
c
c     Now do the filtering
c
      do i= 1, np
        tmp1(i)= data(i)
        tmp2(i)= data(i)
      end do

      call conv_lp( order,f,ftype,data,tmp2,np )

      if ( filter .eq. LP_FLTR ) then
        conv_fltr= 0
        return
      else if ( filter .eq. HP_FLTR ) then
        do i= 1, np
          data(i)= tmp1(i) - data(i)
        end do
        conv_fltr= 0
        return
      end if

      do i= 1, np
        tmp2(i)= tmp1(i) - data(i)
      end do

      call conv_lp( order,2.0*hf/sf,ftype,data,tmp2,np )

      if ( filter .eq. NTCH_FLTR ) then
        do i= 1, np
          data(i)= tmp1(i) - data(i)
        end do
      end if
      conv_fltr= 0
      return
      end
c                    
c
c
      subroutine conv_fltr_err( ierr,logfile )
      implicit none

      include 'convolve.inc'

      integer* 4 ierr
      logical*4 logfile

      if ( ierr .eq. ERR_CFLTR_LF1 ) then
        write(6,100)
        if ( logfile ) write(3,100)
100     format(' Error- Low frequency setting in Convolution filter < 0.0')
                             
      else if ( ierr .eq. ERR_CFLTR_LF2 ) then
        write(6,110)
        if ( logfile ) write(3,110)
110     format(' Error- Low frequency setting in Convolution filter is greater'/
     &         '        than half the sample frequency')
 
      else if ( ierr .eq. ERR_CFLTR_HF1 ) then
        write(6,120)
        if ( logfile ) write(3,120)
120     format(' Error- High frequency setting in Convolution filter < 0.0')

      else if ( ierr .eq. ERR_CFLTR_HF2 ) then                                 
        write(6,130)                                    
        if ( logfile ) write(3,130)
130     format(' Error- High frequency setting in Convolution filter is greater'/
     &         '        than half the sample frequency')

      else if ( ierr .eq. ERR_CFLTR_TYP ) then
        write(6,140)
        if ( logfile ) write(3,140)
140     format(' Error- Invalid filter type for Convolution filter'/
     &         '        Valid filter types: 0 (Butterworth) and 1 (Chebychev)')

      else if ( ierr .eq. ERR_CFLTR_ORD ) then
        write(6,150)
        if ( logfile ) write(3,150)
150     format(' Error- Invalid filter order for Convolution filter'/
     &         '        Orders of 2, 4, 6, 8, and 10 are valid')
      end if
      return
      end
c
c
c
      subroutine conv_lp( order,f,type,din,dout,np )
      implicit none

      integer*4 order,type,np
      real*4    f,din(*),dout(*)

      real*4    c,w,a0,b0,b1,ca,c1,c2,c3,c4,d,hpi,cnst(10),coef(8),ca2
      integer*4 i,j,nt
      
      data cnst / 1.0, 1.3160740,1.5033033,1.6404924,1.7504911,
     &            1.0, 1.2787175,1.4320239,1.5385184,1.6202971 /
      data coef / 0.92103004, 1.3722269, 0.29289321, 0.6776604,
     &            0.0,        1.2900624, 0.17157287, 0.5005791 /
      data hpi / 1.57079633 /

      nt= (order+1)/2
      if ( nt .lt. 1 ) then
        nt= 1
      else if ( nt .gt. 5 ) then
        nt= 5
      end if
      c= cnst(nt+5*type)
      w= coef(type+1)
      a0= coef(type+3)
      b1= coef(type+5)
      b0= coef(type+7)
      d= atan(c*sin(hpi*f)/cos(hpi*f))
      if ( d .gt. hpi ) d= hpi
      ca= sin(w-d)/sin(w+d)
      ca2= ca*ca
      c1= 1.0 + b0*ca2 - b1*ca
      c2= a0*(1.0-ca)*(1.0-ca)/c1
      c3= (b1*(1.0+ca2) - 2.0*ca*(1.0+b0))/c1
      c4=(ca2 + b0 - b1*ca)/c1
c
c     Forward pass- each pass is 2nd order
c
      do j= 1, nt
        dout(1)= din(1)
        dout(2)= din(2)
        do i= 3, np
          dout(i)= c2*(din(i) + 2.0*din(i-1) + din(i-2)) - 
     &             c3*dout(i-1) - c4*dout(i-2)
        end do
        do i= 1, np
          din(i)= dout(i)
        end do
      end do
c
c    Reverse Pass for Phase correction
c 
      do j= 1, nt
        dout(np)= din(np)
        dout(np-1)= din(np-1)
        do i= np-2, 1, -1
          dout(i)= c2*(din(i) + 2.0*din(i+1) + din(i+2)) - 
     &             c3*dout(i+1) - c4*dout(i+2)
        end do
        do i= 1, np
          din(i)= dout(i)
        end do
      end do
      return
      end
