Danijel - 9 months ago 68

C Question

I'm calculating discrete Fourier transform using Numerical Recipes and (to confirm the result) Matlab. Working with real values only.

My Matlab code

`in(1)=0.0;`

in(2)=1.0;

in(3)=2.0;

in(4)=3.0;

in(5)=4.0;

in(6)=5.0;

in(7)=6.0;

in(8)=7.0;

out = fft(in);

gives me

`out =`

28.0000 + 0.0000i

-4.0000 + 9.6569i

-4.0000 + 4.0000i

-4.0000 + 1.6569i

-4.0000 + 0.0000i

-4.0000 - 1.6569i

-4.0000 - 4.0000i

-4.0000 - 9.6569i

What

`data`

`void realft( float data[], unsigned long n, int isign ){...}`

in order to get the same output as Matlab?

From this NR forum link, I found out that input to

`realft`

`in[0] = 0`

Test code:

`#include <stdio.h>`

#define LEN 8

int main()

{

float inout[LEN+1];

inout[0] = 0.0;

inout[1] = 0.0;

inout[2] = 1.0;

inout[3] = 2.0;

inout[4] = 3.0;

inout[5] = 4.0;

inout[6] = 5.0;

inout[7] = 6.0;

inout[8] = 7.0;

realft( inout, LEN, 1 );

for( unsigned int i=0; i<LEN+1; i=i+1)

printf("%15.10f \n",inout[i]);

return 0;

}

Output from the test code is:

`0.00000000`

28.00000000

-4.00000000

-4.00000000

-9.65685463

-4.00000000

-4.00000000

-4.00000000

-1.65685427

which is similar, but different from Matlab.

Here is the Fourier transform code.

`realft`

`four1`

`#include <math.h>`

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(float data[], unsigned long nn, int isign)

{

unsigned long n, mmax, m, j, istep, i;

double wtemp, wr, wpr, wpi, wi, theta; /* Double precision for the trigonometric recurrences. */

float tempr, tempi;

n=nn << 1;

j=1;

for (i=1;i<n;i+=2)

{

/* This is the bit-reversal section of the routine. */

if (j > i)

{

SWAP(data[j],data[i]); /* Exchange the two complex numbers.*/

SWAP(data[j+1],data[i+1]);

}

m=nn;

while (m >= 2 && j > m)

{

j -= m;

m >>= 1;

}

j += m;

}

/* Here begins the Danielson-Lanczos section of the routine. */

mmax=2;

/* Outer loop executed log 2 nn times. */

while (n > mmax)

{

istep=mmax << 1;

theta=isign*(6.28318530717959/mmax); /* Initialize the trigonometric recurrence.*/

wtemp=sin(0.5*theta);

wpr = -2.0*wtemp*wtemp;

wpi=sin(theta);

wr=1.0;

wi=0.0;

/* Here are the two nested inner loops. */

for (m=1;m<mmax;m+=2)

{

for (i=m;i<=n;i+=istep)

{

j=i+mmax; /*This is the Danielson-Lanczos formula: */

tempr=wr*data[j]-wi*data[j+1];

tempi=wr*data[j+1]+wi*data[j];

data[j]=data[i]-tempr;

data[j+1]=data[i+1]-tempi;

data[i] += tempr;

data[i+1] += tempi;

}

wr=(wtemp=wr)*wpr-wi*wpi+wr; /* Trigonometric recurrence. */

wi=wi*wpr+wtemp*wpi+wi;

}

mmax=istep;

}

}

/* Calculates the Fourier transform of a set of n real-valued data points.

Replaces this data (which is stored in array data[1..n] ) by the positive

frequency half of its complex Fourier transform. The real-valued first

and last components of the complex transform are returned as elements

data[1] and data[2] , respectively. n must be a power of 2. This routine

also calculates the inverse transform of a complex data array if it is

the transform of real data. (Result in this case must be multiplied by 2/n.) */

void realft( float* data, unsigned long n, int isign )

{

unsigned long i,i1,i2,i3,i4,np3;

float c1=0.5,c2,h1r,h1i,h2r,h2i;

double wr,wi,wpr,wpi,wtemp,theta; /* Double precision for the trigonometric recurrences.*/

theta=3.141592653589793/(double) (n>>1); /* Initialize the recurrence.*/

if (isign == 1)

{

c2 = -0.5;

four1(data,n>>1,1); /* The forward transform is here. */

}

else

{

c2=0.5; /* Otherwise set up for an inverse transform. */

theta = -theta;

}

wtemp=sin(0.5*theta);

wpr = -2.0*wtemp*wtemp;

wpi=sin(theta);

wr=1.0+wpr;

wi=wpi;

np3=n+3;

for (i=2;i<=(n>>2);i++) /* Case i=1 done separately below.*/

{

i4=1+(i3=np3-(i2=1+(i1=i+i-1)));

h1r=c1*(data[i1]+data[i3]); /* The two separate transforms are separated out of data. */

h1i=c1*(data[i2]-data[i4]);

h2r = -c2*(data[i2]+data[i4]);

h2i=c2*(data[i1]-data[i3]);

data[i1]=h1r+wr*h2r-wi*h2i; /* Here they are recombined to form the true transform of the original real data. */

data[i2]=h1i+wr*h2i+wi*h2r;

data[i3]=h1r-wr*h2r+wi*h2i;

data[i4] = -h1i+wr*h2i+wi*h2r;

wr=(wtemp=wr)*wpr-wi*wpi+wr; /* The recurrence.*/

wi=wi*wpr+wtemp*wpi+wi;

}

if (isign == 1)

{

data[1] = (h1r=data[1])+data[2]; /* Squeeze the first and last data together to get them all within theoriginal array. */

data[2] = h1r-data[2];

}

else

{

data[1]=c1*((h1r=data[1])+data[2]);

data[2]=c1*(h1r-data[2]);

four1(data,n>>1,-1); /* This is the inverse transform for the case isign=-1.*/

}

}

Answer Source

The first aspect that is different between Matlab and Numerical Recipes' implementations, other than the off-by-one indexing that you've already noticed and accounted for, is that they are based on a slightly different definition of the FFT. More specifically, Matlab uses a negative complex exponential for the forward transform while Numerical Recipes uses a positive complex exponential. Correspondingly Numerical Recipes' implementation would give a result that is the complex conjugate of that from Matlab.

The other thing is that those implementations produce a result in a different packing order, and Numerical Recipes' implementation only outputs the non-redundant lower half of the spectrum. Graphically, this mapping could be represented (after changing the sign of the imaginary parts in accordance with the previous point with respect to complex conjugation) by the following graph: