FFT outputs differ between matlab fft function and hand coded FFT
조회 수: 2 (최근 30일)
이전 댓글 표시
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;
}
Any help on this is appreciated.
Thank you in advance.
댓글 수: 0
답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Spectral Measurements에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!