基于DSP的快速傅里叶(FFT)算法

哈尔滨商业大学DSP课程设计报告题目快速傅立叶变换（FFT）算法专业电子信息工程班级 08级02班姓名学号王玉辉200810930172李砚秋200810930062杨兴臻200810930292尤琳200810930052指导教师姜海涛日期2011年12月12日目录1.设计目的 .... 错误！未定义书签。

1.1. 设计目的..................................................................... 错误！未定义书签。

1.2. 使用设备..................................................................... 错误！未定义书签。

2.设计任务与要求错误！未定义书签。

3.原理与分析 .. 错误！未定义书签。

4.实验步骤 .... 错误！未定义书签。

5.软件设计 .... 错误！未定义书签。

6.系统仿真及调试错误！未定义书签。

7.完成结果或效果错误！未定义书签。

8.心得体会 .... 错误！未定义书签。

9.参考文献 .... 错误！未定义书签。

1. 设计目的1.1. 设计目的1．掌握用窗函数法设计FFT 快速傅里叶的原理和方法； 2．熟悉FFT 快速傅里叶特性；3．了解各种窗函数对快速傅里叶特性的影响。

1.2. 使用设备PC 兼容机一台，操作系统为Windows2000(或Windows98，WindowsXP ，以下默认为Windows2000)，安装Code Composer Studio 2.0 软件。

2. 设计任务与要求按原程序仿真完成后，修改参数，观察波形变化。

3. 原理与分析1． FFT 的原理和参数生成公式)()()()()(2121222121k X W k X W r x W W r x k x k N rk N N r kN rk N Nr +=+=∑∑-=-=公式（1）FFT 运算公式FFT 并不是一种新的变换，它是离散傅立叶变换（DFT ）的一种快速算法。

由于我们在计算DFT 时一次复数乘法需用四次实数乘法和二次实数加法；一次复数加法则需二次实数加法。

每运算一个X （k ）需要4N 次复数乘法及2N+2（N-1）=2（2N-1）次实数加法。

所以整个DFT 运算总共需要4N^2 次实数乘法和N*2(2N-1)=2N(2N-1)次实数加法。

如此一来，计算时乘法次数和加法次数都是和N^2 成正比的，当N 很大时，运算量是可观的，因而需要改进对DFT 的算法减少运算速度。

根据傅立叶变换的对称性和周期性，我们可以将DFT 运算中有些项合并。

我们先设序列长度为N=2^L ，L 为整数。

将N=2^L 的序列x(n)(n=0,1,……，N-1)，按N 的奇偶分成两组，也就是说我们将一个N 点的DFT 分解成两个N/2 点的DFT ，他们又重新组合成一个如下式所表达的N 点DFT ：一般来说，输入被假定为连续的。

当输入为纯粹的实数的时候，我们就可以利用左右对称的特性更好的计算DFT 。

我们称这样的RFFT 优化算法是包装算法：首先2N 点实数的连续输入称为“进包”。

其次N 点的FFT 被连续被运行。

最后作为结果产生的N 点的合成输出是“打开”成为最初的与DFT 相符合的2N 点输入。

使用这战略，我们可以划分FFT 的大小，它有一半花费在包装输入O （N ）的操作和打开输出上。

这样的RFFT 算法和一般的FFT 算法同样迅速，计算速度几乎都达到了两次DFT 的连续输入。

下列一部分将描述更多的在TMS320C54x 上算法和运行的细节。

2． 程序流程图4. 实验步骤1．实验准备-设置软件仿真模式。

-启动CCS 。

2．打开工程，浏览程序，工程目录为C:\ICETEK-F2812-AG-EDUlab\DSP281x_examples\lab0503-FFT 3．编译并下载程序 Project->Rebuilt allFile->Load Program->fft.out 4．打开观察窗口：选择菜单View->Graph->Time/Frequency …进行如下图所示设置。

图2-25-1 观察窗口设置1图2-25-2 观察窗口设置2图2-25-3 观察窗口设置35．清除显示：在以上打开的窗口中单击鼠标右键，选择弹出式菜单中“Clear Display”功能。

6．设置断点：在程序FFT.c 中有注释“break point ”的语句上设置软件断点。

7．运行并观察结果⑴选择“Debug”菜单的“Animate ”项，或按F12 键运行程序。

⑵观察“FFT”窗口中时域和频域图形。

8．退出CCS5.软件设计#include "DSP281x_Device.h" // DSP281x Headerfile Include File#include "DSP281x_Examples.h" // DSP281x Examples Include File#include "f2812a.h"#include"math.h"#define PI 3.1415926#define SAMPLENUMBER 128void InitForFFT();void MakeWave();//void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]);int INPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];float fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]){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=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;dataI[xx]=dataR[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=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; 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=dataR[k]; TI=dataI[k]; temp=dataR[k+b];dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p];} /* END for (3) */} /* END for (2) */} /* END for (1) */for ( i=0;i<SAMPLENUMBER/2;i++ ){w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);}} /* END FFT */main(){int i;InitForFFT();MakeWave();for ( i=0;i<SAMPLENUMBER;i++ ){fWaveR[i]=INPUT[i];fWaveI[i]=0.0f;w[i]=0.0f;}FFT(fWaveR,fWaveI);for ( i=0;i<SAMPLENUMBER;i++ ){DATA[i]=w[i];}while ( 1 ); // break point}void InitForFFT(){int i;for ( i=0;i<SAMPLENUMBER;i++ ){sin_tab[i]=sin(PI*2*i/SAMPLENUMBER);cos_tab[i]=cos(PI*2*i/SAMPLENUMBER);}}void MakeWave(){int i;for ( i=0;i<SAMPLENUMBER;i++ ){INPUT[i]=sin(PI*2*i/SAMPLENUMBER*3)*1024;}}6.系统仿真及调试通过观察频域和时域图，程序计算出了测试波形的功率谱，与CCS 计算的FFT 结果相近。

7.完成结果或效果1．当“#define SAMPLENUMBER 128”时运行结果为：2．当“#define SAMPLENUMBER 64”时将程序修改为：……int x0,x1,x2,x3,x4,x5,xx;……x0=x1=x2=x3=x4=x5=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01;x3=(i/8)&0x01; x4=(i/16)&0x01; x5=(i/32)&0x01; xx=x0*32+x1*16+x2*8+x3*4+x4*2+x5;……for ( k=j;k<64;k=k+2*b )……修改设置运行结果为：3．当“#define SAMPLENUMBER 32”时将程序修改为：……int x0,x1,x2,x3,x4,xx;……x0=x1=x2=x3=x4=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01; x4=(i/16)&0x01;xx=x0*16+x1*8+x2*4+x3*2+x4;……for ( k=j;k<32;k=k+2*b )……修改设置运行结果为：4．当“#define SAMPLENUMBER 16”时将程序修改为：……int x0,x1,x2,x3,xx;……x0=x1=x2=x3 =0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;xx=x0*8+x1*4+x2*2+x3;……for ( k=j;k<16;k=k+2*b )……修改设置运行结果为：5．当“#define SAMPLENUMBER 256”时将程序修改为：……int x0,x1,x2,x3,x4,x5,x6,x7,xx;……x0=x1=x2=x3=x4=x5=x6=x7=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; x7=(i/128)&0x01;xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7;……for ( k=j;k<256;k=k+2*b )……运行结果为：6．当“#define SAMPLENUMBER 512”时将程序修改为：……int x0,x1,x2,x3,x4,x5,x6,x7,x8,xx;……x0=x1=x2=x3=x4=x5=x6=x7=x8=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; x7=(i/128)&0x01; x8=(i/256)&0x01;xx=x0*256+x1*128+x2*64+x3*32+x4*16+x5*8+x6*4+x7*2+x8; ……for ( k=j;k<512;k=k+2*b )……修改设置运行结果为：7．当“#define SAMPLENUMBER 620”时将程序修改为：……int x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,xx;……x0=x1=x2=x3=x4=x5=x6=x7=x8=x9=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; x7=(i/128)&0x01; x8=(i/256)&0x01; x9=(i/512)&0x01;xx=x0*512+x1*256+x2*128+x3*64+x4*32+x5*16+x6*8+x7*4+x8*2+x9;……for ( k=j;k<620;k=k+2*b )……修改设置运行结果：结论：经试验当SAMPLENUMBER超过650时，运算超出内存范围，不能编译。