本帖最后由 chaily 于 2014-3-30 13:35 编辑

之前有过不少的玩家发过许多关于FFT的贴了，但关于可随意移植到任意单片机的FFT算法个人感觉并不是很完善，尤其对于初学者和一些玩51的，大神看了不要笑话，嘿嘿。
这里对网上流传的其中一类FFT算法进行了个人测试和总结，在这里分享给大家。所有的资料都在附件里面了，这里贴出主要源代码和详细注释，里面也还有一些本人不清楚的问题，有待大神一起探讨。

/********************************************************************
这里讲讲傅立叶快速变换FFT算法的C例程，不恰当的地方还请探讨。
此例程可以移植到任何的单片机，但前提是单片机内存要够大(如果做32点的应该256字节内存就足够了，256字节内存的51机做64点好像也够)，至于速度快或慢，只要你能承受那些慢速的就无所谓。
这里以128点采样作为例子。要注意的是，采样点数N和采样率Fs还有信号频率F的关系。
FFT算法可以把对波形采样回来的数据进行转换，比如有一个模拟信号过来，假设对其以某个采样率采样了128个点的数据，经过FFT换算以后，就可以获得这段波形里面含有的不同频率的特征，也就是变成频谱。换算的结果是以复数的形式表现的，比如某个频率的信息以a＋bi的形式表示，这里a是实部，b是虚部。而例程中会以数组s16 Fft_Real[128]来表示128个频率的实部，以s16 Fft_Image[128]来表示128个频率的虚部。
采样点数N和采样率Fs还有信号频率F的关系是这样的，在一个模拟信号中，采样率至少应该是这个信号中需要分析的最高频率的2倍，而这个信号中需要分析的最低频率Fl＝Fs/N，也就是采样率除以采样点数，同时这个频率也是fft换算过程的频率分辩率。可能有点晕了，这里举例说明。
假设我们有一个信号里面含有100Hz的信号和含有1000Hz的信号叠加在一起，用示波器是很难分清的，通过对这信号进行采样和FFT转换就可以得到这组信号的频谱特征，因为需要分析的最低信号为100Hz，我们假设现在对信号要进行128点的采样，那么我们至少要能分辨100Hz的信号，采样率可以是100Hzx128点＝12800Hz，也就是我们需要以每秒采集12800点的速度来采集1/100秒的时间(因为只采集128点，所以时间上是1/100秒)，当然也可以使用6400Hz的采样率，这样分辩率是50Hz，这里就以12800Hz为例。
采集回来的数据要怎么处理呢？首先对于FFT算法，需要有数组数据输入，经过FFT换算以后会以数组数据输出。我们可以直接定义两个数组s16 Fft_Real[128]和s16 Fft_Image[128]，我们把采集来的128个数据按一个特殊的顺序输入到s16 Fft_Real[128]这个数组里面(当作实部输入)，而s16 Fft_Image[128]则全部输入0，然后运行fft算法转换。转换过后，s16 Fft_Real[128]里面的128个数会变成128个频率的实部，而s16 Fft_Image[128]里面的128个数会变成128个频率的虚部。哪128个频率呢，实际上只能获得64个频率。在这两个输出的数组里面，数组s16 Fft_Real[128]的第一个数(第一个数是s16 Fft_Real[0])和s16 Fft_Image[128]的第一个数(第一个数是s16 Fft_Image[0])分别是0Hz这个频率(直流分量)的复数形式结果里面的实部和虚部，而第二个数是100Hz的复数形式的结果，第三个数是200Hz的结果，依此类推，第64个数是6300Hz的结果。假如我们想知道信号中是否含有100Hz这个频率存在，那么我们需要查看s16 Fft_Real[1]和s16 Fft_Image[1]这两个数，把实部和虚部分别进行平方求和以后再开平方，就是100Hz信号的模，把模值除以采样点数N的1/2就得到100Hz的信号的峰值。同理，要知道信号里面是否含有1000Hz，我们需要取出s16 Fft_Real[10]和s16 Fft_Image[10]来计算。这128组数据里面，只有前面64组是有用的，后面的都是镜像数据，我们不理会。因为采样率是12800Hz，我们能查看的最高频率是6400Hz，也就是需要拿s16 Fft_Real[64]和s16 Fft_Image[64]出来算，但这两个值很可能不再准确，因为是极限的问题。

以下我们来看看相关的例程，这个例程可以用于移植到任何的单片机中，只要单片机的速度不太慢(太慢你愿意慢慢等也没关系)，内存足够就可以运行。这个例程可以修改为256点和512点或者1024点甚至更高点数的fft，下面会讲到需要修改一些什么参数。增加采样点数可以在采样率不变的情况下增加分辩率，如果点数增加一倍，采样率也提高一倍，那么分辩率不变，但能识别的最高频率将提高一倍，这是采样点数增加的好处，但同时会减慢运算速度。此例程使用stm32f103在72M主频下测试256点fft,运算时间为10mS左右，官方的函数库不知道会不会经过优化而更快，但考虑到许多人不喜欢使用函数库，也考虑这个例程的可移植性，这里直接使用fft函数而不讲函数库的内容。

*********************************************************************/



//以下为FFT傅立叶转换算法用到的相关定义
u8 ADC_Count=0;
/********注意，这里u8就是8位的unsigned char类型数据，在51里面就是unsigned char。以下的数组s16 Fft_Real[128]就是16位signed类型的数组，相当于signed int，而Fft_Real只是一个数组的名字，随便起也可以的，Fft_Image也一样**********/

u8 i,j,k,b,p; //这是用到的一些寄存器，注意这里如果要算256点以上的fft的话，需改成16位的u16。

s16 Temp_Real,Temp_Imag,temp; //中间临时变量，名称也是自己定义的，但要与fft函数里面的对应
u16 TEMP1; //用于求功率的，可不需要
s16 Fft_Real[128]; //fft实部，128数组
s16 Fft_Image[128]; //fft虚部，128数组

/*****以下为放大128倍后的sin正弦函数数组表格，这个表格可以用“正弦波表生成”这个软件来生成，要注意需要做多少点的fft，就需要生成多少点的。数据类型要选择整形，不要选择正整型，生成位数都为8位，这里为了能在普通单片机快速运行，不使用浮点数，使用查表法以达到最快的运算。如果是256点以上的，数组数据类型也是使用s8，但数组数量增加，那就要求更大的内存，注意这里相当于指针用法，如果内存足够的，最好把表格写入内存，那样运行速度快，如果内存资源紧的，就把表格写入程序区，对于51就是一个signed char code table[N]和signed char table[N]的区别，带了code就告诉编译器我这个表格要放在ROM程序存储区，不带code就直接放入内存，其他单片机不同写法自己研究研究  ******/
s8 SIN_TAB[128]={0x00, 0x06, 0x0c, 0x12, 0x18, 0x1e, 0x24, 0x2a,
0x30, 0x36, 0x3b, 0x41, 0x46, 0x4b, 0x50, 0x55,
0x59, 0x5e, 0x62, 0x66, 0x69, 0x6c, 0x70, 0x72,
0x75, 0x77, 0x79, 0x7b, 0x7c, 0x7d, 0x7e, 0x7e,
0x7f, 0x7e, 0x7e, 0x7d, 0x7c, 0x7b, 0x79, 0x77,
0x75, 0x72, 0x70, 0x6c, 0x69, 0x66, 0x62, 0x5e,
0x59, 0x55, 0x50, 0x4b, 0x46, 0x41, 0x3b, 0x36,
0x30, 0x2a, 0x24, 0x1e, 0x18, 0x12, 0x0c, 0x06,
0x00, -0x06, -0x0c, -0x12, -0x18, -0x1e, -0x24, -0x2a,
-0x30, -0x36, -0x3b, -0x41, -0x46, -0x4b, -0x50, -0x55,
-0x59, -0x5e, -0x62, -0x66, -0x69, -0x6c, -0x70, -0x72,
-0x75, -0x77, -0x79, -0x7b, -0x7c, -0x7d, -0x7e, -0x7e,
-0x7f, -0x7e, -0x7e, -0x7d, -0x7c, -0x7b, -0x79, -0x77,
-0x75, -0x72, -0x70, -0x6c, -0x69, -0x66, -0x62, -0x5e,
-0x59, -0x55, -0x50, -0x4b, -0x46, -0x41, -0x3b, -0x36,
-0x30, -0x2a, -0x24, -0x1e, -0x18, -0x12, -0x0c, -0x06};

//以下是放大128倍后的cos余弦函数数组表格，这里注意事项与上面相同，只不过选择余弦来生成
s8 COS_TAB[128]={0x7f, 0x7e, 0x7e, 0x7d, 0x7c, 0x7b, 0x79, 0x77,
0x75, 0x72, 0x70, 0x6c, 0x69, 0x66, 0x62, 0x5e,
0x59, 0x55, 0x50, 0x4b, 0x46, 0x41, 0x3b, 0x36,
0x30, 0x2a, 0x24, 0x1e, 0x18, 0x12, 0x0c, 0x06,
0x00, -0x06, -0x0c, -0x12, -0x18, -0x1e, -0x24, -0x2a,
-0x30, -0x36, -0x3b, -0x41, -0x46, -0x4b, -0x50, -0x55,
-0x59, -0x5e, -0x62, -0x66, -0x69, -0x6c, -0x70, -0x72,
-0x75, -0x77, -0x79, -0x7b, -0x7c, -0x7d, -0x7e, -0x7e,
-0x7f, -0x7e, -0x7e, -0x7d, -0x7c, -0x7b, -0x79, -0x77,
-0x75, -0x72, -0x70, -0x6c, -0x69, -0x66, -0x62, -0x5e,
-0x59, -0x55, -0x50, -0x4b, -0x46, -0x41, -0x3b, -0x36,
-0x30, -0x2a, -0x24, -0x1e, -0x18, -0x12, -0x0c, -0x06,
0x00, 0x06, 0x0c, 0x12, 0x18, 0x1e, 0x24, 0x2a,
0x30, 0x36, 0x3b, 0x41, 0x46, 0x4b, 0x50, 0x55,
0x59, 0x5e, 0x62, 0x66, 0x69, 0x6c, 0x70, 0x72,
0x75, 0x77, 0x79, 0x7b, 0x7c, 0x7d, 0x7e, 0x7e};

/***********以下是采样存储序列表，这个表格可以在FFT_Code_Tables.h这个文件里面找到，更多点的FFT也在里面找，就是里面的unsigned int code BRTable[N]的那些，是用来控制采样点数据输入到实部数组s16 Fft_Real[128]里面的***************/
u8 LIST_TAB[128]={0,64,32,96,16,80,48,112,
8,72,40,104,24,88,56,120,
4,68,36,100,20,84,52,116,
12,76,44,108,28,92,60,124,
2,66,34,98,18,82,50,114,
10,74,42,106,26,90,58,122,
6,70,38,102,22,86,54,118,
14,78,46,110,30,94,62,126,
1,65,33,97,17,81,49,113,
9,73,41,105,25,89,57,121,
5,69,37,101,21,85,53,117,
13,77,45,109,29,93,61,125,
3,67,35,99,19,83,51,115,
11,75,43,107,27,91,59,123,
7,71,39,103,23,87,55,119,
15,79,47,111,31,95,63,127};


//以下为fft算法主体部分
void FFT(void)
{
u8 N=7; //这里因为128是2的7次方，如果是计算256点，则是2的8次方，N就是8，如果是512点则N＝9，如此类推
u16 NUM_FFT=128; //这里要算多少点的fft就赋值多少，值只能是2的N次方
    for( i=1; i<=N; i++)             /* for(1) */
    {
        b=1;
        b <<=(i-1);                 //蝶式运算，用于计算 隔多少行计算。例如第一级 1和2行计算，，，第二级
        for( j=0; j<=b-1; j++)       /* for (2) */
        {
            p=1;
            p <<= (N-i);            
            p = p*j;
            for( k=j; k<NUM_FFT; k=k+2*b)   /* for (3) 基二fft */
            {
                Temp_Real = Fft_Real[k]; Temp_Imag = Fft_Image[k]; temp = Fft_Real[k+b];
                Fft_Real[k] = Fft_Real[k] + ((Fft_Real[k+b]*COS_TAB[p])>>7) + ((Fft_Image[k+b]*SIN_TAB[p])>>7);
                Fft_Image[k] = Fft_Image[k] - ((Fft_Real[k+b]*SIN_TAB[p])>>7) + ((Fft_Image[k+b]*COS_TAB[p])>>7);
                Fft_Real[k+b] = Temp_Real - ((Fft_Real[k+b]*COS_TAB[p])>>7) - ((Fft_Image[k+b]*SIN_TAB[p])>>7);
                Fft_Image[k+b] = Temp_Imag + ((temp*SIN_TAB[p])>>7) - ((Fft_Image[k+b]*COS_TAB[p])>>7);     
                //移位，防止溢出。结果已经是本值的1/64               
              Fft_Real[k] >>= 1;            
                Fft_Image[k] >>= 1;
               Fft_Real[k+b]  >>= 1;                 
                Fft_Image[k+b]  >>= 1;
            }     
        }
    }
//         Fft_Real[0]=Fft_Image[0]=0;  //去掉直流分量，也可以不去掉
//   Fft_Real[63]=Fft_Image[63]=0;
/********以上已经把128点的实部和虚部求完，下一次运算前需要把所有虚部重新清零
要求某个频率点的模，则模值=根号(实部平方+虚部平方)，即sqrt((Fft_Real[n]*Fft_Real[n])+(Fft_Image[n]*Fft_Image[n]))
第n个频率点的值是数组上的Fft_Real[n]和Fft_Image[n]，而Fft_Real[0]是直流分量。Fft_Real[1]是最低频率点，也是最小频率分辩率值，等于采样率/采样点数N
波形峰值大小=模值/(N/2)， N为采样点数
注意，由于上面把求得的值已经移位除法，相当于缩小了64倍(移位7位好像是128倍吧？？后面为什么还要移动一位？这里待高手指点，本人也不是很清楚，这里只做移植总结)
得到的那些实部虚部的结果爱怎么处理怎么处理，可以做音频频谱强度显示啦什么的******************/
}




void Fft_Imagclear(void) //fft虚部清零函数，在运行FFT函数之前需要先运行这个
{
        u8 a;  //注意这里如果是256点以上要改成u16，下面的a<128条件也要相应的修改
        for(a=0;a<128;a++)
                {
                        Fft_Image[a]=0;
                }
}



/********采样程序如下，这是stm32里面的采样函数而已，别的单片机该怎么写自己搞定，反正是采集够点数了，就停止采样，采样的数据输入到实部的时候是按那个LIST_TAB[128]表格顺序的，这里函数已经写好一个模型了，也就是Fft_Real[LIST_TAB[ADC_Count]]=你的采样数据，注意只输入到实部，虚部全部输入为0，所以运行FFT前虚部需要清零******/
void TIM2_IRQHandler(void)   //时间中断服务函数， 这里控制自动定时采样
{
        TIM_ClearFlag(TIM2, TIM_FLAG_Update); //该服务函数定时时间的倒数就是采样率
        //stm32里面所有中断都不会自动清零，需要返回前清零标志位，否则会一直无限产生中断，其他单片机自己考虑实际情况

        Fft_Real[LIST_TAB[ADC_Count]]=ADC_GetConversionValue(ADC1); //这里ADC_GetConversionValue(ADC1)是获得ADC的采样值，这里不同的单片机可能处理不同
                        if(ADC_Count<127) //采样计数条件判断
                                        {ADC_Count++;}
                        else
                                        {TIM_Cmd(TIM2,DISABLE);ADC_Count=128;} //采样够了就把控制中断的定时器关闭，这里写法可以自己发挥
               
}





关于主函数可以大概这样写：
int main(void)
{
        while(1)
           {
                ADC_Count=0;
                TIM_Cmd(TIM2,ENABLE); //打开控制采样的定时器中断的定时器，至于你的ADC是怎么采样的自己搞定
                while(ADC_Count<128); //当采样点达到128点前循环等待，达到128点后跳过，这个点数随着需要的点数修改，需要配合采样时间中断的写法
                FFT(); //运行FFT换算
                Fft_Imagclear(); //运行虚部清零
          }

}



附件里面的内容有：
FFT_Code_Tables.h 文件
FFT变换的实际意义
正弦波表生成软件