如下:
/*******************************************************************************
下面的算法是FFT算法
********************************************************************************/
int **gFFTBitTable = NULL;
const int MaxFastBits = 16; //2^16=65536=64K
bool IsPowerOfTwo (int x)//判断X是否是2的N次方的函数,绝对高效的代码
{
if ( x < 2 )
return false;
if ( x & (x-1) ) /* 判断X是否是2的N次方,高效! */
return false;
return true;
}
int NumberOfBitsNeeded (int PowerOfTwo)
/* 已知PowerOfTwo是2的N次方,该函数就是求n
* 如果是输入参数是4,则返回2;如果是8,返回3,如果是16,返回4
*/
{
int i;
if ( PowerOfTwo < 2 ) {
fprintf(stderr, "Error: FFT called with size %d/n", PowerOfTwo);
exit(1);
}
for ( i=0; ; i++ )
if ( PowerOfTwo & (1 << i) )
return i;
}
int ReverseBits ( int index, int NumBits )
/* 位反转,比如:
* 输入,index=0x80,NumBits=16,则执行后返回值为0x01
* 输入,index=0x82,NumBits=16,则执行后返回值为0x41
*/
{
int i, rev;
for ( i=rev=0; i < NumBits; i++ ) {
rev = (rev << 1) | (index & 1);
index >>= 1;
}
return rev;
}
void InitFFT()
{
gFFTBitTable = new int *[MaxFastBits];
int len=2;
for(int b=1; b<=MaxFastBits; b++) {
gFFTBitTable[b-1] = new int[len];
for(int i=0; i<len; i++)
gFFTBitTable[b-1] = ReverseBits(i, b);
len <<= 1;
}
}
inline int FastReverseBits(int i, int NumBits)
{
if (NumBits <= MaxFastBits)
return gFFTBitTable[NumBits-1];
else
return ReverseBits(i, NumBits);
}
/*
* Complex Fast Fourier Transform
* 复数快速傅里叶变换
*/
void FFT2 (
int NumSamples,
int InverseTransform, //逆变换
float *RealIn,
float *ImagIn,
float *RealOut,
float *ImagOut )
{
int NumBits; /* 保存复数需要的位数 */
int i, j, k, n;
int BlockSize, BlockEnd;
double angle_numerator = 2.0 * M_PI;
float tr, ti; /* 临时变量实部,虚部 */
if ( !IsPowerOfTwo(NumSamples) ) {
fprintf (stderr, "%d 不是2的n次方/n", NumSamples);
exit(1);
}
if (!gFFTBitTable)
InitFFT();
if ( InverseTransform )
angle_numerator = -angle_numerator;
NumBits = NumberOfBitsNeeded ( NumSamples );
/*
** Do simultaneous data copy and bit-reversal ordering into outputs...
*/
if (ImagIn==NULL)
{
for ( i=0; i < NumSamples; i++ ) {
j = FastReverseBits ( i, NumBits );
RealOut[j] = RealIn;
ImagOut[j] = 0;
}
}
else
{
for ( i=0; i < NumSamples; i++ ) {
j = FastReverseBits ( i, NumBits );
RealOut[j] = RealIn;
ImagOut[j] = ImagIn;
}
}
/*
** Do the FFT itself...
*/
BlockEnd = 1;
for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 ) {
double delta_angle = angle_numerator / (double)BlockSize;
float sm2 = sin ( -delta_angle-delta_angle );
float sm1 = sin ( -delta_angle );
float cm2 = cos ( -delta_angle-delta_angle );
float cm1 = cos ( -delta_angle );
float w = cm1+cm1;
float ar0, ar1, ar2, ai0, ai1, ai2;
for ( i=0; i < NumSamples; i += BlockSize ) {
ar2 = cm2;
ar1 = cm1;
ai2 = sm2;
ai1 = sm1;
for ( j=i, n=0; n < BlockEnd; j++, n++ ) {
ar0 = w*ar1 - ar2;
ar2 = ar1;
ar1 = ar0;
ai0 = w*ai1 - ai2;
ai2 = ai1;
ai1 = ai0;
k = j + BlockEnd;
tr = ar0*RealOut[k] - ai0*ImagOut[k];
ti = ar0*ImagOut[k] + ai0*RealOut[k];
RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;
RealOut[j] += tr;
ImagOut[j] += ti;
}
}
BlockEnd = BlockSize;
}
/*
** Need to normalize if inverse transform...
*/
if ( InverseTransform )
{
float denom = (float)NumSamples;
for ( i=0; i < NumSamples; i++ ) {
RealOut /= denom;
ImagOut /= denom;
}
}
}
/*
* Real Fast Fourier Transform
*
* This function was based on the code in Numerical Recipes in C.
* In Num. Rec., the inner loop is based on a single 1-based array
* of interleaved real and imaginary numbers. Because we have two
* separate zero-based arrays, our indices are quite different.
* Here is the correspondence between Num. Rec. indices and our indices:
*
* i1 <-> real
* i2 <-> imag
* i3 <-> real[n/2-i]
* i4 <-> imag[n/2-i]
*/
void RealFFT(int NumSamples,float *RealIn,float *RealOut,float *ImagOut)
{
int Half = NumSamples>>1;
int i;
float theta = M_PI / float(Half);
float *tmpReal = new float[Half];
float *tmpImag = new float[Half];
for(i=0; i<Half; i++) {
tmpReal = RealIn[ i<<1 ];
tmpImag = RealIn[ (i<<1) +1 ];
}
FFT2(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
float wtemp = float(sin(0.5*theta));
float wpr = -2.0 * wtemp * wtemp;
float wpi = float(sin(theta));
float wr = 1.0+wpr;
float wi = wpi;
int i3;
float h1r, h1i, h2r, h2i;
for(i=1; i<Half/2; i++) {
i3 = Half-i;
h1r = 0.5*(RealOut+RealOut[i3]);
h1i = 0.5*(ImagOut-ImagOut[i3]);
h2r = 0.5*(ImagOut+ImagOut[i3]);
h2i = -0.5*(RealOut-RealOut[i3]);
RealOut=h1r+wr*h2r-wi*h2i;
ImagOut=h1i+wr*h2i+wi*h2r;
RealOut[i3]=h1r-wr*h2r+wi*h2i;
ImagOut[i3]=-h1i+wr*h2i+wi*h2r;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
RealOut[0] = (h1r=RealOut[0])+ImagOut[0];
ImagOut[0] = h1r-ImagOut[0];
delete []tmpReal;
delete []tmpImag;
}
/*
* PowerSpectrum
* 功率谱
*
* This function computes the same as RealFFT, above, but
* adds the squares of the real and imaginary part of each
* coefficient, extracting the power and throwing away the
* phase.
*
* For speed, it does not call RealFFT, but duplicates some
* of its code.
* 为了速度,这里没有调用RealFFT,但是复制了他的一些代码
*/
void PowerSpectrum(
int NumSamples,
float *In,
float *Out )
{
int Half = NumSamples/2;
int i;
float theta = M_PI / Half;
float *tmpReal = new float[Half];
float *tmpImag = new float[Half];
float *RealOut = new float[Half];
float *ImagOut = new float[Half];
for(i=0; i<Half; i++) {
tmpReal = In[i<<1];
tmpImag = In[(i<<1) +1];
}
FFT2(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);
float wtemp = float(sin(0.5*theta));
float wpr = -2.0 * wtemp * wtemp;
float wpi = float(sin(theta));
float wr = 1.0+wpr;
float wi = wpi;
int i3;
float h1r, h1i, h2r, h2i, rt, it;
for(i=1; i<Half>>1; i++) {
i3 = Half-i;
h1r = 0.5*(RealOut+RealOut[i3]);
h1i = 0.5*(ImagOut-ImagOut[i3]);
h2r = 0.5*(ImagOut+ImagOut[i3]);
h2i = -0.5*(RealOut-RealOut[i3]);
rt = h1r+wr*h2r-wi*h2i;
it = h1i+wr*h2i+wi*h2r;
Out = rt*rt + it*it;
rt = h1r-wr*h2r+wi*h2i;
it = -h1i+wr*h2i+wi*h2r;
Out[i3] = rt*rt + it*it;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
rt = (h1r=RealOut[0])+ImagOut[0];
it = h1r-ImagOut[0];
Out[0] = rt*rt + it*it;
rt = RealOut[Half>>1];
it = ImagOut[Half>>1];
Out[Half>>1] = rt*rt + it*it;
delete[] tmpReal;
delete[] tmpImag;
delete[] RealOut;
delete[] ImagOut;
}
