fft from numerical recipes
heiko liebhart
i am quite as new to programming in general and to fortran in particular
as to programming ffts.
as i wanted to apply a fouriertransformation to one of my progs for
pulseshaping i started with the subroutine four1 provided in numerical
recipes.
as a start i used a gaussian as input data expecting to get another
gaussian as output data.
so i generated the input values, filled in the zeros in the imag part in
order to use a real function, called the subroutine an off we go for the
crunching.
after rearranging the outputarray in the way according to the book i had
to realize unfortunately the output i got was not really what i suspected.
what i got was something like a double gaussian. i.e. the sign was
alternating. like sign of every second value was inverted oszillating
from plus to minus.
in the next step i took a sin wave as input.
again what i got was not the suspected delta peak at the frequency of
the oszillation but two antsymmetric peaks in the imaginary part.
according to the symmetry mentioned in the book this should be right for
a real and odd funktion like sin but is not fitting to my
understanding of fouriertransformation.
i was wondering whats going wrong. is the trick in the rearranging of
the output or even deeper in the theory of the transformation so that
maybe i would have to use the product of H and H*.
can somebody please give me a hand to get me through this numerical
jungle even step by step.
i would appreciate every hint and assistance.
thanks alot
Beliavsky
wrote:
> greetings
>
> i am quite as new to programming in general and to fortran in particular
> as to programming ffts.
>
> as i wanted to apply a fouriertransformation to one of my progs for
> pulseshaping i started with the subroutine four1 provided in numerical
> recipes.
The conventions for FFT codes can be tricky. When one purchases
Numerical Recipes, one gets a set of "drivers" showing how the codes
are to be used. For the F90 version of NR, the FFT driver is called
xtwofft.f90 . The F77 version probably has something similar. Have you
run the driver, and do its results make sense to you? If you cannot
understand the results of the driver or if you still have problems
with your own code after you have understood it, the best place to ask
your question is the Numerical Recipes forum at http://www.numerical-recipes.com/forum/
.
me...@skyway.usask.ca
>greetings
>
>i am quite as new to programming in general and to fortran in particular
>as to programming ffts.
>
>as i wanted to apply a fouriertransformation to one of my progs for
>pulseshaping i started with the subroutine four1 provided in numerical
>recipes.
>
>as a start i used a gaussian as input data expecting to get another
>gaussian as output data.
>
>so i generated the input values, filled in the zeros in the imag part in
>order to use a real function, called the subroutine an off we go for the
>crunching.
>
>after rearranging the outputarray in the way according to the book i had
mean splitting the spectrum at the Nyquist to get +- frequency.
>to realize unfortunately the output i got was not really what i suspected.
>
>what i got was something like a double gaussian. i.e. the sign was
>alternating. like sign of every second value was inverted oszillating
>from plus to minus.
>
>in the next step i took a sin wave as input.
>again what i got was not the suspected delta peak at the frequency of
>the oszillation but two antsymmetric peaks in the imaginary part.
>according to the symmetry mentioned in the book this should be right for
> a real and odd funktion like sin but is not fitting to my
>understanding of fouriertransformation.
Do you understand the concept of Nyquist frequency ?
If not, that's a good place to start.
>
>i was wondering whats going wrong. is the trick in the rearranging of
>the output or even deeper in the theory of the transformation so that
>maybe i would have to use the product of H and H*.
>
>can somebody please give me a hand to get me through this numerical
>jungle even step by step.
>
>i would appreciate every hint and assistance.
>
>thanks alot
I could be wrong, but it sounds to me as if you have not had
any previous experience with Fourier transforms, or spectra
If so, then that is the place to start, not FFT. FFT is just
a fast method of calculating a Fourier transform - you
could do it as xr = sum x(t) cos(wt) , xi = sum x(t) sin(wt).
Same result, just not as fast.
Chris
glen herrmannsfeldt
(snip)
> as a start i used a gaussian as input data expecting to get another
> gaussian as output data.
Despite the name FFT, it is actually a Fourier series, that is,
periodic with period equal to the length of the input. The result,
then, is the transform of a periodic series of Gaussians, not
just one. I think that explains the effect you see.
If you make a very long input array with a Gaussian close to
the beginning, that is, approximate the lim period --> infinity,
the result should look more like you expect.
If you want to experiment, it is a little easier with an
interpreted language like R. Even more, R will easily make
graphs so you can see the results, and R is free. After you
find out what you actually want to do, write the Fortran code
for production use.
-- glen
serg271
wrote:
> greetings
>
> i am quite as new to programming in general and to fortran in particular
> as to programming ffts.
>
> as i wanted to apply a fouriertransformation to one of my progs for
> pulseshaping i started with the subroutine four1 provided in numerical
> recipes.
>
> as a start i used a gaussian as input data expecting to get another
> gaussian as output data.
>
> so i generated the input values, filled in the zeros in the imag part in
> order to use a real function, called the subroutine an off we go for the
> crunching.
>
> after rearranging the outputarray in the way according to the book i had
> to realize unfortunately the output i got was not really what i suspected.
>
> what i got was something like a double gaussian. i.e. the sign was
> alternating. like sign of every second value was inverted oszillating
> from plus to minus.
>
Looks like you have some problem with output data rearranging, or
understanding it. Could be that in output even are real parts and odd
are imaginary parts ?
If you are looking for alternative FFT source code, you can also try
KissFFT
http://sourceforge.net/projects/kissfft/
It have implementation for real fft and multidimentional fft, which is
especially relevant to image processing.
Harris
news:f03ar3$kt7$1...@news.albasani.net...
Since many things can cause a FFT routine to produce unexpected results,
from invalid input data to element rearrangings in the output, I would
suggest trying some very typical signals (e.g. delta, pulse, chirp, ...)
with your code and compare it with the results from other software (even
Excel has an option in the Data Analysis pack for Fourier analysis). This is
the only way to be sure that at least your data are as you expected.
--
Harris
Martin Brown
wrote:
I suspect you have not followed the implicit FFT convention that the
array element (1) in Fortran and [0] in C contains the DC component
(and is where your Gaussian should have been centred). It is usually
convenient for display purposes to shift the array along by half its
length.
Try FFT of a delta function in various places in the zeroed array.
That will allow you to get a better handle on the actual convention
used. And if you are only interested in real input move to the real to
complex conjugate symmetric variant of the transform as soon as
possible - no point in doing twice the work and storing redundant
information.
Regards,
Martin Brown
jg.camp...@gmail.com
>
> i am quite as new to programming in general and to fortran in particular
> as to programming ffts.
>
> as i wanted to apply a fouriertransformation to one of my progs for
> pulseshaping i started with the subroutine four1 provided in numerical
> recipes.
>
> as a start i used a gaussian as input data expecting to get another
> gaussian as output data.
>
> so i generated the input values, filled in the zeros in the imag part in
> order to use a real function, called the subroutine an off we go for the
> crunching.
>
> after rearranging the outputarray in the way according to the book i had
> to realize unfortunately the output i got was not really what i suspected.
>
Others have pointed out some pitfalls in the use of FFTs.
It may help to attempt to replicate Exs. 3.6.3--5 in:
http://www.jgcampbell.com/ip/ip3.txt
Your FFT is just one of a number of 'fast' algorithms for computing
the DFT. Your results should be the same as for the DFT; however, a
factor of N, 1/N, 1/sqrt(N), could have been introduced.
Later in that chapter, I show use of 'four1' in C; but since you are
in FORTRAN, that may be little use.
Best regards,
Jon C.
Brendan
properly, what you may be seeing is simply the ringing from the
transform. Ringing is the result of having to work with
discrete(digital) data.
me...@skyway.usask.ca
window sidelobes) . It just makes up the req'd
number of pts that FFT needs - e.g. the simple FFT needs 2^n.
Tapering does reduce sidelobes.
Chris
Jon Harrop
> i am quite as new to programming in general and to fortran in particular
> as to programming ffts.
>
> as i wanted to apply a fouriertransformation to one of my progs for
> pulseshaping i started with the subroutine four1 provided in numerical
> recipes.
Forget about Fortran and NR, they both suck. Use something like Mathematica,
Matlab, R, F# or OCaml and you won't suffer from unnecessary problems. If
you must use an unsafe low-level language use C and the FFTW library from
MIT.
--
Dr Jon D Harrop, Flying Frog Consultancy
The F#.NET Journal
http://www.ffconsultancy.com/products/fsharp_journal/?usenet
Jon Harrop
> Since many things can cause a FFT routine to produce unexpected results,
> from invalid input data to element rearrangings in the output, I would
> suggest trying some very typical signals (e.g. delta, pulse, chirp, ...)
> with your code and compare it with the results from other software (even
> Excel has an option in the Data Analysis pack for Fourier analysis). This
> is the only way to be sure that at least your data are as you expected.
Amazingly, Excel only supports 2^n sizes... :-(
gps
>
> > as i wanted to apply a fouriertransformation to one of my progs for
> > pulseshaping i started with the subroutine four1 provided in numerical
> > recipes.
>
> Forget about Fortran and NR, they both suck. Use something like Mathematica,
> Matlab, R, F# or OCaml and you won't suffer from unnecessary problems. If
> you must use an unsafe low-level language use C and the FFTW library from MIT.
A leaping frog squirt at its finest, what do you think is beneath the
covers of those safeways of yours? But, if you actually knew something
useful you could have suggested any number of places housing FFTs eons
before your newbies shed their diapers, say
Singleton, 1968
http://www.netlib.org/go/fft.f
Steven G. Johnson
wrote:
> what i got was something like a double gaussian. i.e. the sign was
> alternating. like sign of every second value was inverted oszillating
> from plus to minus.
The problem was that you used the wrong origin. The correct origin
for nearly every FFT implementation is the *first* (zeroth) element of
the array. That is, if you want a Gaussian exp(-an^2), your input
array should be x[0] = 1 and x[n] = x[N-n] = exp(-an^2) for 0 < n <= N/
2. Notice the implicit periodic boundaries. (Shift this by 1 for a 1-
based array in Fortran.)
(Even then, the output will be only approximately Gaussian, because
the FFT is not computing the Fourier transform, it is computing an
approximate Fourier series called the discrete Fourier transform.)
What it sounds like you did was the common error of assuming the
origin was at the center of the array, and so your Gaussian was peaked
at N/2. According to the shift theorem, this causes your output X[k]
to be multiplied by a linear phase exp(2 pi i (N/2) k / N) = (-1)^k,
hence the sign oscillation.
Google for "discrete Fourier transform" (DFT) to look at the
transformation an FFT is actually computing, and think about how it
compares to the Fourier transform you are familiar with.
> in the next step i took a sin wave as input.
> again what i got was not the suspected delta peak at the frequency of
> the oszillation but two antsymmetric peaks in the imaginary part.
Um, even for an analytical Fourier transform, the Fourier transform of
sin(kx) is a sum of two imaginary-amplitude delta functions arranged
antisymmetrically around the origin. The same is true in the discrete
case. It sounds like you have a bit of review to do.
Regards,
Steven G. Johnson
Jon Harrop
> A leaping frog squirt at its finest, what do you think is beneath the
> covers of those safeways of yours? But, if you actually knew something
> useful you could have suggested any number of places housing FFTs eons
> before your newbies shed their diapers, say
>
> Singleton, 1968
> http://www.netlib.org/go/fft.f
Is it O(log n) for all "n"? Does it give O(sqrt(log n)) error for all "n"?
No. In fact it will die if you give it an array with a length that has a
prime factor above 23. So your code sucks. Ditch it. Stop citing stuff from
the 60's and pretending it is up to date. Use FFTW. No headache...
heiko liebhart
helllo again
i was overwhelmed by all the feedback and would like to thank you all
for your suggestions and the very usefull posting of the page
http://cimss.ssec.wisc.edu/~paulv/fft/fft_comparison.html
pointing out the problem i am encountering
referring to some suggestion i would like to say.
1)although i am quite new to the topic i am aware of the nyquist theorem
and i am aware of the fact that the procedure is not doing an analytical
FT but a discrete FT i.e. finite fourier series.
2) i also posted in the numerical recipes forum first but got about 80
views and one answer. but nothing helpfull at all.
3)as for the choice of the language fortran was the first choice because
the other program is in fortran
rearanging the data is the deal
thanks for all the assistance
glen herrmannsfeldt
(snip)
> 3)as for the choice of the language fortran was the first choice because
> the other program is in fortran
Yes, that is fine, but Fortran is not the best language for
learning the properties of FFT.
If you use R, you can easily graph the results of the FFT and
quickly learn what the transform of different functions looks
like in the visual sense. It is much easier to see in a graph
than a large table of numbers.
-- glen
Harold Stevens
In <18APR07....@skyway.usask.ca> me...@skyway.usask.ca:
[Snip...]
> number of pts that FFT needs - e.g. the simple FFT needs 2^n.
Piling on here, I'd suggest the OP google "nyquist sampling" too. :)
> Tapering does reduce sidelobes.
<Sigh...> As my 15 or so years tweaking radar simulations can atest, but
IMO getting it "right" consistently is truly a black art. :)
--
Regards, Weird (Harold Stevens) * IMPORTANT EMAIL INFO FOLLOWS *
Pardon any bogus email addresses (wookie) in place for spambots.
Really, it's (wyrd) at airmail, dotted with net. DO NOT SPAM IT.
Kids jumping ship? Looking to hire an old-school type? Email me.