快速傅立叶变换FFT算法DSP实验.docx
《快速傅立叶变换FFT算法DSP实验.docx》由会员分享，可在线阅读，更多相关《快速傅立叶变换FFT算法DSP实验.docx（29页珍藏版）》请在冰豆网上搜索。
快速傅立叶变换FFT算法DSP实验
快速傅立叶变换（FFT）算法实验
摘要：
FFT（FastFourierTransformation），即为快速傅里叶变换，是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
这种算法大大减少了变换中的运算量，使得其在数字信号处理中有了广泛的运用。
本实验主要要求掌握在CCS环境下用窗函数法设计FFT快速傅里叶的原理和方法；并且熟悉FFT快速傅里叶特性；以及通过本次试验了解各种窗函数对快速傅里叶特性的影响等。
引言：
快速傅里叶变换FFT是离散傅里叶变换DFT的一种快速算法。
起初DFT的计算在数字信号处理中就非常有用，但由于计算量太大，即使采用计算机也很难对问题进行实时处理，所以并没有得到真正的运用。
1965年J.W.库利和T.W.图基提出快速傅里叶变换，采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少,特别是被变换的抽样点数N越多,FFT算法计算量的节省就越显著。
从此，对快速傅里叶变换（FFT）算法的研究便不断深入，数字信号处理这门新兴学科也随FFT的出现和发展而迅速发展。
根据对序列分解与选取方法的不同而产生了FFT的多种算法，基本算法是基2DIT和基2DIF。
FFT的出现，使信号分析从时域分析向频域分析成为可能，极大地推动了信号分析在各领域的实际应用。
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：
其中，
，
设x（n）为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m），即N点DFT变换大约就需要N^2次运算。
当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）^2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。
这样变换以后，总的运算次数就变成N+2*（N/2）^2=N+N^2/2。
继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。
而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
8点DFT的FFT运算流图
计算离散傅里叶变换的快速方法，有按时间抽取的FFT算法和按频率抽取的FFT算法。
前者是将时域信号序列按偶奇分排，后者是将频域信号序列按偶奇分排。
它们都借助于的两个特点：
一是的周期性；另一是的对称性，这里符号*代表其共轭。
这样，便可以把离散傅里叶变换的计算分成若干步进行，计算效率大为提高。
时间抽取算法 令信号序列的长度为N=2,其中M是正整数，可以将时域信号序列x（n）分解成两部分，一是偶数部分x（2n），另一是奇数部分x（2n+1），其中。
于是信号序列x（n）的离散傅里叶变换可以用两个N/2抽样点的离散傅里叶变换来表示和计算。
一个抽样点数为N的信号序列x（n）的离散傅里叶变换，可以由两个N/2抽样点序列的离散傅里叶变换求出。
依此类推，这种按时间抽取算法是将输入信号序列分成越来越小的子序列进行离散傅里叶变换计算，最后合成为N点的离散傅里叶变换。
① N=2点的离散傅里叶变换的计算全由蝶形运算组成，需要M级运算，每级包括N/2个蝶形运算，总共有个蝶形运算。
所以，总的计算量为次复数乘法运算和Nlog2N次复数加法运算。
② FFT算法按级迭代进行，N抽样点的输入信号具有N个原始数据x0（n），经第一级运算后，得出新的N个数据x1（n），再经过第二级迭代运算，又得到另外N个数据x2（n），依此类推，直至最后的结果x（k）=xM（k）=X（k）在逐级迭代计算中，每个蝶形运算的输出数据存放在原来存贮输入数据的单元中，实行所谓“即位计算”，这样可以节省大量存放中间数据的寄存器。
③ 蝶形运算中加权系数随迭代级数成倍增加。
对于N=8,M=3情况，需进行三级迭代运算。
在第一级迭代中，只用到一种加权系数；蝶形运算的跨度间隔等于1。
在第二级迭代中，用到两种加权系数即、；蝶形运算的跨度间隔等于2。
在第三级迭代中，用到4种不同的加权系数即、、、；蝶形运算的跨度间隔等于4。
可见，每级迭代的不同加权系数的数目比前一级迭代增加一倍；跨度间隔也增大一倍。
④ 输入数据序列x（n）需重新排列为x（0）、x⑷、x⑵、x⑹、x⑴、x⑸、x⑶、x⑺，这是按照二进制数的码位倒置所得到的反序数，例如N=8中数“1”的二进制数为“001”，将其码位倒转变为“100”，即为十进制数“4”。
频率抽取算法 按频率抽取的FFT算法是将频域信号序列X（k）分解为奇偶两部分，但算法仍是由时域信号序列开始逐级运算，同样是把N点分成N/2点计算FFT，可以把直接计算离散傅里叶变换所需的N次乘法缩减到N/2次。
频率信号序列X（2l）是时间信号序列x1（n）+x2（n）的N/2点离散傅里叶变换，频率信号序列X（2l+1）是时间信号序列【x1（n）-x2（n）】的N/2点离散傅里叶变换，因此，N点离散傅里叶变换的计算，通过两次加（减）法和一次乘法，从原来序列获得两个子序列，所以，频率抽取算法也具有蝶形运算形式。
其计算量完全和时间抽取算法一样，即只需次乘法运算和Nlog2N次加（减）法运算。
实际上，频率抽取算法与时间抽取算法的信号流图之间存在着转置关系，如将流图适当变形，可以得出多种几何形状。
除了基2的FFT算法之外，还有基4、基8等高基数的FFT算法以及任意数为基数的FFT算法。
二、硬件框图：
“数字信号处理”实验室教学实验箱结构图
教学实验箱结构图
F2812-A评估板原理框图
F2812-A评估板实物图
DSP教学实验箱的硬件连接:
1．连接电源：
打开实验箱，取出三相电源连接线（如右图），将电源线的一端插入实验箱外部左恻箱壁上的电源插孔中。
确认实验箱面板上电源总开关（位于实验箱底板左上角）处于“关”的位置，连接电源线的另一端至220V交流供电插座上，保证稳固连接。
2．使用电源连接线（如右图，插头是带孔的）连接各模块电源：
确认实验箱总电源断开。
连接ICETEK-CTR板上边插座到实验箱底板上+12V电源插座；ICETEK-CTR板下边插座到实验箱底板上+5V电源插座；如使用PP（并口）型仿真器，则连接仿真器上插座到实验箱底板上+5V电源插座连接DSP评估板模块电源插座到实验箱底板上+5V电源插座。
注意各插
头要插到底，防止虚接或接触不良。
3．连接DSP评估板信号线：
当需要连接信号源输出到A/D输入插座
时，使用信号连接线（如右图）分别连接相应插座。
4、接通电源：
检查实验箱上220V电源插座（箱体左侧）中保险管是否完好，在连接电源线以后，检查各模块供电连线是否正确连接，打开实验箱上的电源总开关（位于实验箱底板左上角），使开关位于“开”的位置，电源开关右侧的指示灯亮。
三、软件流程图：
五、调试过程和步骤：
、软件调试
1．实验准备
-设置软件仿真模式。
-启动CCS。
2．打开工程，浏览程序，工程目录为D:
\dsp\t7\fft\fft.pjt
3．编译并下载程序
4．打开观察窗口：
选择菜单View->Graph->Time/Frequency…进行如下图所示设置。
选择菜单View->Graph->Time/Frequency…进行如下图所示设置。
选择菜单View->Graph->Time/Frequency…进行如下图所示设置。
5．清除显示：
在以上打开的窗口中单击鼠标右键，选择弹出式菜单中“ClearDisplay”功能。
6．设置断点：
在程序FFT.c中有注释“breakpoint”的语句上设置软件断点。
7.运行并观察结果
⑴选择“Debug”菜单的“Animate”项，或按F12键运行程序。
⑵观察“FFT”窗口中时域和频域图形。
注意：
由于实验运算复杂，需要等一会才能看到运行完结果。
8.退出CCS
、硬件连接
1．实验准备：
⑴连接实验设备。
⑵准备信号源进行AD输入。
①取出2根实验箱附带的信号线（如右图，两端均为单声道语音插头）。
②用1根信号线连接实验箱左侧信号源的波形输出A端口和“A/D输入”模块的“ADCIN0”插座注意插头要插牢、到底。
这样，信号源波形输出A的输出波形即可送到ICETEK-F2812A板的AD输入通道0。
③用1根信号线连接实验箱左侧信号源的波形输出B端口和“A/D输入”模块的“ADCIN1”插座注意插头要插牢、到底。
这样，信号源波形输出B的输出波形即可送到ICETEK-F2812A板的AD输入通道1。
④设置波形输出A：
-向内侧按波形频率选择旋钮，直到标有正弦波的指示灯点亮。
-上下调节波形频率选择旋钮，直到标有100-1KHz的指示灯点亮。
-调节幅值调整旋钮，将波形输出A的幅值调到适当位置。
⑤设置波形输出B：
-向内侧按波形频率选择旋钮，直到标有正弦波的指示灯点亮。
-上下调节波形频率选择旋钮，直到标有1K-10KHz的指示灯点亮。
-调节幅值调整旋钮，将波形输出B的幅值调到适当位置。
注意：
由于模数输入信号未经任何转换就进入DSP，所以必须保证输入的模拟信号的幅度在0-3V之间。
必须用示波器检测信号范围，保证最小值0V最大值3V，否则容易损坏DSP芯片的模数采集模块。
2．设置CodeComposerStudio2.21在硬件仿真（Emulator）方式下运行：
3．启动CodeComposerStudio2.21：
选择菜单Debug->ResetCPU。
4．打开工程文件：
工程目录：
D:
\dsp\t8\mixerfft\mixerfft.pjt
5．编译、下载程序。
6．运行程序观察结果：
按CTR控制板的K6键，可以显示A、B两信号源频谱，K7键实现混频显示，按K8实现键A、B两信号源分屏显示。
7．退出CCS：
六、实验结果：
七、结果分析：
通过观察软件仿真结果中FFT”窗口的时域图形和频域图形可知，由上述程序所计算出的测试波形的功率谱与CCS计算出的快速傅里叶变换的结果相近；通过硬件仿真结果可知，按CTR控制板的K6键，可以显示出A、B两信号源频谱，按K7键可以实现混频显示，按K8实现键A、B两信号源分屏显示。
八、结束语：
通过这次实验更加深刻地理解了快速傅里叶变换的物理意义，并且掌握了用窗函数法设计FFT快速傅里叶变换的原理和方法，进一步熟悉了FFT快速傅里叶变换特性，通过对快速傅里叶变换的实际应用了解了各种窗函数对快速傅里叶变换特性的影响。
本科学习期间，通过学习数字信号处理了解到快速傅里叶变换的知识，但当时仅限于对其理论层面的理解，通过这次实验，加强了自己动手能力的锻炼，使自己进一步认识到理论知识与实践结合起来的重要性，从而为自己以后的学习工作提供良好的发展方向。
《DSP原理与应用》这门课程在介绍DSPs芯片特点和应用的基础上，以TI公司C28x系列的TMS320F2812芯片为描述对象，系统的介绍了DSPs芯片的基本特点、硬件结构、工作原理、开发环境和使用方法，其中包括CPU内部结构、时钟和系统控制、存储空间及通用I/O接口、中断管理方式、片内外设、寻址方式和指令系统，以及本实验所用到的集成开发环境CCS、DSP最小系统和相应的软件设计。
通过学习《DSP原理与应用》我基本掌握了DSPs芯片的主要知识体系，并且通过上课老师的实物讲解和课本概念的联系，结合理论知识和实际应用，基本建立了DSP系统的基本概念与逻辑概念、物理概念之间的联系，同时通过这次实验进一步把DSPs的基本概念和原理应用到了实际的DSP系统中。
九、程序附录：
C语言程序代码：
程序1
#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;
}
}
程序2
#include"DSP281x_Device.h"//DSP281xHeaderfileIncludeFile
#include"DSP281x_Examples.h"//DSP281xExamplesIncludeFile
#include"f2812a.h"
#include"LCD.h"
#include"math.h"
#defineADCNUMBER300
//Prototypestatementsforfunctionsfoundwithinthisfile.
interruptvoidadc_isr（void）;
voidDelay（unsignedintnDelay）;
structstruLCDGraphstruGraph,struGraph1;
unsignedintnScreenBuffer[30*128];
//Globalvariablesusedinthisexample:
Uint16LoopCount;
Uint16ConversionCount;
Uint16Voltage1[1024];
Uint16Voltage2[1024];
Uint16Voltage_1,Voltage_2,flage=0;
Uint16nGraphBuf1[ADCNUMBER],nGraphBuf2[ADCNUMBER];
intnGraphBuf3[ADCNUMBER];
intci=0,keyflage,nAD;
Uint16nMixing[1024];
//液晶----------------------------------------------------------
#defineCTRLED（*（unsignedint*）0x108004）//port8004
#defineMCTRKEY（*（unsignedint*）0x108005）//port8005
#defineCTRCLKEY（*（unsignedint*）0x108006）//port8006
#defineCTRSTATUS（*（unsignedint*）0x108000）//port8000
//#definepi3.1415926
intnModeAD;
Uint16ad1,ad2;
#definePI3.1415926
#defineSAMPLENUMBER128
intINPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];
//intnGraphBuf1[ADCNUMBER],nGraphBuf2[ADCNUMBER];
floatfWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER];
floatsin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];
floatfDataR[SAMPLENUMBER],fDataI[SAMPLENUMBER],w[SAMPLENUMBER];
voidFFT（floatdataR[SAMPLENUMBER],floatdataI[SAMPLENUMBER]）;
voidInitForFFT（）;
//voidMakeWave（）;
intj,a=0;
main（）
{intj,uWork,uWork1;
unsignedint*pWork;
InitSysCtrl（）;//初始化cpu
//InitPll（0x5）;
DINT;//关中断
LCDTurnOff（）;
LCDSetScreenBuffer（nScreenBuffer）;
for（uWork=0,pWork=nScreenBuffer;uWork<30*128;uWork++,pWork++）（*pWork）=0;
LCDSetDelay（0）;
LCDTurnOn（）;//打开显示
LCDCLS（）;//清除显示内存
InitXintf（）;
InitPieCtrl（）;//初始化pie寄存器
/*for（j=0;j<1024;j++）
{
Voltage1[j]=0;
Voltage2[j]=0;
}*/
IER=0x0000;//禁止所有的中断
IFR=0x0000;
InitPieVectTable（）;//初始化pie中断向量表
//Interruptsthatareusedinthi
升级会员 