用单片机做FFT

1、这次的任务是要用单片机做FFT(快速离散傅利叶变换),我觉得单片机越来越无敌了. 于是,连续几天被埋在时间与空间的斗争中.因为单片机的硬件资源和条件所限,不能像PC上的MatLab那样方便地进行FFT,我们单片机 (SPCE061A)的时钟频率最快约为50MHz,其RAM为2K字,ROM为32K字.从前学数据结构和算法的时候,千方百计想着用时间换空间还是用 空间换时间,现在才发现,那时是站着说话不腰疼,因为无论拿什么换什么,都是行得通的,损失的效率对于PC应用来说(就我现在的水平),是可以忽略不计 的,毕竟不是ACM(此理论的唯一一次失手就是在某次ACM比赛上). 整体来说,单片机最大的弱点还

2、是在于速度太慢,根本不可能进行实时FFT,而手动启动来进行运算的速度恐怕也慢得让人难以接受.为了加速,我几天苦苦挣扎 于算法的优化上.若不是因为数据精度还有一定的要求(汗),不能采用8Bit数据,我肯定会将两个数据压缩到一个字(16位)中去以节省空间.坚决榨干 单片机的每一滴油!我将所有的冗余计算全部转化为变量表示独立计算然后重复利用,而常量,哪怕是局部的,我也将它写成常数表(虽然还不知道如何在C中将常 数表写入ROM区,再汗一个).现在唯一的希望就在,单片机对浮点数运算的支持能力上.从前碰到浮点数我都会转化为整数,舍不得拿单片机来算浮点 数,而这次计算量实在太大了一点儿. 到目前为止,只在G

3、CC上编译通过,并且结果是正确的(与MabLab对照),接下来还是任重而道远啊 贴出主要部分代码:/*/FFT(快速离散傅利叶变换)/本例为8级256点数据/HUST STI OE0402 POPAPPLE/popapple/2007-8-9 18:10:20*/#include stdio.h#define LEVEL 8 /级数#define LENGTH 256 /点数/用结构体来表示复数typedef struct _JComplexfloat real;float image;JComplex,*pJComplex;typedef unsigned int uint;/WN为复常数码

4、表,为提高计算速度/其值为exp(-j*2*pi/LENGTH)n/n为0.LENGTH/2-1JComplex WNLENGTH/2= 1.0000,-0.0000, 0.9997,-0.0245, 0.9988,-0.0491, 0.9973,-0.0736, 0.9952,-0.0980, 0.9925,-0.1224, 0.9892,-0.1467, 0.9853,-0.1710, 0.9808,-0.1951, 0.9757,-0.2191, 0.9700,-0.2430, 0.9638,-0.2667, 0.9569,-0.2903, 0.9495,-0.3137, 0.9415,

5、-0.3369, 0.9330,-0.3599, 0.9239,-0.3827, 0.9142,-0.4052, 0.9040,-0.4276, 0.8932,-0.4496, 0.8819,-0.4714, 0.8701,-0.4929, 0.8577,-0.5141, 0.8449,-0.5350, 0.8315,-0.5556, 0.8176,-0.5758, 0.8032,-0.5957, 0.7883,-0.6152, 0.7730,-0.6344, 0.7572,-0.6532, 0.7410,-0.6716, 0.7242,-0.6895, 0.7071,-0.7071, 0.6

6、895,-0.7242, 0.6716,-0.7410, 0.6532,-0.7572, 0.6344,-0.7730, 0.6152,-0.7883, 0.5957,-0.8032, 0.5758,-0.8176, 0.5556,-0.8315, 0.5350,-0.8449, 0.5141,-0.8577, 0.4929,-0.8701, 0.4714,-0.8819, 0.4496,-0.8932, 0.4276,-0.9040, 0.4052,-0.9142, 0.3827,-0.9239, 0.3599,-0.9330, 0.3369,-0.9415, 0.3137,-0.9495,

7、 0.2903,-0.9569, 0.2667,-0.9638, 0.2430,-0.9700, 0.2191,-0.9757, 0.1951,-0.9808, 0.1710,-0.9853, 0.1467,-0.9892, 0.1224,-0.9925, 0.0980,-0.9952, 0.0736,-0.9973, 0.0491,-0.9988, 0.0245,-0.9997, 0.0000,-1.0000, -0.0245,-0.9997,-0.0491,-0.9988,-0.0736,-0.9973,-0.0980,-0.9952,-0.1224,-0.9925, -0.1467,-0

8、.9892,-0.1710,-0.9853,-0.1951,-0.9808,-0.2191,-0.9757,-0.2430,-0.9700, -0.2667,-0.9638,-0.2903,-0.9569,-0.3137,-0.9495,-0.3369,-0.9415,-0.3599,-0.9330, -0.3827,-0.9239,-0.4052,-0.9142,-0.4276,-0.9040,-0.4496,-0.8932,-0.4714,-0.8819, -0.4929,-0.8701,-0.5141,-0.8577,-0.5350,-0.8449,-0.5556,-0.8315,-0.

9、5758,-0.8176, -0.5957,-0.8032,-0.6152,-0.7883,-0.6344,-0.7730,-0.6532,-0.7572,-0.6716,-0.7410, -0.6895,-0.7242,-0.7071,-0.7071,-0.7242,-0.6895,-0.7410,-0.6716,-0.7572,-0.6532, -0.7730,-0.6344,-0.7883,-0.6152,-0.8032,-0.5957,-0.8176,-0.5758,-0.8315,-0.5556, -0.8449,-0.5350,-0.8577,-0.5141,-0.8701,-0.

10、4929,-0.8819,-0.4714,-0.8932,-0.4496, -0.9040,-0.4276,-0.9142,-0.4052,-0.9239,-0.3827,-0.9330,-0.3599,-0.9415,-0.3369, -0.9495,-0.3137,-0.9569,-0.2903,-0.9638,-0.2667,-0.9700,-0.2430,-0.9757,-0.2191, -0.9808,-0.1951,-0.9853,-0.1710,-0.9892,-0.1467,-0.9925,-0.1224,-0.9952,-0.0980, -0.9973,-0.0736,-0.

11、9988,-0.0491,-0.9997,-0.0245;/复数乘法,进出参数均为结构体指针void JC_Mul(pJComplex pjc_1,pJComplex pjc_2,pJComplex pjc_r)pjc_r-real = (pjc_1-real)*(pjc_2-real) - (pjc_1-image)*(pjc_2-image);pjc_r-image = (pjc_1-real)*(pjc_2-image) + (pjc_1-image)*(pjc_2-real);/复数加法,进出参数均为结构体指针void JC_Add(pJComplex pjc_1,pJComplex

12、pjc_2,pJComplex pjc_r)pjc_r-real = (pjc_1-real) + (pjc_2-real);pjc_r-image = (pjc_1-image) + (pjc_2-image);/复数减法,进出参数均为结构体指针void JC_Sub(pJComplex pjc_1,pJComplex pjc_2,pJComplex pjc_r)pjc_r-real = (pjc_1-real) - (pjc_2-real);pjc_r-image = (pjc_1-image) - (pjc_2-image);/FFT主程序void popFFT(pJComplex da

13、taIN)uint i,j;uint mask_1,mask_2;uint tempIndex_1,tempIndex_2;uint t,m_1,m_2,n_1;JComplex bufLENGTH/2;for(i=0;iLENGTH;i+) mask_1 = 1; mask_2 = 1(LEVEL-1); tempIndex_1 = 0; for(j=0;jLEVEL;j+) if(i&mask_1) tempIndex_1 |= mask_2; mask_1 = 1; dataINi.real = dataINtempIndex_1.image; dataINtempIndex_1.image = 0;for(t=1;t=LEVEL;t+) m_1 = 11; n_1 = LENGTH/m_1;/蝶形区个数 for(i=0;in_1;i+) for(j=0;jm_2;j+) /前半部索引 tempIndex_1 = (it)+j; /总索引(前半部) tempIndex_2 = tempIndex_1 + m_2;/总索引(后半部) JC_Mul(&dataINtempIndex_2,&WNj*n_1,&bufj); JC_Sub(&dataINtempIndex_1,&bufj,&dataINtempIndex_2); /必须先算后半部 JC_Add(