基于DSP的C程序实验报告快速傅立叶变换FFT算法

1、 目录1 引言2 实验原理3 FFT基本结构 （1）信号流图 （2）软件程序流图4 实验程序5 调试过程与步骤6 实验结果7 结果分析8 遇到的问题及解决办法9 实验体会 实验题目：快速傅立叶变换(FFT)算法1 引言众所周知，FFT 是离散傅立叶变换（DFT）的一种快速算法。由于计算DFT 时一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法则需二次实数加法。每运算一个X（k）需要4N 次复数乘法及2N+2（N-1）=2（2N-1）次实数加法。所以整个DFT运算总共需要4N2 次实数乘法和N*2(2N-1)=2N(2N-1)次实数加法。如此一来，计算时乘法次数和加法次数都是和N2 成正

2、比的，当N 很大时，运算量是可观的，因而需要改进对DFT 的算法提高运算速度。我们观察DFT的系数特性，可发现其对称性和周期性，我们可以将DFT 运算中有些项合并。因此，FFT就孕育而生了。2 实验原理1FFT 的原理和参数生成公式：我们先设序列长度为N=2L，L 为整数。将N=2L 的序列x(n)(n=0,1,，N-1)，按N的奇偶分成两组，也就是说我们将一个N 点的DFT 分解成两个N/2 点的DFT，他们又重新组合成一个如下式所表达的N 点DFT：我们称这样的RFFT 优化算法是包装算法：首先2N 点实数的连续输入称为“进包”。其次N 点的FFT 被连续运行。最后作为结果产生的N 点的合

3、成输出是“打开”成为最初的与DFT 相符合的2N 点输入。三FFT 的基本结构：1.FFT信号流图如下： 整个过程共有log2N次，每次分组间隔为2（L-1）-1=<L<=log2N（1）如上图第一次蝶形运算间隔为一，如第一个和第二个，第三个和第四个，以此类推；第二次间隔为二，如第一个和第三个，第二个和第四个等（2）基本运算单元以下面的蝶形运算为主：计算公式如下： （3）在FFT运算中，旋转因子WmN=cos(2m/N)-jsin(2m/N)，求正弦和余弦函数值的计算量是很大的。（4）本程序采用的输入信号为：1024*sin(2*pi*3*t)，采样频率为1024 2程序流程图：四

4、实验程序#include "DSP281x_Device.h"     / DSP281x Headerfile Include File#include "DSP281x_Examples.h"   / DSP281x Examples Include File#include "f2812a.h"#include"math.h&quo

5、t;#define PI 3.1415926#define SAMPLENUMBER 128void InitForFFT();void MakeWave();/oid FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER);int INPUTSAMPLENUMBER,DATASAMPLENUMBER;float fWaveRSAMPLENUMBER,fWaveISAMPLENUMBER,wSAMPLENUMBER;float 

6、;sin_tabSAMPLENUMBER,cos_tabSAMPLENUMBER;void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER)int x0,x1,x2,x3,x4,x5,x6,xx;int i,j,k,b,p,L;float TR,TI,temp;/* following code invert sequence */for ( i=0;i<SAMPLENUMBER;i+ )x0=

7、x1=x2=x3=x4=x5=x6=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;dataIxx=dataRi;for ( i=0;i<SAMPLENUMBER;i+ )dataRi=dataIi; data

8、Ii=0; /* following code FFT */for ( L=1;L<=7;L+ ) /* for(1) */b=1; i=L-1;while ( i>0 ) b=b*2; i-; /* b= 2(L-1) */for ( j=0;j<=b-1;j+ ) /* for (2) */p=1

9、; i=7-L;while ( i>0 ) /* p=pow(2,7-L)*j; */p=p*2; i-;p=p*j;for ( k=j;k<128;k=k+2*b ) /* for (3) */TR=dataRk; TI=dataIk; temp=dataRk+b;dataRk=dataRk+dataRk+b*cos_tabp+dataIk+b*sin_tabp;dataIk=dataIk-dataRk+b*si

10、n_tabp+dataIk+b*cos_tabp;dataRk+b=TR-dataRk+b*cos_tabp-dataIk+b*sin_tabp;dataIk+b=TI+temp*sin_tabp-dataIk+b*cos_tabp; /* END for (3) */ /* END for (2) */ /* END for (1) */for ( i=0;i<SAMPLENUMBER/2;i+ ) 

11、;wi=sqrt(dataRi*dataRi+dataIi*dataIi); /* END FFT */main()int i;InitForFFT();MakeWave();for ( i=0;i<SAMPLENUMBER;i+ )fWaveRi=INPUTi;fWaveIi=0.0f;wi=0.0f;FFT(fWaveR,fWaveI);for ( i=0;i<SAMPLENUMBER;i+ )DATAi=wi;while ( 1 ); /

12、 break pointvoid InitForFFT()int i;for ( i=0;i<SAMPLENUMBER;i+ )sin_tabi=sin(PI*2*i/SAMPLENUMBER);cos_tabi=cos(PI*2*i/SAMPLENUMBER);void MakeWave()int i;for ( i=0;i<SAMPLENUMBER;i+ )INPUTi=sin(PI*2*i/SAMPLENUMBER*3)*1024;五调试过程与步骤1编译并下

13、载程序。2打开观察窗口：*选择菜单View->Graph->Time/Frequency进行如下图所示设置。图1图2图33清除显示：在以上打开的窗口中单击鼠标右键，选择弹出式菜单中“Clear Display”功能。4设置断点：在程序FFT.c 中有注释“break point”的语句上设置软件断点。图45运行并观察结果。 选择“Debug”菜单的“Animate”项，或按Alt+F5 键运行程序。 观察“Test Wave”窗口中时域图形； 图5 在“Test Wave”窗口中点击右键，选择属性，更改图形显示为FFT。观察频域图形。 图6 观察“FFT”窗口中的由CCS 计算出的

14、正弦波的FFT。 图7六. 实验结果通过观察频域和时域图，程序计算出了测试波形的功率谱，与CCS 计算的FFT 结果相近。七. 结果分析(1)观察图6和图7，可以看到二者波形相似，但横纵坐标均不相同，纵坐标大约是二倍的关系，横坐标大约为142倍。 (2)观察图8，因为两个频率比较相近，因此出现了前两个频谱交叠的现象。八.遇到的问题及解决办法 1.CCS启动失败 经过检查，并没有发现什么问题，一时间无从下手，将USB插头换了一个插口后问题解决，可能原因是电脑的USB插口老化接触不良造成的。 2. 程序编译过程中产生了错误 经过检查，工程中缺少（.cmd）文件，添加后编译成功。 3. 未能出现想要的图像 经过仔细检查，浏览了程序后，发现断点设置有些问题，改正后问题解决，显示正确图像。 九实验体会 1.本次试验详细的介绍了FFT的特点和原理，以及设计方法，通过实验充分展示了FFT快速傅里叶变换相比于普通的算法的明显优点。2. 通过本次试验，我想我们已经熟练的掌握了FFT快速傅里叶变换的