Menu
▾
▴
Re: [Numpy-discussion] real_fft
|
From: <jo...@oh...> - 2002-05-13 09:21:16
|
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: ra...@ph... 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: ban...@so...
> _______________________________________________
> Numpy-discussion mailing list
> Num...@li...
> https://lists.sourceforge.net/lists/listinfo/numpy-discussion
>
>
|
