home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / prgramer / adaptor / examples / f77 / fft.f next >
Text File  |  1993-06-28  |  4KB  |  184 lines
  1.       Program Fast_Fourier_Transform
  2.  
  3. C     Programm zur Berechnung der Potenzen der Einheitswurzeln
  4.  
  5. C     N  Problemgroesse; 2er Potenz
  6.  
  7. C     N2 Anzahl der zu berechnenden Potenzen; N2 = N/2
  8. C     Direkt Fouriersynthese bzw. Fourieranalyse (1/-1)
  9. C     W(N) Feld zur Aufnahme der Potenzen der Einheitswurzeln
  10.  
  11.  
  12.       integer N,N2,k
  13.  
  14.       complex W(:)
  15. cmf$  layout W(:serial)
  16.       complex C(:)
  17.       real F(:), maxf
  18.  
  19.       real Direkt
  20.  
  21.       print * ,'Input k : (N = 2**k)'
  22.       read * ,k
  23.       N = 2**k
  24.       N2 = 2**(k-1)
  25.      
  26.       allocate(W(0:N2-1))
  27.       allocate(C(0:N-1))
  28.       allocate(F(0:N-1))
  29. c
  30. c     random initialization
  31. c
  32.       call cmf_random (f)
  33. c
  34. !HPF$ INDEPENDENT, LOCAL_ACCESS
  35.       do i=0,N-1
  36.         c(i) = f(i)
  37.       end do
  38. c
  39.       Direkt = -1
  40.       call Presort (N,k,C)
  41.       call Roots (W,N2,N,Direkt)
  42.       call Transform(C,W,N,N2,k)
  43. c
  44. c     inverse fft
  45.       Direkt = 1
  46.       call Presort (N,k,C)
  47.       call Roots (W,N2,N,Direkt)
  48.       call Transform(C,W,N,N2,k)
  49. !HPF$ INDEPENDENT, LOCAL_ACCESS
  50.       do i=0,n-1 
  51.         f(i) = abs (c(i)/n - f(i))
  52.       end do
  53.       maxf = 0.0
  54. !HPF$ INDEPENDENT, LOCAL_ACCESS
  55.       do i = 0, N-1
  56.          reduce (maxval, maxf, f(i))
  57.       end do
  58.       print *,'Max = ', maxf
  59.       deallocate(F)
  60.       deallocate(C)
  61.       deallocate(W)
  62.  
  63.       end
  64.  
  65. C     *******************************************************
  66.  
  67. C     -------------------------------------------------------
  68. C     ------------ Umordnung der Speicherplaetze ------------
  69. C     ------------     nach Cooley und Tukey     ------------
  70. C     -------------------------------------------------------
  71.  
  72.       subroutine Presort (N,k,C)
  73.       complex C(0:N-1)
  74.       complex H(0:N-1)
  75.       Integer kl,kn
  76.       integer hi, hj, i, index(0:n-1), ix
  77.  
  78. c     initialization of index
  79.  
  80. !HPF$ INDEPENDENT, LOCAL_ACCESS
  81.       do i = 0, n-1 
  82.          index(i) = i
  83.       end do
  84.  
  85. c     executing bit reverse for all elements of index
  86.  
  87.       do L=1,k
  88.         kl = 2**(L-1)       
  89.         kn = N / (kl + kl)    
  90.  
  91. !HPF$ INDEPENDENT, LOCAL_ACCESS
  92.         do i = 0, n-1 
  93.            ix = index(i)
  94. c          hi = ix / 2
  95.            hi = ishft (ix,-1)
  96. c          hj = hi / kn
  97.            hj = ishft (hi,L-k)
  98. c          hi = mod (hi,kn)
  99.            hi = iand (hi,kn-1)
  100.            index (i) = hi + 2*kn*hj + mod(ix,2)*kn
  101.         end do
  102.        end do
  103.        index = index + 1
  104.        H = C(index-1)
  105. c      call global_get (H, C, index-1)
  106. !HPF$ INDEPENDENT, LOCAL_ACCESS
  107.        do i = 0, n-1 
  108.           C(i) = H(i)
  109.        end do
  110.        end 
  111.  
  112. C     ********************************************************
  113.  
  114.  
  115. C     -------------------------------------------------------
  116. C     ----- Berechnung der Potenzen der Einheitswurzeln -----
  117. C     ----- der Alg. spiegelt eine Rekurrenz erster     -----
  118. C     ----- Ordnung wieder.                             -----
  119. C     -------------------------------------------------------
  120.  
  121.       subroutine Roots (W,N2,N,Direkt)
  122.  
  123.       complex W(0:N2-1),omega
  124. cmf$  layout W(:serial)
  125.       real phi
  126.       real Direkt
  127.  
  128.       phi = Direkt * 8.0 * Atan(1.0) / real(n)
  129.       omega = cmplx(cos(phi),sin(phi))
  130.      
  131.       W(0) = 1
  132.       do i = 1,N2-1
  133.          W(i) = W(i-1) * omega
  134.       end do
  135.       end 
  136.  
  137.  
  138.  
  139. C     ********************************************************
  140.  
  141.  
  142. C    **********************************************************
  143.  
  144.  
  145. C    --------------------------------------------
  146. C    ----- Transformation nach Cooley-Tukey -----
  147. C    --------------------------------------------
  148.  
  149.       subroutine Transform(C,W,N,N2,k)
  150.  
  151. C     Feld c repraesentiert die Fourierkoeffizienten
  152. C     Feld w beinhaltet die Potenzen der Einheitswurzeln
  153.  
  154.       complex C(0:N-1),W(0:N2-1),H(0:N-1)
  155. cmf$  layout W(:serial)
  156.       complex t
  157.  
  158.       integer kl,kn,kl1,hi,index(0:N-1),ni
  159.       
  160.       do L=1,k
  161.         kl = 2**(L-1)
  162.         kn = N2/kl
  163.         kl1 = 2 * kl
  164. !HPF$ INDEPENDENT, LOCAL_ACCESS
  165.         do i = 0, n-1
  166.            index(i) = ieor (i,kl)
  167.         end do
  168.         index = index + 1
  169.         H = C(index-1)
  170. !HPF$ INDEPENDENT, LOCAL_ACCESS
  171.         do i = 0, n-1 
  172. c          hi = mod (i,kl1)
  173.            hi = iand (i,kl1-1)
  174. c          ni = mod (i,kl)
  175.            ni = iand (i,kl-1)
  176.            if (hi .ge. kl) then
  177.               C(i) = H(i) - C(i) * w(ni*kn)
  178.             else
  179.               C(i) = C(i) + w(ni*kn) * H(i)
  180.            end if
  181.         end do
  182.       end do
  183.       end 
  184.