快速傅里叶变换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]