Numeric:FFT inverse_real_fft doesn't trevni!
Author 
Message 
mhus.. #1 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
The FFT module with the LLNL Numeric stuff does not seem to perform an inverse on the result of real_fft. (At least not in an intuitive manner.) An N point real_fft returns N/2 +1 complex numbers. From experience, not documentation, I guess that element 0 is real and the DC component, and element N/2 is real and the Nyquist value. But, trying inverse_real_fft does not seem to give the number of values or the values I expect. I thought that calling inverse_real_fft on the same array as real_fft produced would give back the original array, possibly scaled by N. Not so. I pasted in a sample of a session. Por favor, how can I invert the real_fft call? % python Python 1.5.1 (#1, Sep 3 1998, 22:51:17) [GCC 2.7.2.3] on linuxi386 Copyright 19911995 Stichting Mathematisch Centrum, Amsterdam Quote: >>> from Numeric import * >>> import FFT >>> a = arange(8) >>> fra = FFT.real_fft(a) >>> fra
array([ 28. +0.j, 9.65685425+4.j, 4. 4.j, 1.656854254.j, 4. +0.j]) Quote: >>> fra = FFT.real_fft(a,8) >>> fra
array([ 28. +0.j, 9.65685425+4.j, 4. 4.j, 1.656854254.j, 4. +0.j]) Quote: >>> aa = FFT.inverse_real_fft(fra) >>> aa
array([1.93137085+0.j , 4.350530040.53206627j, 0.247033333.08992042j]) Quote: >>> aa = FFT.inverse_real_fft(fra,8) >>> aa
array([1.20710678+0.j , 2.207106782.62132034j, 0.79289322+0.79289322j, 0.207106781.79289322j, 1.62132034+0.j ])
undocumenteddementedly yours, Michael Huster
== Sent via Deja.com http://www.***.com/ Share what you know. Learn what you don't.

Mon, 29 Oct 2001 03:00:00 GMT 


Christian Tisme #2 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
Quote:
> The FFT module with the LLNL Numeric stuff does not seem to perform an > inverse on the result of real_fft. (At least not in an intuitive > manner.) An N point real_fft returns N/2 +1 complex numbers. From > experience, not documentation, I guess that element 0 is real and the DC > component, and element N/2 is real and the Nyquist value. > But, trying inverse_real_fft does not seem to give the number of values > or the values I expect. I thought that calling inverse_real_fft on the > same array as real_fft produced would give back the original array, > possibly scaled by N. Not so. I pasted in a sample of a session.
No, that's a misconception. The real fft is just a speedup for cases where your data has only real coefficients. Providing an inverse function would be useless, since after any modifications of your complex result, you will usually not expect a real result. The inverse_real_fft is just a reverse FFT on real data, likewise assuming that you are in the frequency domain with real data. (sounds strange? Exactly as strange as for the time domain, that's what causes all the Nyquist trouble :) ciao  chris 
Applied Biometrics GmbH : Have a break! Take a ride on Python's KaiserinAugustaAllee 101 : *Starship* http://starship.python.net 10553 Berlin : PGP key > http://wwwkeys.pgp.net PGP Fingerprint E182 71C7 1A9D 66E9 9D15 D3CC D4D7 93E2 1FAE F6DF we're tired of banana software  shipped green, ripens at home

Mon, 29 Oct 2001 03:00:00 GMT 


Warren B. Foc #3 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
Quote:
>>I thought that calling inverse_real_fft on the >> same array as real_fft produced would give back the original array, >> possibly scaled by N. Not so. I pasted in a sample of a session. >No, that's a misconception. >The real fft is just a speedup for cases where your data >has only real coefficients. Providing an inverse function >would be useless, since after any modifications of your >complex result, you will usually not expect a real result.
That may be the intention of the FFT module, but it is not the intention of the FFTPACK library on which it is based. I have written working fortran programs using FFTPACK, under the assumption that the rffti function behaves as Michael expected inverse_real_fft to. I find that operation extremely useful for doing things like convolution, interpolation, or differentiation in the Fourier domain. Warren Focke
 "Telopolies running the Internet would be almost exactly like foxes guarding the henhouse, except that foxes are smart and agile." Bob Metcalfe

Mon, 29 Oct 2001 03:00:00 GMT 


Christian Tisme #4 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
Quote:
> >>I thought that calling inverse_real_fft on the > >> same array as real_fft produced would give back the original array, > >> possibly scaled by N. Not so. I pasted in a sample of a session. > >No, that's a misconception. > >The real fft is just a speedup for cases where your data > >has only real coefficients. Providing an inverse function > >would be useless, since after any modifications of your > >complex result, you will usually not expect a real result. > That may be the intention of the FFT module, but it is not the > intention of the FFTPACK library on which it is based. I have written > working FORTRAN programs using FFTPACK, under the assumption that the > rffti function behaves as Michael expected inverse_real_fft to. > I find that operation extremely useful for doing things like > convolution, interpolation, or differentiation in the Fourier domain.
Hmm. So maybe there is a bug in the FFT module? Since you have experience using FFTPACK, it might be helpful if you had a look into FFT.py if this is correct. I'd really be interested to see what's the inverse of what, and couldn't figure it out yet. ciao  chris 
Applied Biometrics GmbH : Have a break! Take a ride on Python's KaiserinAugustaAllee 101 : *Starship* http://starship.python.net 10553 Berlin : PGP key > http://wwwkeys.pgp.net PGP Fingerprint E182 71C7 1A9D 66E9 9D15 D3CC D4D7 93E2 1FAE F6DF we're tired of banana software  shipped green, ripens at home

Tue, 30 Oct 2001 03:00:00 GMT 


Steven G. Johns #5 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
Quote:
> No, that's a misconception. > The real fft is just a speedup for cases where your data > has only real coefficients. Providing an inverse function > would be useless, since after any modifications of your > complex result, you will usually not expect a real result.
That's rather too pessimistic, actually. The output y[] of a real FFT has Hermitian symmetry: y[Nk] = y[k]*. Any modification which preserves this symmetry will continue to produce a real result after an inverse transform. For example, if you multiply the results of two real transforms, the result still has Hermitian symmetry. Thus, the extremely important case of fast real convolutions/correlations benefits from a fast complex>real transform. Most filters, e.g. a bandpass filter, also preserve this symmetry (since you usually want to do the same thing to both positive and negative frequencies). (Another useful case comes in solving differential equations in physical systems with inversion symmetry.) An efficient realcomplex FFT will output only half of the spectrum (the other half being redundant), and the inverse transform will take as input the same half of the spectrum. Thus, the Hermitian symmetry is maintained "automatically" since you can only modify half the spectrum, and the other half is assumed to be symmetric. Cordially, Steven G. Johnson

Wed, 31 Oct 2001 03:00:00 GMT 


Christian Tisme #6 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
... Quote: > >Providing an inverse function > >would be useless, since after any modifications of your > >complex result, you will usually not expect a real result. > Not, true. Since the returned coeficients are only the positive frequencies, > I should be able to do operations such as bandpass, correlate, etc just on > these coefficients then call a *proper* inverse_real_fft and get back a > purely real timeseries.
Yes, sure. Sorry that I forgot about these, after I used FFT for many years (but also many years ago). Quote: > This is the way we have implemented it in C using FFTPACK. FFTPACK has such > a routine which will take the results of real_fft, and invert them. The > wrapped FFTPACK routine in the FFT module just is not the one I've used in > the past.
Well, can you provide the call which you use now? I think this is useful for others. As a remark, the Hartley transform can be used for exactly the operations which you mentioned above, with the same efficiency as real fft, but in a more symmetric manner. No need to have inverse versions at all to mimick the symmetry which exists in the complex transform only. While the DFT is based upon terms like F(v) = (1/N) sum over f(t) exp(J 2pi v t /N) and f(t) = sum over F(v) exp(2 pi v t /N) the discrete Hartley transform (DHT) uses the cas function: H(v) = (1/N) sum over f(t) cas(2 pi v t /N) and f(v) = sum over H(v) cas(2 pi v t /N) where cas is defined as cas(phi) = cos(phi)+sin(phi) DHT can be implemented directly or in terms of DFT. Many properties match each others, and you always just deal with real numbers and minimzed brain convolution. For further reading, see Ronald N. Bracewell, "The Hartley Transform", Oxford Press 1986. Maybe I try to add that to FFT.py in some future. ciao  chris 
Applied Biometrics GmbH : Have a break! Take a ride on Python's KaiserinAugustaAllee 101 : *Starship* http://starship.python.net 10553 Berlin : PGP key > http://wwwkeys.pgp.net PGP Fingerprint E182 71C7 1A9D 66E9 9D15 D3CC D4D7 93E2 1FAE F6DF we're tired of banana software  shipped green, ripens at home

Sat, 03 Nov 2001 03:00:00 GMT 


Michael Huste #7 / 7

Numeric:FFT inverse_real_fft doesn't trevni!
Quote:
> > The FFT module with the LLNL Numeric stuff does not seem to perform an > > inverse on the result of real_fft. (At least not in an intuitive > > manner.) An N point real_fft returns N/2 +1 complex numbers. From > > experience, not documentation, I guess that element 0 is real and the DC > > component, and element N/2 is real and the Nyquist value. > > But, trying inverse_real_fft does not seem to give the number of values > > or the values I expect. I thought that calling inverse_real_fft on the > > same array as real_fft produced would give back the original array, > > possibly scaled by N. Not so. I pasted in a sample of a session. >Christian Tismer replied: >No, that's a misconception. >The real fft is just a speedup for cases where your data >has only real coefficients.
Correct. An N point real series into real_fft returns N/2+1 coefficients. But the first (DC) and last (Nyquist) are purely real, thus the number of independent coefficients are the same. The complex coefficients are the positive frequencies *ONLY*. There IS an algorithm for the inverse with the same speedup. Quote: >Providing an inverse function >would be useless, since after any modifications of your >complex result, you will usually not expect a real result.
Not, true. Since the returned coeficients are only the positive frequencies, I should be able to do operations such as bandpass, correlate, etc just on these coefficients then call a *proper* inverse_real_fft and get back a purely real timeseries. This is the way we have implemented it in C using FFTPACK. FFTPACK has such a routine which will take the results of real_fft, and invert them. The wrapped FFTPACK routine in the FFT module just is not the one I've used in the past. There is the SAME efficiency on the inverse as there is on the forward. We had separate, wrapped calls to set up the work arrays, ezffti_, then a wrapped call to the do the real forward fft, ezfftf_, and finally a wrapped call to ezfftb_ for the inverse. With the wrapping we got the efficiency of the real fft both ways, and the inverse of the forward fft equals (within numerical error) the original timeseries. Quote: >The inverse_real_fft is just a reverse FFT on real data, >likewise assuming that you are in the frequency domain >with real data.
It does not make sense to have purely real frequency domain data. Inverse transforming a purely real frequency domain data results in a minimum time waveform. E.g. a filtered delta function. The analogy between a purely real TIME series (like a pressure) does not go over into purely real frequency domain series. Quote: >(sounds strange? Exactly as strange as for the time domain, >that's what causes all the Nyquist trouble :)
Different issue. _______________________________________________________________ Get Free Email and Do More On The Web. Visit http://www.msn.com

Sat, 03 Nov 2001 03:00:00 GMT 


