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

Answer

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:

enter image description here

Comments