CODEWITHSUNDEEP
fortranfftfftpack
websealevel
websealevel

Reputation: 91

Trouble with Fourier Transform using fftpack5.1

I have an issue using the FFTPACK5.1 in Fortran 90 which contains subroutines to compute discrete Fourier transforms. I manage to install it and use the routines but when I'm checking if everything is ok with a simple sine wave with a frequency A I get a non zero coefficient not at A (in frequency space, in the spectrum) but at 2A. There is a shift in the spectrum and I don't understand why. I'm almost sure (but I have doubts) that I compute correctly the frequency axis steps:

With N being the number of points of my original sine wave, and Fech my sample frequency I compute the frequency axis steps as df(i)=Fech(i-1)/N.

I'm using the rfft1f routine, so if someone has an experience with it and knows my problem I will be really greatful to understand what is wrong here.

Here is my code:

! n: number of samples in the discret signal
integer ( kind = 4 ), parameter :: n = 4096
real, parameter :: deuxpi=6.283185307
!frequence is the frequence of the signal
!fech is the frequence of sampling
real :: frequence,fech 
integer :: kk
! r is the signal i want to process
! t is the built time and f is the built frequency
real ( kind = 4 ) r(n),t(n),f(n)


!Arrays routines need to work (internal recipe):
real ( kind = 4 ), allocatable, dimension ( : ) :: work
real ( kind = 4 ), allocatable, dimension ( : ) :: wsave

!size of arrays wsave and work for internal recipe 
lensav = n + int ( log ( real ( n, kind = 4 ) ) / log ( 2.0E+00 ) ) + 4
lenwrk = n
allocate ( work(1:lenwrk) )
allocate ( wsave(1:lensav) )


! initializes rttft1f, wsave array
call rfft1i ( n, wsave, lensav, ier )


frequence=0.5
fech=20

! I built here the signal
do kk=1,n
t (kk) = (kk-1) / (fech)
f (kk) = fech*(kk-1) / n
r (kk) = sin( deuxpi * t(kk) * frequence  ) 
end do


!and I call the rfft1f to build the Discrete Fourier Transform:
call rfft1f ( n, inc, r, lenr, wsave, lensav, work, lenwrk, ier )

!I get back r which contains now the fourier coefficients and plot it

I'm expecting with a simple sine wave to have a dirac at the frequency at 0.5 (cf code) but instead I get a dirac at 1., in the frequency domain. Besides, the spectrum looks odd... Here is what I get:

Spectrum

Upvotes: 2

Views: 952

Answers (1)

SleuthEye
SleuthEye

Reputation: 14579

As is typical of routines computing discrete fourier transform of real-valued sequences, the resulting complex valued spectrum is returned for only half the spectrum. To fit the values into the original N-element array, only the real-part of the first value (which is always real) is returned. Similarly in the case of even values of n, the real-part of the n/2-th complex value (which is also always real) is returned. For all other complex values, the real and imaginary parts are interleaved.

So for even n, the corresponding frequency is given by:

r(1)       -> 0
r(2), r(3) -> Delta
r(4), r(5) -> 2*Delta
...
r(n)       -> (n/2)*Delta

And for odd n:

r(1)        -> 0
r(2), r(3)  -> Delta
r(4), r(5)  -> 2*Delta
...
r(n-1),r(n) -> ((n-1)/2) * Delta

where Delta is given as fech/n.

To plot the complex values, you'd probably want to plot either the real (indices 1,2,4,6,...) & imaginary (indices 3,5,7,...) parts as two separate graphs, or the amplitude & phase (again as two separate graphs).

Upvotes: 2

Related Questions