christy 发表于 2015-11-3 17:05

快速傅里叶变换FFT

本代码实现 Cooley-Tukey 蝶形算法计算复数域的一维快速傅里叶变换。非迭代方式,节约内存,执行速度快。要求数据个数为2的整幂次方,符合语法规范,可直接使用。
如下代码及示范数据,输出为:

(363.000000000000,0.000000000000000E+000)
(-52.9411231676211,-65.4558436870575)
(-15.0000000874228,2.00000000000000)
(14.9411242803314,14.5441577965562)
(31.0000000000000,0.000000000000000E+000)
(14.9411266645322,-14.5441563129425)
(-14.9999999125772,-2.00000000000000)
(-52.9411277772425,65.4558422034438)

(36.0000000000000,0.000000000000000E+000)
(21.0000007843665,-4.589695601353583E-007)
(33.0000000000000,0.000000000000000E+000)
(44.0000003562839,-6.556710155648600E-008)
(55.0000000000000,0.000000000000000E+000)
(62.9999992156335,4.589695601353583E-007)
(73.0000000000000,0.000000000000000E+000)
(37.9999996437161,6.556710155648600E-008)

Module FFT_Mod
Implicit None
Integer , parameter :: DP = Selected_Real_Kind( 15 )
Integer , parameter :: FFT_Forward = 1
Integer , parameter :: FFT_Inverse = -1
contains

Subroutine fcFFT( x , forback )
    !//Subroutine FFT , Cooley-Tukey , radix-2
    !// www.fcode.cn
    Real(Kind=DP) , parameter :: PI = 3.141592654_DP
    Complex(Kind=DP) , Intent(INOUT) :: x(:)
    Integer , Intent(IN) :: forback
    Integer :: n
    integer :: i , j , k , ncur , ntmp , itmp
    real(Kind=DP) :: e
    complex(Kind=DP) :: ctmp
    n = size(x)
    ncur = n
    Do
      ntmp = ncur
      e = 2.0_DP * PI / ncur
      ncur = ncur / 2
      if ( ncur < 1 ) exit
      Do j = 1 , ncur
      Do i = j , n , ntmp
          itmp = i + ncur
          ctmp = x(i) - x(itmp)
          x(i) = x(i) + x(itmp)
          x(itmp) = ctmp * exp( forback * cmplx( 0.0_DP , e*(j-1) ) )
      End Do   
      End Do
    End Do
    j = 1
    Do i = 1, n - 1
      If ( i < j ) then
      ctmp = x(j)
      x(j) = x(i)
      x(i) = ctmp
      End If
      k = n/2
      Do while( k < j )
      j = j - k
      k = k / 2
      End Do
      j = j + k
    End Do
    Return
End Subroutine fcFFT

End Module FFT_Mod

Program www_fcode_cn
use FFT_Mod
Implicit None
integer :: i
integer , parameter :: N = 8
complex(Kind=DP) :: x(N) =
Call fcFFT( x , FFT_Forward )
Do i = 1 , N
    Write (*, *) x(i)
End Do
Call fcFFT( x , FFT_Inverse )
Do i = 1 , N
    Write (*, *) x(i)/N
End Do
End Program www_fcode_cn
页: [1]
查看完整版本: 快速傅里叶变换FFT