[Numpy-discussion] real_fft

[Scott Ransom]

Hi Johan,

> I'm not expecting anybody to look at the whole programs, so I have just cut
> out the important part (however, the complete source is included at the
> bottom of this mail.  The program is creating colored noise)

As a side note, this is a pretty neat function.  Can you give
me a reference for it?  I would like to know exactly what is
going on...

> The problem is (most likely) that the C program uses a library called
> "Numerical Recipes".  In this library there is a function called realft().
> I don't know these FFT functions all that well, but there is a distinct
> difference from the NR (Numerical Recipies) realft() and the FFT.real_fft()
> function of Numerical Python.   This is best illustrated by an example:
> Assume a list/array of 1024 integers.  If you put this array through
> FFT.real_fft() you get a 513 long array as a result.  The NR realft() gives
> a 2048 long array.  Now,  since C cannot deal with complex numbers it has to
> use two entries for each number.  Still the array is much larger than the
> Numpy version.
> 
> Anybody know why ?

Well, the FFT module returns an array that contains only the
positive frequencies (from 0 freq (i.e. the DC value) to the
Nyquist) from the real-valued FFT.  This is N/2+1 values if N
is the number of points in the real-valued FFT.

The Numerical Recipes (NR) routine should return N/2 values
(although you actually get _N_ floats back instead of N/2
complex numbers -- these are the real and imaginary
components).  NR also packs the Nyquist and DC components into
the first two floats (i.e. the first complex number).  This way
you still get all the information, but in the same number of
floats as the original array.  If this is confusing, I
recommend reading the section on how NR packs its FFT arrays.
The book can be found on-line at:

http://www.ulib.org/webRoot/Books/Numerical_Recipes/bookcpdf.html

> wfb and hfb are equal length. Is this a legal way to convolve using
> FFT.real_fft() ?
> 
> wfb = FFT.real_fft(wfb)
> hfb = FFT.real_fft(hfb)
> 
> for i in range(0, len(wfb)):
>        wfb[i] = wfb[i] * hfb[i]
> 
> wfb = FFT.inverse_real_fft(wfb)

Yes.  But it is not very efficient because of the for loop.  I
have modified your code to make it array-based (i.e. using the
wonderful features of Numeric).  Notice that all the for loops
are gone (or at least hidden somewhere underneath in the C
implementation...).  I _think_ it does what you want it to do,
and the only not-quite-so-intuitive thing is the
umath.multiply.accumulate call, which performs th
recurrence-like multiplications in the hfb for-loop.

Cheers,

Scott

---------------

import umath, Numeric, FFT, RandomArray, sys

def noiseGen(points, Qd, b):
    mhb = -b/2.0
    Qd = umath.sqrt(Qd)  # Deviation of the noise
    hfb = Numeric.ones(points, 'd')
    indices = Numeric.arange(len(hfb)-1, typecode='d')
    hfb[1:] = (mhb+indices)/(indices+1.0)
    hfb = umath.multiply.accumulate(hfb)
    wfb = Qd*RandomArray.standard_normal(points)
    return FFT.inverse_real_fft(FFT.real_fft(wfb)*FFT.real_fft(hfb))

if __name__ == '__main__':
    X = noiseGen(2**5, 1, -2)
    for x in X:  print x

---------------

-- 
Scott M. Ransom              Address:  McGill Univ. Physics Dept.
Phone:  (514) 398-6492                 3600 University St., Rm 338
email:  ransom at physics.mcgill.ca       Montreal, QC  Canada H3A 2T8 
GPG Fingerprint: 06A9 9553 78BE 16DB 407B  FFCA 9BFA B6FF FFD3 2989