[Numpy-discussion] real_fft

Johan Fredrik Øhman johanfo at ohman.no
Mon May 13 02:22:02 EDT 2002


First, I would like to say I have solved the problem myself, but thanks to
those who tried to help me :).  The error was related to the way the NR
realft() function interpreted the parameters.



realft(data[],length,1)
I thought length was how long the data array is.  However the routine reads
and uses data
TWICE the length !! :-/



Now, back to the side note.  Yes, colored noise is quite interesting.  It
has applications especially in clock simulations.

I have based my routine on the document "Discrete simulation of power low
noise", IEEE International Frequency Control Symposium 27-29

This is a document not found on the internet, but if you ask your university
or library they might have it.   In the appendix of this document there is a
sample C program which simulates power law noise.


Now that my python port of this is working you can have the code :)
The parameters to the noiseGen are:
points = number of points to generate
X      = the array in which to add the noise¨
Qd    = The noise variance
b       = the power law variable.   f^(b + 2).  Where f is the frequency.
           i.e b = 0 gives what is called white phase noise
                b=-1  gives white flicker noise
                b=-2  white frequency noise (phase random walk noise)
                b=-3  flicker frequency noise (pink noise)
                b=-4  random walk frequency noise
Non-integer values of b is also allowed, .i.e -2.5

############################################################################
##
# MODULE NAME:   Colored Noise Module
# AUTHOR: Johan Fredrik Øhman
#
# This module should simulate the colored noise found in clocks
#
############################################################################
##
import math, FFT, RNG, sys, emath, Numeric

__gausian   = RNG.NormalDistribution(0, 1)
gausian = RNG.CreateGenerator(0, __gausian)

def noiseGen(points, X, Qd, b):
    mhb = -b / 2.0
    Qd = math.sqrt(Qd)  # Deviation of the noise
    hfb = [0] * (points * 2)
    wfb = [0] * (points * 2)
    hfb[0] = 1.0
    wfb[0] = Qd * gausian.ranf()
    for i in range(1, len(wfb)/2):
        # Generate hk coefficients
        hfb[i] = hfb[i-1]/float(i) * (i-1 + mhb)
        # Fille wk with white noise
        wfb[i] = Qd * gausian.ranf()

    hfb = FFT.real_fft(hfb)
    wfb = FFT.real_fft(wfb)

    # Multiply the complex vectors
    # Convolation
    wfb = wfb * hfb
    wfb = FFT.inverse_real_fft(wfb)

    for i in range(0,len(wfb)/2):
        X[i] += wfb[i]


if __name__ == '__main__':
    X = [0] * (2**7)
    noiseGen(2**7, X, 1, -4)
    c = 0
    for i in X:
        print c ,i
        c += 1


> 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
>
> _______________________________________________________________
>
> Have big pipes? SourceForge.net is looking for download mirrors. We supply
> the hardware. You get the recognition. Email Us: bandwidth at sourceforge.net
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>





More information about the NumPy-Discussion mailing list