[-玩转FFT算法-]-FFT算法单片机任你移植-个人总结篇

2019-12-09 20:00发布

本帖最后由 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变换的实际意义
正弦波表生成软件
友情提示: 此问题已得到解决,问题已经关闭,关闭后问题禁止继续编辑,回答。
99条回答
zhengensong
2019-12-16 12:58
本帖最后由 zhengensong 于 2014-4-4 08:20 编辑
chaily 发表于 2014-4-3 12:00
我也不清楚,所以这部分也是我的疑问,哈~


楼主看看余弦的表格中的第一个值COS(0)本应该等于1的,但是他取的是127,8位的有符号整数是从-128到+127的,这样他后边在计算完成之后只要用到COS和SIN的地方都要除以128的,这样也是为了防止溢出的,如果COS和SIN的值,最大的是1,如果取浮点数据,运算特别慢不说,占内存,结果也不精确,他这样扩大倍数计算,计算完成之后再把倍数降下来,这样就使得数据在整数范围内有较快的准确的运算了,后边的移位,也就是又除以了2,我猜测可能是怕得到的结果溢出,他在这里除以2,得到的相当于实部和虚部都除以2,这样就相当于实部和虚部的平方和开放除以2,如果计算模值,他只需要再除以32就可以了,不过他的表格中的数据128个他就用了64个我真心猜不到为啥,如果楼主能验证的话,还请楼主帮忙验证一下能否把正弦和余弦两个数组中的64-127删除,看看结果是否正确,我能想到的就这么多,希望对大家有帮助,不足之处希望大家再次纠正……

一周热门 更多>