Digital Signal Processing Reference
In-Depth Information
means of sequential calculation of the butterfly stages with the required twiddle
factors computed first.
/*============================================================
FastFT() - radix-2 decimation-in-time Fast Fourier Transform
Prototype: void FastFT(short Length, short Bits, complex *X);
Arguments: Length - length of FFT
Bits - Length = 2^Bits
X - ptr to array of complex numbers
Return: none (conversion done in place)
============================================================*/
void FastFT(short Length, short Bits, complex *X)
{ long k, i, p, q, /* counters and indices */
M; /* stage of FFT */
complex A, B, V, W, Tmp; /* complex constants */
short
n,
/* counter */
RIndex;
/* reversed ndex */
/* Swap the values in x */
for(n = 0; n < Length; n++)
{ RIndex = RevBits(n,Bits); /* get reversed index */
if (RIndex < n) continue; /* only need to swap half */
Tmp = *(X+n); *(X+n) = *(X+RIndex); *(X+RIndex) = Tmp;
}
/* Implement the butterfly in stages */
M = 2;
while (M <= Length) /* while loop controls stages of FFT */
{ W.re = cos(2*PI / M);
/* twiddle factors */
W.im = sin(2*PI / M);
V.re = 1;
V.im = 0;
for(k = 0; k < M/2; k++) /* index through stages */
{ for(i = 0; i < Length; i += M)
{ p = k + i;
q = p + M / 2;
A = X[p];
B = cmul(X[q],V);
X[p] = cadd(A,B);
/* butterfly operations */
X[q] = csub(A,B);
}
V = cmul(V,W);
}
M = 2 * M;
}
}
Listing 9.1
FastFT
function.
These two functions can now be used to compute the FFT of any input
sequence. All that is required is that the sequence be adjusted to a length that is a
power of two by padding the end of the sequence with zeros. The function can
then be provided with the pointer to the sequence, the length of the sequence, and
the number of bits in the indices. After the function has completed its work, the
array
X
will hold the complex valued FFT. These values will be stored as the real
and imaginary values of the transform, but can be better interpreted as the
