基于DSP的快速傅里叶FFT算法.docx
《基于DSP的快速傅里叶FFT算法.docx》由会员分享，可在线阅读，更多相关《基于DSP的快速傅里叶FFT算法.docx（19页珍藏版）》请在冰豆网上搜索。
基于DSP的快速傅里叶FFT算法
1.设计目的
1.1.设计目的
1．掌握用窗函数法设计FFT快速傅里叶的原理和方法；
2．熟悉FFT快速傅里叶特性；
3．了解各种窗函数对快速傅里叶特性的影响。
1.2.使用设备
PC兼容机一台，操作系统为Windows2000（或Windows98，WindowsXP，以下默认为Windows2000），安装CodeComposerStudio2.0软件。
2.设计任务与要求
按原程序仿真完成后，修改参数，观察波形变化。
3.原理与分析
1．FFT的原理和参数生成公式
公式
（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->Rebuiltall
File->LoadProgram->fft.out
4．打开观察窗口：
选择菜单View->Graph->Time/Frequency…进行如下图所示设置。
图2-25-1观察窗口设置1
图2-25-2观察窗口设置2
图2-25-3观察窗口设置3
5．清除显示：
在以上打开的窗口中单击鼠标右键，选择弹出式菜单中“ClearDisplay”功能。
6．设置断点：
在程序FFT.c中有注释“breakpoint”的语句上设置软件断点。
7．运行并观察结果
⑴选择“Debug”菜单的“Animate”项，或按F12键运行程序。
⑵观察“FFT”窗口中时域和频域图形。
8．退出CCS
5.软件设计
#include"DSP281x_Device.h"//DSP281xHeaderfileIncludeFile
#include"DSP281x_Examples.h"//DSP281xExamplesIncludeFile
#include"f2812a.h"
#include"math.h"
#definePI3.1415926
#defineSAMPLENUMBER128
voidInitForFFT（）;
voidMakeWave（）;
//voidFFT（floatdataR[SAMPLENUMBER],floatdataI[SAMPLENUMBER]）;
intINPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];
floatfWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];
floatsin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];
voidFFT（floatdataR[SAMPLENUMBER],floatdataI[SAMPLENUMBER]）
{
intx0,x1,x2,x3,x4,x5,x6,xx;
inti,j,k,b,p,L;
floatTR,TI,temp;
/**********followingcodeinvertsequence************/
for（i=0;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{
dataR[i]=dataI[i];dataI[i]=0;
}
/**************followingcodeFFT*******************/
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];
}/*ENDfor（3）*/
}/*ENDfor
（2）*/
}/*ENDfor
（1）*/
for（i=0;i{
w[i]=sqrt（dataR[i]*dataR[i]+dataI[i]*dataI[i]）;
}
}/*ENDFFT*/
main（）
{
inti;
InitForFFT（）;
MakeWave（）;
for（i=0;i{
fWaveR[i]=INPUT[i];
fWaveI[i]=0.0f;
w[i]=0.0f;
}
FFT（fWaveR,fWaveI）;
for（i=0;i{
DATA[i]=w[i];
}
while
（1）;//breakpoint
}
voidInitForFFT（）
{
inti;
for（i=0;i{
sin_tab[i]=sin（PI*2*i/SAMPLENUMBER）;
cos_tab[i]=cos（PI*2*i/SAMPLENUMBER）;
}
}
voidMakeWave（）
{
inti;
for（i=0;i{
INPUT[i]=sin（PI*2*i/SAMPLENUMBER*3）*1024;
}
}
6.系统仿真及调试
通过观察频域和时域图，程序计算出了测试波形的功率谱，与CCS计算的FFT结果相近。
7.完成结果或效果
1．当“#defineSAMPLENUMBER128”时
运行结果为：
2．当“#defineSAMPLENUMBER64”时
将程序修改为：
……
intx0,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．当“#defineSAMPLENUMBER32”时
将程序修改为：
……
intx0,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．当“#defineSAMPLENUMBER16”时
将程序修改为：
……
intx0,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．当“#defineSAMPLENUMBER256”时
将程序修改为：
……
intx0,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．当“#defineSAMPLENUMBER512”时
将程序修改为：
……
intx0,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．当“#defineSAMPLENUMBER620”时
将程序修改为：
……
intx0,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时，运算超出内存范围，不能编译。
8.心得体会
通过本次实验巩固了理论知识，加强了对本课程的理解，更加熟悉了DSP的操作及应用。
9.参考文献
（1）高海林，钱满义。
DSP技术及其应用。
清华大学出版社，北京交通大学出版社。
2009
（2）张雄伟，曹铁勇，陈亮，杨吉斌等。
DSP芯片的原理与开发应用（第4版）。
电子工业出版社，2009陈金鹰。
DSP技术及应用。
机械工业出版社，2004
（3）张彦编。
DSP技术及应用实验指导书.郑州轻工业学院计算机及通信工程学院。
升级会员 