void RADIX_2_DIT_FFT(double *fftBuffer, long fftFrameSize, long sign)
/*
FFT routine, (C)1996 S.M.Bernsee. Sign = -1 is FFT, 1 is iFFT (inverse)
Fills fftBuffer[0...2*fftFrameSize-1] with the Fourier transform of the
time domain data in fftBuffer[0...2*fftFrameSize-1]. The FFT array takes
and returns the cosine and sine parts in an interleaved manner, ie.
fftBuffer[0] = cosPart[0], fftBuffer[1] = sinPart[0], asf. fftFrameSize
must be a power of 2. It expects a complex input signal (see footnote 2),
ie. when working with 'common' audio signals our input signal has to be
passed as {in[0],0.,in[1],0.,in[2],0.,...} asf. In that case, the transform
of the frequencies of interest is in fftBuffer[0...fftFrameSize].
2-radix DIT(Decimation In Time)-FFT algorithm.
*/
{
double wr, wi, arg, *p1, *p2, temp;
double tr, ti, ur, ui, *p1r, *p1i, *p2r, *p2i;
long i, bitm, j, le, le2, k;
// re-order x(n) or X[n] in bit reverse order.
// both x(n) and X[n] are complex value, a complex value occupies contineous
// two memory, the first is real part, and the second is imaginary part.
for (i = 2; i < 2*fftFrameSize-2; i += 2) {
for (bitm = 2, j = 0; bitm < 2*fftFrameSize; bitm <<= 1) {
if (i & bitm) j++;
j <<= 1;
}
// swap two point according bit reverse order
if (i < j) {
p1 = fftBuffer+i; p2 = fftBuffer+j;
temp = *p1; *(p1++) = *p2;
*(p2++) = temp; temp = *p1;
*p1 = *p2; *p2 = temp;
}
}
// After re-order in bit reverse rule, then the output will be in normal order.
for (k = 0, le = 2; k < (long)(log(fftFrameSize)/log(2.)+.5); k++) {
le <<= 1;
le2 = le>>1;
ur = 1.0; // cos0
ui = 0.0; // sin0
arg = M_PI / (le2>>1); // 2*PI/(2^r)
wr = cos(arg); // wr = cos(arg)
wi = sign*sin(arg); // wi = sin(sign*arg), if forward FFT, then sign = -1, otherwise sign = 1
for (j = 0; j < le2; j += 2) {
p1r = fftBuffer+j; p1i = p1r+1; // the first point of one butterfly calculation: P1
p2r = p1r+le2; p2i = p2r+1; // the second point of one butterfly calculation: P2
for (i = j; i < 2*fftFrameSize; i += le) {
tr = *p2r * ur - *p2i * ui; // real part of: P2 * WN^r
ti = *p2r * ui + *p2i * ur; // imag part of: P2 * WN^r
*p2r = *p1r - tr; *p2i = *p1i - ti;
*p1r += tr; *p1i += ti;
p1r += le; p1i += le;
p2r += le; p2i += le;
}
tr = ur*wr - ui*wi; // cos(A+arg) = cosAcos(arg) - sinAsin(arg)
ui = ur*wi + ui*wr; // sin(A+arg) = sinAcos(arg) + cosAsin(arg)
ur = tr;
}
}
}
|
楼主
