FFT outputs differ between matlab fft function and hand coded FFT

2 views (last 30 days)
Hello,
I am working on writing a C-code that does the same function as the matlab built-in FFT function. For the implementation, I perform the following :
  1. Bit reversal of the input arrays.
  2. Decimation in time implementation of FFT.
The problem that I am facing is that when I pass a set of inputs as a numbers, the output is the same as matlab fft output. However if I pass an input such as the one below, there is difference between the C- code output and Matlab fft output.:
for(i=0; i<N; i++)
{
(X[i]).real = cos(i*5*2*(pi/N));//i+1;
(X[i]).imag = 0.0;
printf("%4d%15.5f\t%15.5f\n",i,(X[i]).real, (X[i]).imag);
}
Here is a snippet of the code : (Reference : http://dsp-book.narod.ru/CSD.pdf)
/*
/****************************************************************/
/* Function fft(complex *X, int M) */
/* */
/* This is an elementary, complex, radix 2, decimation in */
/* time FFT. The computations are performed "in place" and */
/* the output overwrites the input array. */
/****************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <FFT_header.h>
#include <math.h>
void fft(complex *X, int M)
/* X is an array of N = 2**M complex points. */
{
complex temp1; /* temporary storage complex variable */
complex W; /* e**(-j 2 pi/ N) */
complex U; /* Twiddle factor W**k */
int i,j,k; /* loop indexes */
int id; /* Index of lower point in butterfly */
int N = 1 << M; /* Number of points for FFT */ /* 1 is left shifted*/
int N2 = N/2;
int L; /* FFT stage */
int LE; /* Number of points in sub DFT at stage L, */
/* and offset to next DFT in stage */
int LE1; /* Number of butterflies in one DFT at*/
/* stage L. Also is offset to lower */
/* point in butterfly at stage L */
float pi = 3.141592653589793;//3.1415926535897;
/*==============================================================*/
/* Rearrange input array in bit-reversed order */
/* */
/* The index j is the bit reversed value of i. Since 0 -> 0 */
/* and N-1 -> N-1 under bit-reversal, these two reversals are */
/* skipped. */
j = 0;
for(i=1; i<(N-1); i++)
{
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Increment bit-reversed counter for j by adding 1 to msb and */
/* propagating carries from left to right. */
k = N2; /* k is 1 in msb, 0 elsewhere */
/*--------------------------------------------------------------*/
/* Propagate carry from left to right */
while(k<=j) /* Propagate carry if bit is 1 */
{
j = j - k; /* Bit tested is 1, so clear it. */
k = k/2; /* Set up 1 for next bit to right. */
}
j = j+k; /* Change 1st 0 from left to 1 */
/*--------------------------------------------------------------*/
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
/* Swap samples at locations i and j if not previously swapped.*/
if(i<j) /* Test if previously swapped. */
{
temp1.real = (X[j]).real;
temp1.imag = (X[j]).imag;
(X[j]).real = (X[i]).real;
(X[j]).imag = (X[i]).imag;
(X[i]).real = temp1.real;
(X[i]).imag = temp1.imag;
}
/*++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/
}
/*==============================================================*/
/* Do M stages of butterflies */
for(L=1; L<= M; L++)
{
LE = 1 << L; /* LE = 2**L = points in sub DFT */
LE1 = LE/2; /* Number of butterflies in sub-DFT */
U.real = 1.0;
U.imag = 0.0; /* U = 1 + j 0 */
W.real = cos(pi/LE1);
W.imag = - sin(pi/LE1); /* W = e**(-j 2 pi/LE) */
/*--------------------------------------------------------------*/
/* Do butterflies for L-th stage */
for(j=0; j<LE1; j++) /* Do the LE1 butterflies per sub DFT*/
{
/*..............................................................*/
/* Compute butterflies that use same W**k */
for(i=j; i<N; i += LE)
{
id = i + LE1; /* Index of lower point in butterfly */
temp1.real = (X[id]).real*U.real - (X[id]).imag*U.imag;
temp1.imag = (X[id]).imag*U.real + (X[id]).real*U.imag;
(X[id]).real = (X[i]).real - temp1.real;
(X[id]).imag = (X[i]).imag - temp1.imag;
(X[i]).real = (X[i]).real + temp1.real;
(X[i]).imag = (X[i]).imag + temp1.imag;
}
/*..............................................................*/
/* Recursively compute W**k as W*W**(k-1) = W*U */
temp1.real = U.real*W.real - U.imag*W.imag;
U.imag = U.real*W.imag + U.imag*W.real;
U.real = temp1.real;
/*..............................................................*/
}
/*--------------------------------------------------------------*/
}
return;
}
Any help on this is appreciated.
Thank you in advance.

Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!