关于FFT算法的原理这里就不多说了，具体参考有关书籍。
N点DFT的运算量
|
|
复数乘法 |
复数加法 |
| 一个X(k) |
N |
N-1 |
| N个X(k)（N点DFT） |
N*N |
N(N-1) |
N点FFT的运算量
|
|
复数乘法 |
复数加法 |
| N个X(k) |
(N/2)*log2N |
N*log2N |
| 参考《数字信号处理》（西电版）P115的蝶形运算公式，可以推导出用计算机程序实现的公式。由于word中输入公式太麻烦了，这里先将推导过程写在纸上，然后拍照上传到这里。 |
|
|
根据上面的推导结果，就可以用程序来实现FFT了。看下面的例程。
#include "math.h"
#define PI 3.14159
#define N 128
//real：实数
//imaginary：虚数
#pragma DATA_SECTION(Input,"ExRamMem");//这是DSP处理器(CCSv4开发环境)中的一种操作内存的方法
int Input[N] = {0.0};//用于存放待处理数据
#pragma DATA_SECTION(SinTab,"ExRamMem");
float SinTab[N] = {0.0};//用于FFT运算的正弦表
#pragma DATA_SECTION(CosTab,"ExRamMem");
float CosTab[N] = {0.0};//用于FFT运算的余弦表
#pragma DATA_SECTION(Real,"ExRamMem");
float Real[N] = {0.0};//待处理数据的实部
#pragma DATA_SECTION(Imag,"ExRamMem");
float Imag[N] = {0.0};//待处理数据的虚部
#pragma DATA_SECTION(Mag,"ExRamMem");
float Mag[N] = {0.0};//FFT的幅度谱
//函数功能：生成一个正弦波信号
void GenSin(void)
{
int i;
for ( i=0;i<N;i++ )
{
INPUT[i]=sin(PI*2*i/N*3)*1024;
}
}
//函数功能：计算用于FFT运算的正、余弦表。
//(me)注意：旋转因子和正、余弦表不是同一个概念。
void TwiddleCmpt(void)
{
int i;
for ( i=0; i<N; i++ )
{
SinTab[i] = sin(PI*2*i/N);
CosTab[i] = cos(PI*2*i/N);
}
}
//输入参数：real：待处理的数据的实部
// imag：待处理的数据的虚部
//返回参数：暂时就返回OK
//函数功能：计算实部real和虚部imag的N点FFT。
//注意：因为FFT的运算为了节约内存，而将下一级的运算结果直接覆盖上一级，所以下面的程序中要采取temp_R、temp_I、temp_R_kb暂时存放一下。
int FFT(float real[N],float imag[N])
{
//调试发现：定义long型变量与定义int型变量时的计算速度居然不一样，定义long型变量的计算速度要慢一点，即使是参与运算的相关变量都是long型的。但是没办法当FFT点数较大时，int型变量的范围是不够用的。
long i,j=0,k,b,p,Nv2,Nm1,bm1,bx2; //Nv2表示N/2，Nm2表示N-1，bm1表示b-1，bx2表示b乘2，PIx2表示PI乘2，定义这五个数是为了减少运算时间。
float temp,temp_R,temp_I,temp_R_kb;
long L,M;
// f=N;
// for(m=1;(f=f/2)!=1;m++) //判断N是否能够被2整除，判断N是2的几次幂。
// {}
Nv2=N/2;
Nm1=N-1;
M=log(N)/log(2);
//这里的倒序采用雷德算法。
//=================== following code invert sequence ===================
for(i=0; i<Nm1; i++)
{
if(i < j)
{
temp = real[j];
real[j] = real[i];
real[i] = temp;
}
k = Nv2;
while(k <= j)
{
j = j - k;
k = k/2;
}
j = j + k;
}
//===================following code FFT ===================
for (L=1; L<=M; L++ )//for(1)：控制运算级数
{
b = 1;
i = L - 1;
while (i > 0)
{
b=b*2; i--;
} // b= 2^(L-1)，蝶尖距
bm1 = b-1; //在进行第二层循环前，先将bm1和bx2算出来。
bx2 = b*2;
for (j=0; j<=bm1; j++ )//for(2)：同一级中的蝶群数。
{
p=1; i=M-L;
while (i > 0)
{
p = p*2; i--;
}
p = p*j; // p=pow(2,M-L)*j;
for (k=j; k<N; k=k+bx2)//for(3)：一个蝶群中有多少只蝴蝶
{
temp_R = real[k]; temp_I = imag[k]; temp_R_kb = real[k+b];
real[k] = real[k] + real[k+b] * CosTab[p] + imag[k+b] * SinTab[p];
imag[k] = imag[k] - real[k+b] * SinTab[p] + real[k+b] * CosTab[p];
real[k+b]= temp_R - real[k+b] * CosTab[p] - real[k+b] * SinTab[p]; //因为此时的real[k]已被上一条语句覆盖，但这里还是想用之前的real[k]，所以才事先将real[k]存储在temp_R中，temp_I和temp_R_kb与此类似。
imag[k+b]= temp_I + temp_R_kb * SinTab[p] - real[k+b] * CosTab[p];
} // END for (3)
} // END for (2)
} // END for (1)
for ( i=0; i<Nv2; i++ )//计算幅度谱
{
Mag[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]);//FFT的输出是由实部和虚部组成，要想看频谱，还求实部和虚部平方和再开根号，这样得到就才是幅度谱。
}
asm(" NOP");
return(OK);
}
int fft(void)
{
Uint16 i;
GenSin();
TwiddleCmpt();
for(i=0; i<N; i++)
{
Real[i] = Input[i];//把采集到的N个数据付给FFT的实部。
Imag[i] = 0.0;//虚部置0
Mag[i] = 0.0;//幅度置0
}
FFT(Real, Imag);
return(OK);
}| 正弦波生成函数GenSin()生成波形如下 |
| |
| 经过fft处理后得到的幅度谱如下 注意： 1. 这里得到的幅度谱，其横坐标只是下标值，若想将横坐标转换为频率，还需要进行计算，公式是：f=n*Fs/N，其中n是对应处理结果Mag的下标值，Fs是采样率，N是fft点数，f即是n点对应的频率。 2. 如下图这个处理结果，给人的感觉是只有Mag[3]的值较大，其他都是0，其实有很多点的值不是等于0的，只是他们的值相对Mag[3]太小，所以显得他们是0。 3. 这里的信号源是交流信号，即不含直流分量，所Mag的第0点的值是0。而一般情况下AD的采样结果都是正值，也就是采样信号中会含有直流分量，这时候在第0点也会有值，只不过该点对应的频率是0。 |
| |
| CCSv4软件自带的频谱分析工具分析的fft处理结果如下 |
| |
以N=128(2^7)为例，给出倒序前后的序号。
| 倒序前序号 |
倒序后序号 |
||
| 十进制 |
二进制 |
二进制 |
十进制 |
| 0 |
000 0 000 |
000 0 000 |
0 |
| 1 |
000 0 001 |
100 0 000 |
64 |
| 2 |
000 0 010 |
010 0 000 |
32 |
| 3 |
000 0 011 |
110 0 000 |
96 |
| 4 |
000 0 100 |
001 0 000 |
16 |
| …… |
…… |
…… |
…… |
| 126 |
111 1 110 |
011 1 111 |
63 |
| 127 |
111 1 111 |
111 1 111 |
127 |
| 倒序除了上面的雷德算法也外，可以采用下面这种更为直接的方法进行，不过这种方法的效率显然没有雷德算法高。 |
|
