首页 > 代码 > 行业代码 > 正文

代码

阅读排行

快速傅里叶变换FFT
2015-05-23 15:59:21   来源:魔啸天龙   评论:0 点击:

本代码实现 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) = [36.d0,21.d0,33.d0,44.d0,55.d0,63.d0,73.d0,38.d0 ]
  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

相关热词搜索:傅里叶

上一篇:经纬度与WGS84坐标转换
下一篇:mfix-gui两相流计算软件

分享到: 收藏