Numeric:FFT inverse_real_fft doesn't trevni! 
Author Message
 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 linux-i386
Copyright 1991-1995 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.65685425-4.j,   4.        +0.j])
Quote:
>>> fra = FFT.real_fft(a,8)
>>> fra

array([ 28.        +0.j,  -9.65685425+4.j,  -4.        -4.j,
              1.65685425-4.j,   4.        +0.j])
Quote:
>>> aa = FFT.inverse_real_fft(fra)
>>> aa

array([-1.93137085+0.j        , -4.35053004-0.53206627j,
             0.24703333-3.08992042j])
Quote:
>>> aa = FFT.inverse_real_fft(fra,8)
>>> aa

array([-1.20710678+0.j        , -2.20710678-2.62132034j,
            -0.79289322+0.79289322j, -0.20710678-1.79289322j,
            -1.62132034+0.j        ])

undocumented-dementedly 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  
 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
Kaiserin-Augusta-Allee 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  
 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  
 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
Kaiserin-Augusta-Allee 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  
 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[N-k] = 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
band-pass 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 real-complex 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  
 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
Kaiserin-Augusta-Allee 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  
 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  
 
 [ 7 post ] 

 Relevant Pages 

1. simple patch so ni doesn't break Numeric

2. Pyrex - Numeric demo doesn't build in 0.7.2

3. tkdesk on fedora core 5 -- doesn't run, doesn't compile

4. BLT installs but doesn't work, and TkTable doesn't build

5. problem with FFT module from Numeric 20.2.0

6. problem with FFT from Numeric 20.2.0

7. Compiler (people Retired) API's ;FFT's;FHT's

8. '-' numeric keyboard in locator field

9. My model doesn't 'walk'

10. 'file mkdir' doesn't fail

11. FFT's in FORTH?

12. FFT asm code in Dr. Dobb's Journal

 

 
Powered by phpBB® Forum Software