实验八-快速傅立叶变换(FFT)实验

1实验七快速傅立叶变换（FFT）实验一实验目的1．熟悉CCS集成开发环境；2．了解FFT的算法原理和基本性质；3．熟悉DSP中cmd文件的作用及对它的修改；4．学习用FFT对连续信号和时域信号进行频谱分析的方法；5．利用DSPLIB中现有的库函数；6．了解DSP处理FFT算法的特殊寻址方式；7．熟悉对FFT的调试方法。二实验内容本实验要求使用FFT变换对一个时域信号进行频谱分析，同时进行IFFT。这里用到时域信号可以是来源于信号发生器输入到CODEC输入端，也可以是通过其他工具计算获取的数据表。本实验使用Matlab语言实现对FFT算法的仿真，然后将结果和DSP分析的结果进行比较，其中原始数据也直接来自Matlab。三实验原理一个N点序列][kx的DFT][mX,以及IDFT分别定义为：1,,1,0,][][10NmWkxmXkmNNk1,,1,0,][1][10NkWmXNkxkmNNm如果利用上式直接计算DFT,对于每一个固定的m,需要计算N次复数乘法，N-1次加法,对于N个不同的m,共需计算N的2次方复数乘法,N*(N-1)次复数加法.显然,随着N的增加,运算量将急剧增加,快速傅里叶算法有效提高计算速度（本例使用基2FFT快速算法）,利用FFT算法只需(N/2)logN次运算。四知识要点.1、CMD文件的功能及编写2、一种特殊的寻址方式：间接寻址间接寻址是按照存放在某个辅助寄存器的16位地址寻址的。C54x的8个辅助寄存器（AR0—AR7）都可以用来寻址64K字数据存储空间中的任何一个存储单元。3、TMS320C54xDSPLIB中关于FFT变换的一些函数的调用（SPRA480B.pdf）利用DSPLIB库时，在主程序中要包含头文件：54xdsp.lib4、FFT在CCS集成开发环境下的相关头文件#includetype.h//定义数据类型的头文件#includemath.h//数学函数的头文件，如sqrt.#includetms320.h//定义数据类型的头文件#includedsplib.h//DSPLIB库文件2五实验程序说明1、实验主要函数/***************************正变换*************************************/cbrev(x,x,NX/2);//倒序rfft(x,64,0);//实数FFT变换//求频谱由于FFT程序计算得到的数据只是频谱的实部和虚部，不包含计算幅度谱的//成分（所以描述DSP的参数中给出计算N点FFT的时间，是指不含计算幅度谱的时间），//因此要得到幅度频谱，必须另外增加程序语句来实现。for(i=0;iNX;i=i+2){p=x[i];//实部q=x[i+1];//虚部n=p*p+q*q;//实部平方加虚部平方f[m]=sqrt(n);//开方后存到f数组中m++;}/**************************反变换************************************/unpacki(x,NX);//还原cbrev(x,x,NX/2);//倒序rifft(x,NX,1);//实数IFFT变换2、各个函数的说明（详情见SPRA480B.pdf）（1）voidcbrev（DATA*x,DATA*r,unshortn）功能：为了FFT/IFFT得到一个正确顺序的变换结果，对他们的输入数据进行倒序。入口参数：x[2*n]x是一个2*n项的一维数组，数组中数据定义为短整型（16位有符号整型）。数组x是作为输入数据，函数对他的数据进行倒序。r[2*n]r是一个2*n项的一维数组，数组中数据定义为短整型（16位有符号整形）。数组r是作为输出数据，函数对x倒序后的结果存到r中。n定义为数组中复数的个数（两个实数表示一个复数），即为数组大小的1/2。函数的使用：函数是对复数进行倒序的，即把数组x中的数据认为是复数。有两个相邻的实数表示一个复数，偶地址为复数的实部，奇地址为复数的虚部。如下式，函数对X[0]+j*X[1]，X[2]+j*X[3]，………X[2n]+j*X[2n+1]…………X[2*N-2]+j*X[2*N-1]这些数据进行倒序。倒序后的结果也是按复数的实部、虚部依次存到r数组中的。3注意：数组中的元素个数必须为偶数。倒序时采用间接寻址，所以数组的首地址的末log（n）+1必须为0。（2）voidcfft(x,n,scale)原理及源程序说明：功能：对复数进行FFT变换。各项参数：x[2*n]x是一个2*n项的一维数组，数组中数据定义为短整形（16位有符号整形）。数组x既作为输入数据，又存放变换后的输出数据。n定义为数组中复数的个数（两个实数表示一个复数），，即为数组大小的1/2。Scale变换系数，如果为0，变换后结果乘以1/nx；否则结果乘以1。函数的使用：函数cfft(x,n,scale)是经过以下俩个宏定义而来的：#definedummy(x,n,scale)cfft##n(x,scale)#definecfft(x,n,scale)dummy(x,n,scale)原始函数为cfft##n(x,scale)，n可取值为16，32，64，128，256，512，1024。函数Cfft（）要求输入数据为倒序，即经过cbrev（）处理之后的数据。同cbrev（）一样，cfft（）也是对X[0]+j*X[1]，X[2]+j*X[3]，………X[2n]+j*X[2n+1]…………X[2*N-2]+j*X[2*N-1]进行的FFT变换，结果按实部/虚部存放。注意：数组中的元素个数必须为偶数。数组的首地址的末log（n）+1必须为0。4（3）rfft(x,n,scale)；实数FFT变换：原理及源程序说明：实数FFT变换涉及到下面两个算法。（1）、利用Ｎ点复序列的FFT算法计算两个Ｎ点实序列FFT设x[k]和y[k]都是N点实序列，X[m]和Y[m]分别表示他们对应的N点DFT。设h[k]=x[k]+jy[k]，已知求得h[k]的DFT为H[k]，根据DFT的性质可得X[m]=1/2{H[m]+H*[(-m)N]}Y[m]=1/(2j){H[m]-H*[(-m)N]}（2）、利用N点复序列计算2N点是序列FFT。（如rfft）设y[k]是一个长度为2N的实序列，Ｙ[M]是DFT。定义如下两个数组x[k]=y[2k]，h[k]=y[2k+1]，由上一个算法可以得到Y[m]=Ｘ[m]+W2NmH[m]Y[m+N]=Ｘ[m]-W2NmH[m]功能：对实数进行FFT变换。各项参数：x[n]x是一个n（n必须为偶数）项的一维数组，数组中数据定义为短整型（16位有符号整型）。数组x既作为输入数据，又存放变换后的输出数据。n定义为数组中实数的个数，即等于数组大小。Scale变换系数，如果为0，变换后结果乘以1/nx；否则结果乘以1。函数的使用：实数fft变换实际上也是把数组中的数据当成n/2个复数进行cfft变换，之后再调用一个调整的函数unpack()。所以可以把rfft(x,n,scale)看成cfft（x，n/2，scale）+unpack（x，n）。其它与cfft（）一样。N点实序列的频谱是N点复序列，需要2N个存储空间（实部、虚部分别占相邻两个存储空间），但由于实序列的频谱存在共轭对称的关系，已知前N/2点复序列，就可以通过共轭对称性求的后N/2点复序列，因此只要求N个存储空间（存放前N/2个复序列）就可以存放频谱。（4）cifft(x,n,scale)/rifft(x,2*n,scale)；复数iFFT/riFFT变换：功能：对复数进行IFFT变换。各项参数：x[2*n]x是一个2*n项的一维数组，数组中数据定义为短整型（16位有符号整形）。数组x既作为输入数据，又存放变换后的输出数据。5n定义为数组中复数的个数（两个实数表示一个复数），，即为数组大小的1/2。Scale变换系数，如果为0，变换后结果乘以1/nx；否则结果乘以1。函数的使用：函数cifft(x,n,scale)与函数rifft(x,2*n,scale)其实是一个函数，实现同样的功能，使用同cfft（）一样。如果要进行实数fft变换（变换结果实数），则还需调用一个unpacki(x,n)函数。（5）unpacki(x,n)函数功能：对rfft变换后的结果进行变换，为了rifft（）得到原始实数的值。各项参数：x[n]x是一个n（n必须为偶数）项的一维数组，数组中数据定义为短整型（16位有符号整形）。数组x既作为输入数据，又存放变换后的输出数据。n定义为数组中实数的个数，即等于数组大小。函数的使用：可以把这个函数看成unpack（）函数的逆变换，具体原理同上。五、数据测试与实验调试原始数据：为了便于观察，我们对已知信号或已知频谱进行分析。可以利用matlab算出一个63阶低通滤波器单位脉冲响应h[k]的64点值存放到数组x[]。利用N点复序列计算2N点实序列频谱（rfft）之前，2N（即64）点实序列被认为是N点复序列。例如：x[NX]={13,-32,-31,22,52,-16,-84,-9,117,55,-142,-128,147,225,-115,-339,30,455,123,-551,-359,597,691,-549,-1143,341,1775,176,-2825,-1578,5900,13543,13543,5900,-1578,-2825,176,1775,341,-1143,-549,691,597,-359,-551,123,455,30,-339,-115,225,147,-128,-142,55,117,-9,-84,-16,52,22,-31,-32,13}表示一个复数：13+j*（-32）61、实域的波形和理论的频域波形如下：实验的时域波形数据的理论频谱（未经FFT运算）72、倒序：调用函数：cbrev(x,x,NX/2);函数运行的结果如下：在CCS界面下，点击View/WatchWindow，然后在Watch1中的name中键入观察的变量（如变量f）。3、实数fft变换：调用函数rfft(x,64,0);注意：正如五（1）所述，快速傅立叶变换FFT或CFFT运算的结果是频谱的实部和虚部数值，没有作取模运算，因此计算完后，如果需要信号频谱的幅度响应，还要对结果进行取模运算。在CCS界面下，点击View/Graph/Time\Frequency，在图形属性对话框中，如下填写：8图形如下：数据的实际幅度频谱（经过FFT运算）94、实数IFFT变换需要调用三个函数：unpacki(x,NX);解锁cbrev(x,x,NX/2);倒序rifft(x,NX,1);IFFT变换图形如下：频谱经过IFFT后恢复原时域信号的波形10六、数据段的分析1、比较.cmd文件和.map文件，了解各段分配情况。fft_test.cmd文件如下：/*C5402DSKDSPMemoryMap*/MEMORY{PAGE0:VECS:origin=0080h,length=0080h/*InternalProgramRAM*/PRAM:origin=7600h,length=8000h/*InternalProgramRAM*/PAGE1:STACK:origin=1180h,length=0560h/*StackMemorySpace*/EXRAM:origin=1ff1h,length=4000h/*ExternalDataRAM*/}/*DSPMemoryAllocation*/SECTIONS{.cinitPRAMPAGE0.textPRAMPAGE0.vectorsVECSPAGE0.stackSTACKPAGE1.constEXRAMPAGE1.bssEXRAMPAGE1}fft_test.map文件（部分）如下：MEMORYCONFIGURATIONnameor