algorithm - Matlab FFT 和 Fortran Rosetta 代码 FFT 的差异

Question

我复制了resettacode的例子

http://rosettacode.org/wiki/Fast_Fourier_transform#Fortran

module fft_mod
  implicit none
  integer,       parameter :: dp=selected_real_kind(15,300)
  real(kind=dp), parameter :: pi=3.141592653589793238460_dp
contains

  ! In place Cooley-Tukey FFT
  recursive subroutine fft(x)
    complex(kind=dp), dimension(:), intent(inout)  :: x
    complex(kind=dp)                               :: t
    integer                                        :: N
    integer                                        :: i
    complex(kind=dp), dimension(:), allocatable    :: even, odd

    N=size(x)

    if(N .le. 1) return

    allocate(odd((N+1)/2))
    allocate(even(N/2))

    ! divide
    odd =x(1:N:2)
    even=x(2:N:2)

    ! conquer
    call fft(odd)
    call fft(even)

    ! combine
    do i=1,N/2
       t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
       x(i)     = odd(i) + t
       x(i+N/2) = odd(i) - t
    end do

    deallocate(odd)
    deallocate(even)

  end subroutine fft

end module fft_mod

program test
  use fft_mod
  implicit none
  complex(kind=dp), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0, 

0.0, 0.0, 0.0/)
  integer :: i

  call fft(data)

  do i=1,8
     write(*,'("(", F20.15, ",", F20.15, "i )")') data(i)
  end do

end program test

然后运行代码。照原样，结果与 matlab R2008a 的简单输入相匹配:

fft([1 1 1 1 0 0 0 0]')

ans =

   4.0000          
   1.0000 - 2.4142i
        0          
   1.0000 - 0.4142i
        0          
   1.0000 + 0.4142i
        0          
   1.0000 + 2.4142i

但是使用不同的数据集，例如(在 Matlab 中):

fft((1:10)')

ans =

  55.0000          
  -5.0000 +15.3884i
  -5.0000 + 6.8819i
  -5.0000 + 3.6327i
  -5.0000 + 1.6246i
  -5.0000          
  -5.0000 - 1.6246i
  -5.0000 - 3.6327i
  -5.0000 - 6.8819i
  -5.0000 -15.3884i

Fortran rosettacode 产量

(  55.000000000000000,   0.000000000000000i )
(   9.854101966249685,   4.081740616606390i )
(   6.854101966249685,  -5.706339097770921i )
(   0.381966011250106,  -9.510565162951536i )
(   0.909830056250527,  -5.877852522924733i )
(  -5.000000000000000,   0.000000000000000i )
(  -2.326237921249264,   3.526711513754838i )
(   3.145898033750315,   5.706339097770921i )
(  12.090169943749473,   1.902113032590309i )
(  17.090169943749473,   5.877852522924733i )

初始化数据的 Rosettacode 行是

 complex(kind=dp), dimension(10) :: data = (/1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0/)

并且 do-loop 被扩展到

  do i=1,10

所以我的问题是:数据是否有问题？或者其中一个是错误的？ Matlab 声称使用相同的算法，所以我不确定。

Answer 1

Rosetta 仅提供基数为 2 的蝴蝶，因此只能处理 2^N 的 FFT 长度。对于长度为 10 的 DFT，您需要一个额外的“蝴蝶”用于基数 5(混合基数 DFT)。

您可以尝试执行长度为 16 的 FFT 并对最后六个条目进行零填充...