// fft.c
#include <math.h>
void fft_init (double *pFFTCoeffs,
int FFTSize)
{
int i;
// Generate 3/4 Sine look-up-table
for (i = 0; i < ((3 * FFTSize) >> 2); i++) {
*pFFTCoeffs++ = sin ((((double)i) * 6.283185307179586476925286766559) / ((double)FFTSize));
}
} // End of fft_init()
void rfft (double RealData[],
double ImagData[],
const double *pFFTCoeffs,
int FFTSize,
int Log2Size)
{
const double *pFFTSineCoeffs;
const double *pFFTCosineCoeffs;
int i, j, k, Stride, BflyCounter, g, h, HalfFFTSize;
int Angle, AngleInc; // Angle step thro sin & cos tables
double RealTemp, ImagTemp, Cos, Sin;
pFFTSineCoeffs = pFFTCoeffs;
pFFTCosineCoeffs = pFFTCoeffs + (FFTSize >> 2);
HalfFFTSize = FFTSize >> 1;
Stride = j = HalfFFTSize;
Angle = h = 0;
// First stage
for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++) {
RealTemp = RealData[h] - RealData[j];
RealData[h] = RealData[h] + RealData[j];
ImagData[h] = 0.0; // Clear imaginary part
RealData[j] = pFFTCosineCoeffs[Angle] * RealTemp;
ImagData[j] = -pFFTSineCoeffs[Angle] * RealTemp;
Angle++;
h++;
j++;
}
AngleInc = 2;
for (i = 1; i < (Log2Size - 1); i++) { // Middle stages
k = Stride;
Stride >>= 1;
Angle = 0;
for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++) {
Cos = pFFTCosineCoeffs[Angle];
Sin = pFFTSineCoeffs[Angle];
Angle += AngleInc;
h = BflyCounter;
j = h + Stride;
for (g = k; g <= FFTSize; g += k, h += k, j += k) {
RealTemp = RealData[h] - RealData[j];
ImagTemp = ImagData[h] - ImagData[j];
RealData[h] = RealData[h] + RealData[j];
ImagData[h] = ImagData[h] + ImagData[j];
RealData[j] = Cos * RealTemp + Sin * ImagTemp;
ImagData[j] = Cos * ImagTemp - Sin * RealTemp;
}
}
AngleInc <<= 1;
}
// Final stage
for (h = 0, j = 1; h < FFTSize; h += 2, j += 2) {
RealTemp = RealData[h] - RealData[j];
ImagTemp = ImagData[h] - ImagData[j];
RealData[h] = RealData[h] + RealData[j];
ImagData[h] = ImagData[h] + ImagData[j];
RealData[j] = RealTemp; /* Cos = 1, sin = 0 */
ImagData[j] = ImagTemp;
}
// Reorder scrambled data
for (j = 0, i = 0; j < FFTSize; j++) {
if (j < i) {
RealTemp = RealData[i]; RealData[i] = RealData[j]; RealData[j] = RealTemp;
ImagTemp = ImagData[i]; ImagData[i] = ImagData[j]; ImagData[j] = ImagTemp;
}
k = HalfFFTSize;
while ((k <= i) && (k >= 1)) {
i -= k;
k >>= 1;
}
i += k;
}
} // End of rfft()
void cfft (double RealData[],
double ImagData[],
const double *pFFTCoeffs,
int FFTSize,
int Log2Size)
{
int i, j, k, Stride, BflyCounter, g, h;
int Angle, AngleInc; // Angle step thro sin & cos tables
double ImagTemp, RealTemp, Cos, Sin;
const double *pFFTSineCoeffs = pFFTCoeffs;
const double *pFFTCosineCoeffs = pFFTCoeffs + (FFTSize >> 2);
for (i = 0, Stride = FFTSize, AngleInc = 1; i < Log2Size; i++) {
k = Stride;
Stride >>= 1;
Angle = 0;
for (BflyCounter = 0; BflyCounter < Stride; BflyCounter++) {
Cos = pFFTCosineCoeffs[Angle];
Sin = pFFTSineCoeffs[Angle];
Angle += AngleInc;
h = BflyCounter;
j = h + Stride;
for (g = k; g <= FFTSize; g += k, h += k, j += k) {
RealTemp = RealData[h] - RealData[j];
ImagTemp = ImagData[h] - ImagData[j];
RealData[h] = RealData[h] + RealData[j];
ImagData[h] = ImagData[h] + ImagData[j];
RealData[j] = Cos * RealTemp + Sin * ImagTemp;
ImagData[j] = Cos * ImagTemp - Sin * RealTemp;
}
}
AngleInc <<= 1;
}
// Reorder scrambled data
for (j = 0, i = 0; j < FFTSize; j++) {
if (j < i) {
RealTemp = RealData[i]; RealData[i] = RealData[j]; RealData[j] = RealTemp;
ImagTemp = ImagData[i]; ImagData[i] = ImagData[j]; ImagData[j] = ImagTemp;
}
k = FFTSize >> 1;
while ((k <= i) && (k >= 1)) {
i -= k;
k >>= 1;
}
i += k;
}
} // End of cfft()
These functions, and more, are available in the SigLib DSP Library.
If you have found this solution useful then please do hit the Google (+1) button so that others may be able to find it as well.
Numerix-DSP Libraries : http://www.numerix-dsp.com/eval/
