Hi there. I have a problem related to the fft routine in numpy. I'm trying to port a short C program to python, however it has turned out to be a little more complicated than initally thought. 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) 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 ? 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) -- JFØ ==================== PYTHON ==================== import math, FFT, RNG, sys __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 wfb = [0] * points hfb[0] = 1.0 wfb[0] = Qd * gausian.ranf() for i in range(1, len(wfb)): # Generate hk coefficients hfb[i] = hfb[i-1]/float(i) * (i-1 + mhb) # Fille wk with white noise wfb[i] = Qd * gausian.ranf() wfb = FFT.real_fft(wfb) hfb = FFT.real_fft(hfb) print hfb # Multiply the complex vectors # Convolation for i in range(0, len(wfb)): wfb[i] = wfb[i] * hfb[i] wfb = FFT.inverse_real_fft(wfb) for i in range(0,len(wfb)): X[i] += wfb[i] if __name__ == '__main__': X = [0] * (2**10) noiseGen(2**10, X, 1, -2) for i in X: print i ========================== C Program ========================== #include <stdlib.h> #include <stdio.h> #include <string.h> #include <math.h> #include "c/numrec.h" #include "c/NRutil.h" void f_beta(int n_pts, float X[], float Q_d, float b, int *idum){ int i,nn; float *hfb,*wfb; float mhb,wr,wi; nn = n_pts * 2; mhb = -b / 2.0; Q_d = sqrt(Q_d); hfb = vector(1,nn); wfb = vector(1,nn); hfb[1] = 1.0; wfb[1] = Q_d * gasdev(idum); for(i=2; i<=n_pts; i++){ hfb[i] = hfb[i-1]*(mhb+(float)(i-2))/((float)(i-1)); wfb[i] = Q_d * gasdev(idum); } for(i=n_pts;i <= nn; i++){ hfb[i] = 0.0; wfb[i] = 0.0; } realft(hfb, n_pts,1); realft(wfb, n_pts,1); wfb[1]=wfb[1]*hfb[1]; wfb[2]=wfb[2]*hfb[2]; for(i=3; i<= nn; i += 2){ wr = wfb[i]; wi = wfb[i+1]; wfb[i] = wr*hfb[i]-wi*hfb[i+1]; wfb[i+1] = wr*hfb[i+1]+wi*hfb[i]; } realft(wfb, n_pts, -1); for(i=1; i <=n_pts;i++){ X[i] += wfb[i]/((float)n_pts); } free_vector(hfb,1,nn); free_vector(wfb,1,nn); } int main(){ int length22; int i; int idum = -210310212; float *X; X = vector(1,1024); f_beta(1024, X, 1, -2, &idum); for(i=1;i<=1024;i++){ printf("%f\n",X[i]); } return 0; } -- Johan Fredrik Øhman
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@physics.mcgill.ca Montreal, QC Canada H3A 2T8 GPG Fingerprint: 06A9 9553 78BE 16DB 407B FFCA 9BFA B6FF FFD3 2989
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
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 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@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@sourceforge.net _______________________________________________ Numpy-discussion mailing list Numpy-discussion@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/numpy-discussion
participants (2)
-
Johan Fredrik Øhman
-
Scott Ransom