[Numpy-discussion] real_fft

12 May 2002

      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