      SUBROUTINE FLTRS( FNUM, NPF, WTS )
C**********************************************************************
C
C                         PGM=FLTRS
C
C                         PROGRAMMER =RBLUKAS
C
C                         9/78
C
C**********************************************************************
C
C   FILTERS IS USED TO COMPUTE THE FILTER WEIGHTS FOR DIFFERENT
C   FILTER LENGTHS AND SHAPES OF A SYMMETRIC MOVING AVERAGE FILTER.
C   THE SHAPE OF THE FILTER IS DETERMINED BY THE FILTER #, FNUM.
C
C   FILTER SHAPES ARE: 
C
C                1 - BARTLETT (TRIANGULAR)
C   FNUM =       2 - TUKEY (FULL COSINE / -PI TO PI)
C                3 - PARZEN (HALF COSINE / -PI/2 TO PI/2)
C
C
C   *****NOTE*****
C   THE NUMBER OF POINTS IN THE FILTER, NPTS, MUST BE ODD.
C
C***********************************************************************
      INTEGER FNUM
      DIMENSION WTS(1)
C
C                  CHECK FOR ODD FILTER LENGTH
C
      IF( MOD(NPF,2) .NE. 1 ) THEN
        WRITE(3,*) 'FLTRS - FILTER LENGTH MUST BE ODD'
        STOP
      END IF
      NWTS = (NPF+1)/2
      PI = 4.0*ATAN(1.0)
      SUM = 0.0
      GO TO (100,200,300),FNUM
      WRITE(6,10)FNUM
      STOP
C
C                    BARTLETT FILTER
C
  100 DO 110 I = 2, NWTS
        K = NWTS - I + 1
        WTS(I) = FLOAT(K)
110     SUM = SUM + 2.0*WTS(I)
      WTS(1) = FLOAT(NWTS)
      SUM = SUM + WTS(1)
      DO 120 I = 1, NWTS
120     WTS(I) = WTS(I)/SUM
      RETURN
C
C                 TUKEY FILTER
C
  200 ARG = 2.*PI/FLOAT(NPF)
      DO 210 I = 2, NWTS
        WTS(I) = 1.0 + COS( FLOAT(I-1)*ARG )
210     SUM = SUM + 2. * WTS(I)
      WTS(1) = 2.0
      SUM = SUM + WTS(1)
      DO 220 I = 1, NWTS
220     WTS(I) = WTS(I)/SUM
      RETURN
C
C                 PARZEN FILTER
C
  300 ARG = PI/FLOAT(NPF)
      DO 310 I = 2, NWTS
        WTS(I) = COS( FLOAT(I-1)*ARG )
310     SUM = SUM + 2.*WTS(I) 
      WTS(1) = 1.0
      SUM = SUM + 1.0
      DO 320 I = 1, NWTS
320     WTS(I) = WTS(I)/SUM
      RETURN
   10 FORMAT(1H0,25X,'FNUM = ',I1,' - IS NOT A VALID FILTER NUMBER')
      END 
      SUBROUTINE BPFLTR( FL, FH, NPF, DT, WTS )
C***********************************************************************
C
C                           PGM = BPFLTR
C
C                           PROGRAMMER = R.B.LUKAS
C
C                           7/26/82
C
C***********************************************************************
C
C  BPFLTR IS A BAND PASS FILTER BASED ON THE METHODS OF CHAPTER 6 OF
C  DIGITAL FILTERS BY R. W. HAMMING. FL AND FH ARE THE LOWER AND UPPER
C  CUTOFF FREQUENCIES DESIRED. NPF IS THE NUMBER OF POINTS IN THE FILTER
C  DT IS THE SAMPLING INTERVAL IN THE SAME TIME UNITS AS FL AND FH.
C  THE WEIGHTS ARE RETURNED IN THE VECTOR WTS, WITH THE CENTRAL WEIGHT I
C  THE FIRST ELEMENT, AND THE REST OF THE WEIGHTS ON THE RIGHT HAND SIDE
C  OF THE SYMMETRIC FILTER.
C
C***********************************************************************
      REAL FL, FH, DT
      DIMENSION WTS(1)
      SINK(X) = SIN(X)/X
      NWTS = (NPF+1)/2
      PI = 4.0 * ATAN( 1.00 ) 
      TWOPI = 2.0 * PI
      ARG = PI / FLOAT( NWTS-1 )
      AL = TWOPI*DT*FL
      AH = TWOPI*DT*FH
      WTS(1) = 2.0D0 * ( FH - FL )
      SUM = WTS(1)
C
C                  COMPUTE IDEAL FILTER WEIGHTS WITH
C                  LANCZOS WINDOWING TO REDUCE GIBBS
C                  PHENOMENON 
C
      DO 10 I = 2, NWTS
        FK = FLOAT( I-1 )
        T = SIN( AH*FK ) - SIN( AL*FK ) 
        WTS(I) = SINK(ARG*FK) * T / ( PI*FK*DT )
10      SUM = SUM + 2.0*WTS(I)
C
C                  NORMALIZE
C
      DO 20 K = 1, NWTS
 20     WTS(K) = WTS(K)/SUM
      RETURN
      END 
