FFT outputs differ between matlab fft function and hand coded FFT

Question

A VP 2019-1-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/442186-fft-outputs-differ-between-matlab-fft-function-and-hand-coded-fft

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 :

Bit reversal of the input arrays.
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;
}