[分享] FFT算法的实现

关于FFT算法全国大学生电子设计竞赛连续出了两年了,07年的音频信号分析仪,09年的音频均衡器也可以用FFT去做.国内的学生最不擅长还是算法，所有网上都找不到相关的赛后优秀论文，所以我在这里给出我已经实现并验证的思路，算法参考了网上给出的一些设计思路和代码，并进行了局部的修改和优化。这里完全套用了FFT算法公式，并不是最为优化的，因为实际的物理量并没有虚部，于是我们强制虚部为零，这样计算出的结果是关于Fs/2对称的，我们可以把一半的实部放虚部去做n/2点的FFT，这样时间上可以减少一半。所以有兴趣的网友可以在我的基础上继续研究优化算法。我也可以利用业余时间给你支持（通过站内邮箱或回帖给我留下联系方法）。
基础的知识网上都有,我就不再讲了,对这个感兴趣的应该都了解一些了.比如采样频率,频率分辨率,以及离散傅立叶变换到快速傅立叶变换的推倒等等,就说说我算法的实现吧.（fs=20kHz，n=512）
首先是倒序，如果用的不是高速的DSP，建议你用一张表去实现，我现在用的是ARM，所以用了一张倒位的表格。我做的是512点的FFT，所以看起来表格挺大的，但是现在ARM的FLASH都挺大，没有任何关系。

FFT的两个公式的实现网上较多的做法是用三级循环的方式,我也采用了这个思路,需要说明的是现在的控制器都有硬件乘法器,里面的乘法计算最好用汇编代码启动硬件乘法器去实现,否则编译器可能会利用软件去实现乘法计算.

for(stage=0; stage级的FFT算法

{

for(nb_index=0; nb_index

{

int tf_index = 0; // The twiddle factor index

for(sb_index=0; sb_index

{

s16 resultMulReCos;


s16 resultMulImCos;

s16 resultMulReSin;

s16 resultMulImSin;


s16 b_index = a_index+s_of_b; // 2nd fft data index

                          

resultMulReCos=MUL_1(cosLUT[tf_index],x_n_re[b_index]); // MUL_1就是汇编代码

resultMulReSin=MUL_1(sinLUT[tf_index],x_n_re[b_index]);

resultMulImCos=MUL_1(cosLUT[tf_index],x_n_im[b_index]);

resultMulImSin=MUL_1(sinLUT[tf_index],x_n_im[b_index]);
         

x_n_re[b_index] = x_n_re[a_index]-resultMulReCos+resultMulImSin;

x_n_im[b_index] = x_n_im[a_index]-resultMulReSin-resultMulImCos;

x_n_re[a_index] = x_n_re[a_index]+resultMulReCos-resultMulImSin;

x_n_im[a_index] = x_n_im[a_index]+resultMulReSin+resultMulImCos;
            

if (((sb_index+1) & (s_of_b-1)) == 0)

a_index = a_index_ref;

else

a_index++;


tf_index += n_of_b;

}

a_index
= ((s_of_b<<1) + a_index) & N_MINUS_1;

a_index_ref = a_index;

}

n_of_b >>= 1;

s_of_b <<= 1;

}
最后一步就是把FFT的结果的实部和虚部换算到实际的物理量
x_n_re =sqrt(x_n_re*x_n_re+x_n_im*x_n_im);
需要注意的是他和实际的物理量之间还有系数关系。当i=0，即为直流分量的时候理论为x_n_re /n，其他分量为x_n_re /(N/2)，但实际上我们一般不关注直流分量，所以就不管他了。有些网友反应直流分量计算不准确，应该就是这个原因了。
来看看效果吧：320×240的TFT液晶，界面不美观，只是为了把频谱显示出来，信号来自麦克风，源来自电脑。

回复 8楼 122402902 的帖子

好牛啊！

回复 15楼 wangjiafu1985 的帖子