      PROGRAM FILTER

C     Low-pass filter of a time series
C     Time series must be evenly sampled!

      INTEGER INFIL, OUTFIL
      PARAMETER ( INFIL = 10, OUTFIL = 11, MAXNT = 10000,
     &   PI = 3.141592, ALPHA = 0.54 )
      REAL TF( MAXNT ), SF( MAXNT, 51 ), SD( 51 )
      REAL TFIN( MAXNT ), SFIN( MAXNT, 51 )

C     *** Read in time series ***

      WRITE( *, * ) 'FCUT = '
      READ ( 5, * ) FCUT
      READ( INFIL, * ) NSD, ( SD( IS ), IS = 1, NSD )
      NT = 0
      DELTAT = 1.0E6
      DO IT = 1, 10000
         READ( INFIL, *, END = 1500 )
     &      TF( IT ), ( SF( IT, IS ), IS = 1, NSD )
         IF ( IT .GE. 2 )
     &      DELTAT = MIN( TF( IT ) - TF( IT - 1 ), DELTAT )
         NT = NT + 1
      END DO
 1500 DELTAF = 1.0 / ( TF( NT ) - TF( 1 ) )

C     *** Forward transform ***

C     ------ Obtain an N which is a power of 2
      N = ( TF( NT ) - TF( 1 ) ) / DELTAT
      N = 2**( INT( LOG10( REAL( N ) ) / 0.301 ) +1 )

      IF ( N .GT. 4096 ) THEN
         N = 4096
         DELTAT = ( TF( NT ) - TF( 1 ) ) / N
      ENDIF

      LMIN = FCUT / DELTAF
      LMIN = MIN( LMIN, N )
      WRITE( *, * ) LMIN, N

      DO IS = 1, NSD
         CALL INTERP( TF, DELTAT, SF(1,IS), NT, TFIN, SFIN(1,IS), N )

         CALL RFFT(   SFIN( 1, IS ), N, 1 )
C        ------ Scale
         CALL RFFTSC( SFIN( 1, IS ), N, 0, 1 )
C        ------ Unpack
         CALL RFFTSC( SFIN( 1, IS ), N, 3, 0 )

C        Zero out high-freq components

         DO I = 1, N + 1

            IF ( I .LE. LMIN ) THEN
               WEIGHT = ALPHA + ( 1.0 - ALPHA )*COS(PI*(I-1)/(LMIN-1))
            ELSE
               WEIGHT = 0.0
            ENDIF

            SFIN( 2 * I - 1, IS ) = WEIGHT * SFIN( 2 * I - 1, IS )
            SFIN( 2 * I,     IS ) = WEIGHT * SFIN( 2 * I,     IS )
         END DO

C        *** Transform back ***

C        ------ Pack
         CALL RFFTSC( SFIN( 1, IS ), N, -3, 0 )
         CALL RFFT(   SFIN( 1, IS ), N, -1 )
      END DO

C     *** Write out time series ***

      WRITE( OUTFIL, * ) NSD, ( SD( IS ), IS = 1, NSD )
      DO IT = 1, N
         WRITE( OUTFIL, * )
     &      TFIN( IT ), ( SFIN( IT, IS ), IS = 1, NSD )
      END DO

      STOP
      END
C**********************************************************************C
      SUBROUTINE INTERP( WNO, DELWNO, G, NWNO, WNOIN, GIN, NT )

C     Produces an evenly sampled kernel from input G

      REAL WNO( NWNO ), WNOIN( NT ), G( NWNO ), GIN( NT )

C     *** Take steps of DELWNO along the k-axis ***

      II = 1

      DO I = 1, NT
         WNOIN(I) = WNO( 1 ) + ( I - 1 ) * DELWNO

C        ------ Desirable/possible to advance interpolation window?
 5000    IF ( WNOIN( I ) .GT. WNO( II + 1 ) ) THEN
            II = II + 1
            IF ( II .LT. NWNO - 1 ) GO TO 5000
         ENDIF

C        *** Linear interpolation ***

         IF ( WNOIN( I ) .LT. WNO( NWNO ) ) THEN
            T = ( WNOIN( I )    - WNO( II ) )
     &        / ( WNO( II + 1 ) - WNO( II ) )
            GIN( I ) = G( II ) + T * ( G( II + 1 ) - G( II ) )
         ELSE
            GIN( I ) = 0.0
         ENDIF
      END DO

      RETURN
      END
