#include <stdlib.h>
#include <stdio.h>
#include <math.h> #define M_PI 3.14159265 int **gFFTBitTable = NULL;
const int MaxFastBits = 16; int IsPowerOfTwo(int x)
{
if (x < 2)
return false; if (x & (x - 1)) /* Thanks to 'byang' for this cute trick! */
return false; return true;
}
//FFT size
int NumberOfBitsNeeded(int PowerOfTwo)
{
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; //返回FFT size 2的幂
} int ReverseBits(int index, int NumBits)
{
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][i] = ReverseBits(i, b); len <<= 1;
}
} inline int FastReverseBits(int i, int NumBits)
{
if (NumBits <= MaxFastBits)
return gFFTBitTable[NumBits - 1][i];
else
return ReverseBits(i, NumBits);
} /*
* Complex Fast Fourier Transform
*/ void FFT(int NumSamples,
bool InverseTransform,
float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
{
int NumBits; /* Number of bits needed to store indices */
int i, j, k, n;
int BlockSize, BlockEnd; double angle_numerator = 2.0 * M_PI;
float tr, ti; /* temp real, temp imaginary */ if (!IsPowerOfTwo(NumSamples)) {
fprintf(stderr, "%d is not a power of two\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...
*/ for (i = 0; i < NumSamples; i++) {
j = FastReverseBits(i, NumBits);
RealOut[j] = RealIn[i];
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
} /*
** Do the FFT itself...
*/ BlockEnd = 1;
for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) { double delta_angle = angle_numerator / (double) BlockSize; float sm2 = sin(-2 * delta_angle);
float sm1 = sin(-delta_angle);
float cm2 = cos(-2 * delta_angle);
float cm1 = cos(-delta_angle);
float w = 2 * 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[i] /= denom;
ImagOut[i] /= 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[i]
* i2 <-> imag[i]
* i3 <-> real[n/2-i]
* i4 <-> imag[n/2-i]
*/ void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
{
int Half = NumSamples / 2;
int i; float theta = M_PI / Half; float *tmpReal = new float[Half];
float *tmpImag = new float[Half]; for (i = 0; i < Half; i++) {
tmpReal[i] = RealIn[2 * i];
tmpImag[i] = RealIn[2 * i + 1];
} FFT(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[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]); RealOut[i] = h1r + wr * h2r - wi * h2i;
ImagOut[i] = 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;
}
如下为一个使用了FFT变换的功率谱计算例子。
/*
* 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.
*/ 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[i] = In[2 * i];
tmpImag[i] = In[2 * i + 1];
} FFT(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 / 2; i++) { i3 = Half - i; h1r = 0.5 * (RealOut[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]); rt = h1r + wr * h2r - wi * h2i;
it = h1i + wr * h2i + wi * h2r; Out[i] = 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 / 2];
it = ImagOut[Half / 2];
Out[Half / 2] = rt * rt + it * it; delete[]tmpReal;
delete[]tmpImag;
delete[]RealOut;
delete[]ImagOut;
}
#include <stdio.h>
#include <math.h> #define M_PI 3.14159265 int **gFFTBitTable = NULL;
const int MaxFastBits = 16; int IsPowerOfTwo(int x)
{
if (x < 2)
return false; if (x & (x - 1)) /* Thanks to 'byang' for this cute trick! */
return false; return true;
}
//FFT size
int NumberOfBitsNeeded(int PowerOfTwo)
{
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; //返回FFT size 2的幂
} int ReverseBits(int index, int NumBits)
{
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][i] = ReverseBits(i, b); len <<= 1;
}
} inline int FastReverseBits(int i, int NumBits)
{
if (NumBits <= MaxFastBits)
return gFFTBitTable[NumBits - 1][i];
else
return ReverseBits(i, NumBits);
} /*
* Complex Fast Fourier Transform
*/ void FFT(int NumSamples,
bool InverseTransform,
float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
{
int NumBits; /* Number of bits needed to store indices */
int i, j, k, n;
int BlockSize, BlockEnd; double angle_numerator = 2.0 * M_PI;
float tr, ti; /* temp real, temp imaginary */ if (!IsPowerOfTwo(NumSamples)) {
fprintf(stderr, "%d is not a power of two\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...
*/ for (i = 0; i < NumSamples; i++) {
j = FastReverseBits(i, NumBits);
RealOut[j] = RealIn[i];
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
} /*
** Do the FFT itself...
*/ BlockEnd = 1;
for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) { double delta_angle = angle_numerator / (double) BlockSize; float sm2 = sin(-2 * delta_angle);
float sm1 = sin(-delta_angle);
float cm2 = cos(-2 * delta_angle);
float cm1 = cos(-delta_angle);
float w = 2 * 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[i] /= denom;
ImagOut[i] /= 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[i]
* i2 <-> imag[i]
* i3 <-> real[n/2-i]
* i4 <-> imag[n/2-i]
*/ void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
{
int Half = NumSamples / 2;
int i; float theta = M_PI / Half; float *tmpReal = new float[Half];
float *tmpImag = new float[Half]; for (i = 0; i < Half; i++) {
tmpReal[i] = RealIn[2 * i];
tmpImag[i] = RealIn[2 * i + 1];
} FFT(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[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]); RealOut[i] = h1r + wr * h2r - wi * h2i;
ImagOut[i] = 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;
}
如下为一个使用了FFT变换的功率谱计算例子。
/*
* 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.
*/ 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[i] = In[2 * i];
tmpImag[i] = In[2 * i + 1];
} FFT(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 / 2; i++) { i3 = Half - i; h1r = 0.5 * (RealOut[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]); rt = h1r + wr * h2r - wi * h2i;
it = h1i + wr * h2i + wi * h2r; Out[i] = 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 / 2];
it = ImagOut[Half / 2];
Out[Half / 2] = rt * rt + it * it; delete[]tmpReal;
delete[]tmpImag;
delete[]RealOut;
delete[]ImagOut;
}
原理介绍
一个叠加了白噪声的正弦波序列，序列长度为1024。根据数字信号处理理论可知，如果两点之间的采样间隔为1ms的话，对该序列进行傅立叶变换时的基频△f=1/1024*1000HZ，同时由上图可知为了滤除叠加在低频正弦波上的高频噪声，需要先将上述序列信号进行时域到频域的转换，然后在得到的频域序列中去除高频部分。
由于有用的正弦波信号为△f，那么需要滤除f>△f的信号。可得K=f/△f。具体的做法是将变换成频域后的序列中的0-K和1024-K到1024这两段序列保留，中间的K到1024-K这段序列使之为0。
具体的FFT操作有FFT变换和FFT逆变换。对于上述这两种变换又有复述形式和指数形式。
对于复数形式，需要将测量得到的1024点数据作为时域数据的实部，然后另外构建一个全部为0的1024点大小的虚部数组。
将上述时域的实部数组和虚部数组最为参数输入复数形式的FFT变换，然后得到代表频域信号的实部数组和虚部数组。
对该实部数组和虚部数组中K到1024-K之间的部分全部取0，然后得到新的代表频域信号的实部数组和虚部数组，然后带入复述形式的FFT逆变换。
变换后得到代表滤波后的时域信号的实部数组和虚部数组，然后将两个数组中对应得实部和虚部取模，计算后的模值根据实部的正负号取相应的正负号。计算得到的1024个模值既是滤波后的结果。
http://www.a1vbcode.com/downloads/fft.zip