
      SUBROUTINE Convolve (filtcoef, x, y, filtLen, aLen)
      INCLUDE 'STDHDR.FOR'
      REAL filtcoef(0:maxv), x(0:maxv), y(0:maxv)
      INTEGER n, m, filtLen, aLen
      REAL summ

      DO n = 0, aLen - 1
        summ = 0.0
        IF (n .GE. filtLen - 1) THEN
          DO m = 0, filtLen - 1
	    summ = summ + filtcoef(m) * x(n - m)
          END DO
        END IF
        y(n) = summ
      END DO
      END !SUB


      SUBROUTINE FFTSwap( x,y)
      REAL x,y
      REAL temp
        temp = x
        x = y
        y = temp
      END

      REAL FUNCTION ExactBlackman (n, j)
      INCLUDE 'stdhdr.for'

      REAL n,j
      ExactBlackman=0.42-0.5*
     +   COS(2.0*pi*j/(n-1))+0.08*COS(4.0*pi*j/(n-1))
      END !FUNCTION



      SUBROUTINE FFT2DCalc (xreal, yimag, c, r, flag)
      INCLUDE 'STDHDR.FOR'

      INTEGER c,r,flag
      REAL xreal(0:maxr,0: maxc), yimag(0: maxr, 0: maxc)

      INTEGER i, a , iErr
      REAL xVector[ALLOCATABLE](:),yVector[ALLOCATABLE](:)
      ALLOCATE(xVector(0:maxv),yVector(0:maxv),STAT=iErr)
      CALL CheckMem(iErr)

       DO j = 0, r - 1
         DO a = 0, c - 1
           xVector(a) = xreal(j, a)
           yVector(a) = yimag(j, a)
         END DO
         IF (flag .EQ. 0) THEN
           CALL FFTCalc(xVector, yVector, c)
         ELSE
           CALL FFTInvCalc(xVector, yVector, c)
         END IF
         DO a = 0, c - 1
           xreal(j, a) = xVector(a)
           yimag(j, a) = yVector(a)
         END DO
       END DO
       DO i = 0, c - 1
         DO a = 0, r - 1
           xVector(a) = xreal(a, i)
           yVector(a) = yimag(a, i)
         END DO
         IF (flag .EQ. 0) THEN
           CALL FFTCalc(xVector, yVector, r)
         ELSE
           CALL FFTInvCalc(xVector, yVector, r)
         END IF
         DO a = 0, r - 1
          xreal(a, i) = xVector(a)
          yimag(a, i) = yVector(a)
         END DO
       END DO

      DEALLOCATE(xVector,yVector,STAT=iErr)
      CALL CheckDealloc(iErr)
      END  !SUBROUTINE



      SUBROUTINE FFTCalc (xreal, yimag, numdat)
      INCLUDE 'stdhdr.for'
      INTEGER numdat
      REAL xreal(0:maxv),yimag(0:maxv)

      CALL FFTSolve(xreal, yimag, numdat, 0)

      END  !SUBROUTINE


      SUBROUTINE FFTInvCalc (xreal, yimag, numdat)
      INCLUDE 'stdhdr.for'

      INTEGER numdat
      REAL xreal(0:maxv),yimag(0:maxv)

      CALL FFTSolve(xreal, yimag, numdat, 1)

      END  !SUBROUTINE



      SUBROUTINE FFTSolve (xreal, yimag, numdat, flag)
      INCLUDE 'stdhdr.for'

      INTEGER numdat,flag
      REAL xreal(0:maxv),yimag(0:maxv)

      INTEGER j, i, k, b, iErr
      REAL harm, prodreal, prodimag, sign
      INTEGER maxpower, cntr, arg
      INTEGER pnt0, pnt1
      REAL cosary[ALLOCATABLE](:),sinary[ALLOCATABLE](:)
      ALLOCATE(cosary(0:maxv),sinary(0:maxv),STAT=iErr)
      CALL CheckMem(iErr)

      j = 0
      IF (flag .NE. 0) THEN
        sign = 1.0
        DO i = 0, numdat - 1
          xreal(i) = xreal(i) / Real( numdat)
          yimag(i) = yimag(i) / Real( numdat)
        END DO
      ELSE
        sign = -1.0
      END IF
      DO i = 0, numdat - 2
        IF (i .LT. j) THEN
          CALL FFTSwap( xreal(i), xreal(j))
          CALL FFTSwap( yimag(i), yimag(j))
        END IF
        k = numdat / 2
        DO WHILE ( k .LE. j)
          j = j - k
          k = k / 2
        END DO
        j = j + k
      END DO
      maxpower = 0
      i = numdat
      DO WHILE (i .NE. 1)
        maxpower = maxpower + 1
        i = i / 2
      END DO
      harm = 2 * pi / Real( numdat)
      DO i = 0, numdat - 1
        sinary(i) = sign * SIN(harm * i)
       cosary(i) = COS(harm * i)
      END DO
      a = 2
      b = 1
      DO cntr = 1, maxpower
        pnt0 = numdat / a
        pnt1 = 0
        DO k = 0, b - 1
          i = k
          DO WHILE (i .LT. numdat)
            arg = i + b
            IF (k .EQ. 0) THEN
              prodreal = xreal(arg)
              prodimag = yimag(arg)
            ELSE
             prodreal = xreal(arg)*cosary(pnt1)-yimag(arg)*sinary(pnt1)
             prodimag = xreal(arg)*sinary(pnt1)+yimag(arg)*cosary(pnt1)
           END IF
           xreal(arg) = xreal(i) - prodreal
           yimag(arg) = yimag(i) - prodimag
           xreal(i) = xreal(i) + prodreal
           yimag(i) = yimag(i) + prodimag
           i = i + a
          END DO
          pnt1 = pnt1 + pnt0
        END DO
        a = 2 * a
        b = b * 2
      END DO
      DEALLOCATE(cosary,sinary,STAT=iErr)
      CALL CheckDealloc(iErr)
      END  !SUBROUTINE



      REAL FUNCTION Hamming (n, j)
      REAL n,j
      INCLUDE 'STDHDR.FOR'
        Hamming = .54 - .46 * COS(2.0 * pi * j / (n-1))
      END !FUNCTION



      REAL FUNCTION Hanning (n, j)
      INCLUDE 'STDHDR.FOR'
      REAL n,j
        Hanning = .5 * (1 - COS(2.0 * pi * j / (n-1)))
      END !FUNCTION



      REAL FUNCTION Parzen (n, j)
      REAL n,j
        Parzen = 1 - ABS((j - .5 * (n - 1)) / (.5 * (n + 1)))
      END !FUNCTION



      SUBROUTINE PowerSpectrumCalc (xreal, yimag, numdat, delta)
      INCLUDE 'STDHDR.FOR'

      INTEGER numdat, i
      REAL xreal(0:maxv),yimag(0:maxv), delta
      REAL numdatr, normal, timespan

      numdatr = numdat * 1.0
      normal = 1.0 / (numdatr ** 2)

      CALL FFTSolve(xreal, yimag, numdat, 0)
      xreal(0) = (SquareAndSum(xreal(0), yimag(0))) * normal
      DO i = 1, (numdat / 2) - 1
        xreal(i) = normal * (SquareAndSum(xreal(i), yimag(i))
     +           + SquareAndSum(xreal(numdat - i), yimag(numdat - i)))
      END DO
      i = numdat / 2
      xreal(i) = SquareAndSum(xreal(i), yimag(i)) * normal
      timespan = numdat * delta
      DO i = 0, (numdat / 2)
        yimag(i) = REAL(i) / timespan
      END DO

      END  !SUBROUTINE



      REAL FUNCTION SquareAndSum (a, b)
      REAL a,b
        SquareAndSum = a ** 2 + b ** 2
      END !FUNCTION



      REAL FUNCTION Welch (n, j)
      REAL n,j
        Welch = 1 - ((j - .5 * (n - 1)) / (.5 * (n + 1))) ** 2
      END !FUNCTION






      SUBROUTINE WindowData (x, numdat, wind)
      INCLUDE 'STDHDR.FOR'

      INTEGER numdat,wind
      REAL x(0:maxv)

      REAL  multiplier, numdatr, ir
      INTEGER i

      numdatr = REAL(numdat)
      DO i = 0, numdat - 1
        ir = REAL(i)
        SELECT CASE (wind)
           CASE (0)
                 multiplier = 1.0
           CASE (1)
               multiplier = Parzen(numdatr, ir)
           CASE (2)
                 multiplier = Hanning(numdatr, ir)
           CASE (3)
                  multiplier = Welch(numdatr, ir)
           CASE (4)
                  multiplier = Hamming(numdatr, ir)
           CASE (5)
                  multiplier = ExactBlackman(numdatr, ir)
           CASE DEFAULT
                  multiplier = 1.0
        END SELECT
        x(i) = multiplier * x(i)
      END DO

      END  !SUBROUTINE

      SUBROUTINE WindowFFTData (xreal, yimag, numdat, wind)
      INCLUDE 'STDHDR.FOR'
      REAL xreal(0:maxv), yimag(0:maxv)
      INTEGER numdat, wind

      CALL WindowData(xreal, numdat, wind)
      CALL WindowData(yimag, numdat, wind)

      END !SUBROUTINE


      SUBROUTINE FIRFreqSample (filtcoef, fp, n, dc, filtType, win)
      INCLUDE 'STDHDR.FOR'
      REAL filtcoef(0:maxv),fp
      INTEGER np, i, j, m, m1, n2, n,dc,filtType,win
      REAL lowbin, highbin, xt, am, q
      REAL a[ALLOCATABLE](:)
      ALLOCATE(a(0:maxv),STAT=iErr)
      Call CheckMem(iErr)

      m = (n + 1) / 2
      am = (n + 1) / 2.0
      m1 = n / 2 + 1
      q = 2.0 * pi / n
      n2 = n / 2
      IF (dc .EQ. 1) THEN
         np = INT(n * fp + .5)
      ELSE
          np = INT(n * fp + 1.0)
      END IF
      IF (filtType .EQ. 1) THEN
         lowbin = 0.0
         highbin = 1.0
      ELSE
         lowbin = 1.0
         highbin = 0.0
      END IF
      DO j = 0, np
        a(j) = lowbin
      END DO

      DO j = np + 1, m1
        a(j) = highbin
      END DO

      IF (dc .EQ. 0) THEN
        DO j = 1, m
	  xt = a(1) / 2.0
	  DO i = 2, m
	    xt = xt + a(i) * COS(q * (am - j) * (i - 1))
	  END DO
	  filtcoef(j - 1) = 2.0 * xt / (n * 1.0)
	  filtcoef(n - j) = filtcoef(j - 1)
        END DO
      ELSE
        DO j = 1, m
	  xt = 0.0
	  DO i = 1, n2
	    xt = xt + a(i) * COS(q * (am - j) * (i - .5))
	  END DO
	  IF (am .NE. m) xt = xt + a(m) * COS(pi * (am - j)) / 2.0
	  filtcoef(j - 1) = 2.0 * xt / n
	  filtcoef(n - j) = filtcoef(j - 1)
        END DO
      END IF
      ! Modify filter by Hanning Window if win = 1
      IF (win .GT. 0)  CALL WindowData(filtcoef, n, win)
      DEALLOCATE(a,STAT=iErr)
      Call CheckDealloc(iErr)
      END !SUBROUTINE


      SUBROUTINE FreqResponse (filtcoef, a, n, k)
      INCLUDE 'STDHDR.FOR'
      REAL filtcoef(0:maxv), a(0:maxv)
      INTEGER i, j, m, n2, n,k
      REAL att, am, q

      q = pi / k
      am = (n + 1) / 2.0
      m = (n + 1) / 2
      n2 = n / 2
      DO j = 1, k + 1
        att = 0.0
        IF (am .EQ. m) THEN
          att = filtcoef(m - 1) / 2.0
        END IF
        DO i = 1, n2
         att = att + filtcoef(i - 1) * COS(q * (am - i) * (j - 1))
        END DO
        a(j - 1) = 2.0 * att
      END DO
      END !SUBROUTINE
