CODEWITHSUNDEEP
algorithmmatlabfortranfft
Charles
Charles

Reputation: 987

Discrepency in Matlab FFT and Fortran Rosetta Code FFT

I copied the example of resettacode's from

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

and ran the code. As is, results match a simple input to 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

But using a different data set, for example (in 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

The Fortran rosettacode yields

(  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 )

The Rosettacode line were the data initialized was

 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/)

And the do-loop was extended to

  do i=1,10

So my question is: is the data ill conditioned? Or is one of them wrong? Matlab claims to use the same algorithm, so I'm stuck not being sure.

Upvotes: 1

Views: 586

Answers (1)

Alexander Vogt
Alexander Vogt

Reputation: 18098

The Rosetta only provides butterflies for a radix of two, and can, therefore, only handle FFT lengths of 2^N. For a DFT of length ten, you need an additional "butterfly" for the radix five (mixed radix DFT).

You can try perform an FFT of length 16 and to zero-pad the last six entries...

Upvotes: 4

Related Questions