home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l292 / 2.ddi / FFT.FOR < prev    next >
Encoding:
Text File  |  1990-02-20  |  9.3 KB  |  381 lines
  1.  
  2.       SUBROUTINE Convolve (filtcoef, x, y, filtLen, aLen)
  3.       INCLUDE 'STDHDR.FOR'
  4.       REAL filtcoef(0:maxv), x(0:maxv), y(0:maxv)
  5.       INTEGER n, m, filtLen, aLen
  6.       REAL summ
  7.  
  8.       DO n = 0, aLen - 1
  9.         summ = 0.0
  10.         IF (n .GE. filtLen - 1) THEN
  11.           DO m = 0, filtLen - 1
  12.         summ = summ + filtcoef(m) * x(n - m)
  13.           END DO
  14.         END IF
  15.         y(n) = summ
  16.       END DO
  17.       END !SUB
  18.  
  19.  
  20.       SUBROUTINE FFTSwap( x,y)
  21.       REAL x,y
  22.       REAL temp
  23.         temp = x
  24.         x = y
  25.         y = temp
  26.       END
  27.  
  28.       REAL FUNCTION ExactBlackman (n, j)
  29.       INCLUDE 'stdhdr.for'
  30.  
  31.       REAL n,j
  32.       ExactBlackman=0.42-0.5*
  33.      +   COS(2.0*pi*j/(n-1))+0.08*COS(4.0*pi*j/(n-1))
  34.       END !FUNCTION
  35.  
  36.  
  37.  
  38.       SUBROUTINE FFT2DCalc (xreal, yimag, c, r, flag)
  39.       INCLUDE 'STDHDR.FOR'
  40.  
  41.       INTEGER c,r,flag
  42.       REAL xreal(0:maxr,0: maxc), yimag(0: maxr, 0: maxc)
  43.  
  44.       INTEGER i, a , iErr
  45.       REAL xVector[ALLOCATABLE](:),yVector[ALLOCATABLE](:)
  46.       ALLOCATE(xVector(0:maxv),yVector(0:maxv),STAT=iErr)
  47.       CALL CheckMem(iErr)
  48.  
  49.        DO j = 0, r - 1
  50.          DO a = 0, c - 1
  51.            xVector(a) = xreal(j, a)
  52.            yVector(a) = yimag(j, a)
  53.          END DO
  54.          IF (flag .EQ. 0) THEN
  55.            CALL FFTCalc(xVector, yVector, c)
  56.          ELSE
  57.            CALL FFTInvCalc(xVector, yVector, c)
  58.          END IF
  59.          DO a = 0, c - 1
  60.            xreal(j, a) = xVector(a)
  61.            yimag(j, a) = yVector(a)
  62.          END DO
  63.        END DO
  64.        DO i = 0, c - 1
  65.          DO a = 0, r - 1
  66.            xVector(a) = xreal(a, i)
  67.            yVector(a) = yimag(a, i)
  68.          END DO
  69.          IF (flag .EQ. 0) THEN
  70.            CALL FFTCalc(xVector, yVector, r)
  71.          ELSE
  72.            CALL FFTInvCalc(xVector, yVector, r)
  73.          END IF
  74.          DO a = 0, r - 1
  75.           xreal(a, i) = xVector(a)
  76.           yimag(a, i) = yVector(a)
  77.          END DO
  78.        END DO
  79.  
  80.       DEALLOCATE(xVector,yVector,STAT=iErr)
  81.       CALL CheckDealloc(iErr)
  82.       END  !SUBROUTINE
  83.  
  84.  
  85.  
  86.       SUBROUTINE FFTCalc (xreal, yimag, numdat)
  87.       INCLUDE 'stdhdr.for'
  88.       INTEGER numdat
  89.       REAL xreal(0:maxv),yimag(0:maxv)
  90.  
  91.       CALL FFTSolve(xreal, yimag, numdat, 0)
  92.  
  93.       END  !SUBROUTINE
  94.  
  95.  
  96.       SUBROUTINE FFTInvCalc (xreal, yimag, numdat)
  97.       INCLUDE 'stdhdr.for'
  98.  
  99.       INTEGER numdat
  100.       REAL xreal(0:maxv),yimag(0:maxv)
  101.  
  102.       CALL FFTSolve(xreal, yimag, numdat, 1)
  103.  
  104.       END  !SUBROUTINE
  105.  
  106.  
  107.  
  108.       SUBROUTINE FFTSolve (xreal, yimag, numdat, flag)
  109.       INCLUDE 'stdhdr.for'
  110.  
  111.       INTEGER numdat,flag
  112.       REAL xreal(0:maxv),yimag(0:maxv)
  113.  
  114.       INTEGER j, i, k, b, iErr
  115.       REAL harm, prodreal, prodimag, sign
  116.       INTEGER maxpower, cntr, arg
  117.       INTEGER pnt0, pnt1
  118.       REAL cosary[ALLOCATABLE](:),sinary[ALLOCATABLE](:)
  119.       ALLOCATE(cosary(0:maxv),sinary(0:maxv),STAT=iErr)
  120.       CALL CheckMem(iErr)
  121.  
  122.       j = 0
  123.       IF (flag .NE. 0) THEN
  124.         sign = 1.0
  125.         DO i = 0, numdat - 1
  126.           xreal(i) = xreal(i) / Real( numdat)
  127.           yimag(i) = yimag(i) / Real( numdat)
  128.         END DO
  129.       ELSE
  130.         sign = -1.0
  131.       END IF
  132.       DO i = 0, numdat - 2
  133.         IF (i .LT. j) THEN
  134.           CALL FFTSwap( xreal(i), xreal(j))
  135.           CALL FFTSwap( yimag(i), yimag(j))
  136.         END IF
  137.         k = numdat / 2
  138.         DO WHILE ( k .LE. j)
  139.           j = j - k
  140.           k = k / 2
  141.         END DO
  142.         j = j + k
  143.       END DO
  144.       maxpower = 0
  145.       i = numdat
  146.       DO WHILE (i .NE. 1)
  147.         maxpower = maxpower + 1
  148.         i = i / 2
  149.       END DO
  150.       harm = 2 * pi / Real( numdat)
  151.       DO i = 0, numdat - 1
  152.         sinary(i) = sign * SIN(harm * i)
  153.        cosary(i) = COS(harm * i)
  154.       END DO
  155.       a = 2
  156.       b = 1
  157.       DO cntr = 1, maxpower
  158.         pnt0 = numdat / a
  159.         pnt1 = 0
  160.         DO k = 0, b - 1
  161.           i = k
  162.           DO WHILE (i .LT. numdat)
  163.             arg = i + b
  164.             IF (k .EQ. 0) THEN
  165.               prodreal = xreal(arg)
  166.               prodimag = yimag(arg)
  167.             ELSE
  168.              prodreal = xreal(arg)*cosary(pnt1)-yimag(arg)*sinary(pnt1)
  169.              prodimag = xreal(arg)*sinary(pnt1)+yimag(arg)*cosary(pnt1)
  170.            END IF
  171.            xreal(arg) = xreal(i) - prodreal
  172.            yimag(arg) = yimag(i) - prodimag
  173.            xreal(i) = xreal(i) + prodreal
  174.            yimag(i) = yimag(i) + prodimag
  175.            i = i + a
  176.           END DO
  177.           pnt1 = pnt1 + pnt0
  178.         END DO
  179.         a = 2 * a
  180.         b = b * 2
  181.       END DO
  182.       DEALLOCATE(cosary,sinary,STAT=iErr)
  183.       CALL CheckDealloc(iErr)
  184.       END  !SUBROUTINE
  185.  
  186.  
  187.  
  188.       REAL FUNCTION Hamming (n, j)
  189.       REAL n,j
  190.       INCLUDE 'STDHDR.FOR'
  191.         Hamming = .54 - .46 * COS(2.0 * pi * j / (n-1))
  192.       END !FUNCTION
  193.  
  194.  
  195.  
  196.       REAL FUNCTION Hanning (n, j)
  197.       INCLUDE 'STDHDR.FOR'
  198.       REAL n,j
  199.         Hanning = .5 * (1 - COS(2.0 * pi * j / (n-1)))
  200.       END !FUNCTION
  201.  
  202.  
  203.  
  204.       REAL FUNCTION Parzen (n, j)
  205.       REAL n,j
  206.         Parzen = 1 - ABS((j - .5 * (n - 1)) / (.5 * (n + 1)))
  207.       END !FUNCTION
  208.  
  209.  
  210.  
  211.       SUBROUTINE PowerSpectrumCalc (xreal, yimag, numdat, delta)
  212.       INCLUDE 'STDHDR.FOR'
  213.  
  214.       INTEGER numdat, i
  215.       REAL xreal(0:maxv),yimag(0:maxv), delta
  216.       REAL numdatr, normal, timespan
  217.  
  218.       numdatr = numdat * 1.0
  219.       normal = 1.0 / (numdatr ** 2)
  220.  
  221.       CALL FFTSolve(xreal, yimag, numdat, 0)
  222.       xreal(0) = (SquareAndSum(xreal(0), yimag(0))) * normal
  223.       DO i = 1, (numdat / 2) - 1
  224.         xreal(i) = normal * (SquareAndSum(xreal(i), yimag(i))
  225.      +           + SquareAndSum(xreal(numdat - i), yimag(numdat - i)))
  226.       END DO
  227.       i = numdat / 2
  228.       xreal(i) = SquareAndSum(xreal(i), yimag(i)) * normal
  229.       timespan = numdat * delta
  230.       DO i = 0, (numdat / 2)
  231.         yimag(i) = REAL(i) / timespan
  232.       END DO
  233.  
  234.       END  !SUBROUTINE
  235.  
  236.  
  237.  
  238.       REAL FUNCTION SquareAndSum (a, b)
  239.       REAL a,b
  240.         SquareAndSum = a ** 2 + b ** 2
  241.       END !FUNCTION
  242.  
  243.  
  244.  
  245.       REAL FUNCTION Welch (n, j)
  246.       REAL n,j
  247.         Welch = 1 - ((j - .5 * (n - 1)) / (.5 * (n + 1))) ** 2
  248.       END !FUNCTION
  249.  
  250.  
  251.  
  252.  
  253.  
  254.  
  255.       SUBROUTINE WindowData (x, numdat, wind)
  256.       INCLUDE 'STDHDR.FOR'
  257.  
  258.       INTEGER numdat,wind
  259.       REAL x(0:maxv)
  260.  
  261.       REAL  multiplier, numdatr, ir
  262.       INTEGER i
  263.  
  264.       numdatr = REAL(numdat)
  265.       DO i = 0, numdat - 1
  266.         ir = REAL(i)
  267.         SELECT CASE (wind)
  268.            CASE (0)
  269.                  multiplier = 1.0
  270.            CASE (1)
  271.                multiplier = Parzen(numdatr, ir)
  272.            CASE (2)
  273.                  multiplier = Hanning(numdatr, ir)
  274.            CASE (3)
  275.                   multiplier = Welch(numdatr, ir)
  276.            CASE (4)
  277.                   multiplier = Hamming(numdatr, ir)
  278.            CASE (5)
  279.                   multiplier = ExactBlackman(numdatr, ir)
  280.            CASE DEFAULT
  281.                   multiplier = 1.0
  282.         END SELECT
  283.         x(i) = multiplier * x(i)
  284.       END DO
  285.  
  286.       END  !SUBROUTINE
  287.  
  288.       SUBROUTINE WindowFFTData (xreal, yimag, numdat, wind)
  289.       INCLUDE 'STDHDR.FOR'
  290.       REAL xreal(0:maxv), yimag(0:maxv)
  291.       INTEGER numdat, wind
  292.  
  293.       CALL WindowData(xreal, numdat, wind)
  294.       CALL WindowData(yimag, numdat, wind)
  295.  
  296.       END !SUBROUTINE
  297.  
  298.  
  299.       SUBROUTINE FIRFreqSample (filtcoef, fp, n, dc, filtType, win)
  300.       INCLUDE 'STDHDR.FOR'
  301.       REAL filtcoef(0:maxv),fp
  302.       INTEGER np, i, j, m, m1, n2, n,dc,filtType,win
  303.       REAL lowbin, highbin, xt, am, q
  304.       REAL a[ALLOCATABLE](:)
  305.       ALLOCATE(a(0:maxv),STAT=iErr)
  306.       Call CheckMem(iErr)
  307.  
  308.       m = (n + 1) / 2
  309.       am = (n + 1) / 2.0
  310.       m1 = n / 2 + 1
  311.       q = 2.0 * pi / n
  312.       n2 = n / 2
  313.       IF (dc .EQ. 1) THEN
  314.          np = INT(n * fp + .5)
  315.       ELSE
  316.           np = INT(n * fp + 1.0)
  317.       END IF
  318.       IF (filtType .EQ. 1) THEN
  319.          lowbin = 0.0
  320.          highbin = 1.0
  321.       ELSE
  322.          lowbin = 1.0
  323.          highbin = 0.0
  324.       END IF
  325.       DO j = 0, np
  326.         a(j) = lowbin
  327.       END DO
  328.  
  329.       DO j = np + 1, m1
  330.         a(j) = highbin
  331.       END DO
  332.  
  333.       IF (dc .EQ. 0) THEN
  334.         DO j = 1, m
  335.       xt = a(1) / 2.0
  336.       DO i = 2, m
  337.         xt = xt + a(i) * COS(q * (am - j) * (i - 1))
  338.       END DO
  339.       filtcoef(j - 1) = 2.0 * xt / (n * 1.0)
  340.       filtcoef(n - j) = filtcoef(j - 1)
  341.         END DO
  342.       ELSE
  343.         DO j = 1, m
  344.       xt = 0.0
  345.       DO i = 1, n2
  346.         xt = xt + a(i) * COS(q * (am - j) * (i - .5))
  347.       END DO
  348.       IF (am .NE. m) xt = xt + a(m) * COS(pi * (am - j)) / 2.0
  349.       filtcoef(j - 1) = 2.0 * xt / n
  350.       filtcoef(n - j) = filtcoef(j - 1)
  351.         END DO
  352.       END IF
  353.       ! Modify filter by Hanning Window if win = 1
  354.       IF (win .GT. 0)  CALL WindowData(filtcoef, n, win)
  355.       DEALLOCATE(a,STAT=iErr)
  356.       Call CheckDealloc(iErr)
  357.       END !SUBROUTINE
  358.  
  359.  
  360.       SUBROUTINE FreqResponse (filtcoef, a, n, k)
  361.       INCLUDE 'STDHDR.FOR'
  362.       REAL filtcoef(0:maxv), a(0:maxv)
  363.       INTEGER i, j, m, n2, n,k
  364.       REAL att, am, q
  365.  
  366.       q = pi / k
  367.       am = (n + 1) / 2.0
  368.       m = (n + 1) / 2
  369.       n2 = n / 2
  370.       DO j = 1, k + 1
  371.         att = 0.0
  372.         IF (am .EQ. m) THEN
  373.           att = filtcoef(m - 1) / 2.0
  374.         END IF
  375.         DO i = 1, n2
  376.          att = att + filtcoef(i - 1) * COS(q * (am - i) * (j - 1))
  377.         END DO
  378.         a(j - 1) = 2.0 * att
  379.       END DO
  380.       END !SUBROUTINE
  381.