Danijel - 1 month ago 17
C Question

# How to compare Real Fourier Transform implementation from Numerical Recipes to Matlab fft?

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`
input do I need to send to Numerical recipes

``````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`
needs to be shifted one place, so I am using
`in[0] = 0`
, and consequently input is 1 element larger than N.

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`
is taken from Numerical Recipes (uses
`four1`
routine):

``````#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.*/
}
}
``````