用单片机做fft.txt

网上搜集的FFT算法C程序集

这次的任务是要用单片机做FFT(快速离散傅利叶变换),我觉得单片机越来越无敌了... 

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

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

贴出主要部分代码:

/************************************************************
//FFT(快速离散傅利叶变换)
//本例为8级256点数据
//HUST STI OE0402 POPAPPLE
//[e-mail]
//2007-8-9 18:10:20
************************************************************/

#include "stdio.h"

#define LEVEL 8 //级数
#define LENGTH 256 //点数

//用结构体来表示复数
typedef struct _JComplex
{
float real;
float image;
}JComplex,*pJComplex;

typedef unsigned int uint;


//WN为复常数码表,为提高计算速度
//其值为exp(-j*2*pi/LENGTH)^n
//n为[0..LENGTH/2-1]
JComplex WN[LENGTH/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,-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.6895,-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},{ 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.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.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.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.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 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 dataIN)
{
uint i,j;
uint mask_1,mask_2;
uint tempIndex_1,tempIndex_2;
uint t,m_1,m_2,n_1;
JComplex buf[LENGTH/2];

for(i=0;i<LENGTH;i++)
{
    mask_1 = 1;
    mask_2 = 1<<(LEVEL-1);
    tempIndex_1 = 0;
    for(j=0;j<LEVEL;j++)
    {
      if(i&mask_1)
        tempIndex_1 |= mask_2;
      mask_1 <<= 1;
      mask_2 >>= 1;
    }
    dataIN[i].real = dataIN[tempIndex_1].image;
    dataIN[tempIndex_1].image = 0;
}
for(t=1;t<=LEVEL;t++)
{
    m_1 = 1<<t;//蝶形区内数据个数
    m_2 = m_1>>1;
    n_1 = LENGTH/m_1;//蝶形区个数
    for(i=0;i<n_1;i++)
    {
      for(j=0;j<m_2;j++) //前半部索引
      {
        tempIndex_1 = (i<<t)+j; //总索引(前半部)
        tempIndex_2 = tempIndex_1 + m_2;//总索引(后半部)
        JC_Mul(&dataIN[tempIndex_2],&WN[j*n_1],&buf[j]);
        JC_Sub(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_2]); //必须先算后半部
        JC_Add(&dataIN[tempIndex_1],&buf[j],&dataIN[tempIndex_1]); //再算前半部
      }
    }
}
return;
}

//主程序,测试
int main()
{
int i;
JComplex dataIN[LENGTH];
for(i=0;i<LENGTH;i++)
    dataIN[i].image = i;
popFFT(dataIN);
for(i=0;i<LENGTH;i++)
    printf("%f\t%f\n",dataIN[i].real,dataIN[i].image);
}