2048点的FFT运算

南京航空航天大学研究生DSP综合实验实验报告（2014级）题目：2048点FFT变换的DSP实现学生：蔡静学号：SX14011112015年6月一、实验目的：1、加深对DFT算法原理和基本性质的理解；2、熟悉FFT算法的原理和FFT子程序的算法流程和应用；3、学习用FFT对连续信号和时域信号进行频谱分析的方法；4、学习DSP中FFT的设计和编程思想；5、学习修改FFT程序使之能够进行2048点的FFT运算；二、实验内容：基于时间抽取的8~1024点FFT的程序及输入、输出说明，然后调试出该程序，并对输入方波信号进行FFT变换，检验输出信号是否是一sinc波形。同时更改程序，使其能进行2048点FFT运算。三、实验设计原理：1.离散傅里叶变换DFT:对于长度为N的有限长序列x(n)，它的离散傅里叶变换(DFT)为在x（n）为复数序列的情况下，对于某个K值，直接按（1）式计算x(k)只需N次复数乘法和（N-1）次复数加法。因此对所有N歌K值，共需要N*N次复数乘法和N(N-1)次复数加法。对于一些相当大的点数来说（N=1024），直接计算它的NFT所需的计算量是很大的，因此DFT的运算应用受到很大的限制。2.快速傅里叶变换旋转因子有如下特性，对称性：周期性：FFT算法就是利用旋转因子的周期性和对称性来减少运算量的。FFT算法就是将长序列的DFT分解为短的DFT。由此我们可以看到，一个点DFT已分解成两个N/2点的DFT，这两个N/2点的DFT再按照式（17-2）组合成一个N点DFT。由于X1(k)和X2(k)只有N/2个点，即k=0,1,2,3```,N/2-1,因此利用周期性可以进行进一步减少运算量。假定序列x(n)的点数N是2的幂次，按照FFT算法可将其分为偶序列和奇序列。偶序列：x(2r)=x1(r)奇序列：x(2r+1)=x2(r)其中：r=0,1,2,3,4.......,N/2-1则x(n)的DFT表示为由于对称性。因此N点DFT可分为两部分：前半部分：后半部分：由上式可以看出只要求出0~N/2-1区间x1(k)和x2(k）的值，就可以求出0~N-1区间x(k)的N点值。以同样方式进行抽取，可以求得N/4点的DFT,重复抽取过程，就可以使N点的DFT用上组2点的DFT来计算，这样就可以大大减少运算量。在基数为2的FFT中，设N=2^M,共有M级运算，每级运算有N/2个2点FFT蝶形运算，因此N点FFT总共有MN/2个蝶形运算。图1.蝶形运算原理图例如基数为2的FFT当N=8时，共需要3级，12个基2DFTFFT的蝶形运算，信号流程如图所示：图2.基数为2N=8时的FFT蝶形运算流程图四、FFT算法的实现过程：1.FFT运算序列的存储器分配2.运算实现（1）输入数据的比特反转：将输入数据进行位码倒置以便在整个运算后的输出序列是一个自然序列。（2）实现N点复数FFT：N点FFT算法的实现可分为三个功能模块，即第一级蝶形运算，第二级蝶形运算，第三季蝶形运算至log2N级蝶形运算。（3）功率谱计算，即计算（4）输出FFT结果：3.汇编程序程序主体由bit-rev位码倒置程序，用来实现输入数据的比特反转。fft:FFT算法子程序，用来完成N点FFT运算。在运算过程中为避免结果溢出，对对每个蝶形运算结果右移一位。Fft子程序分唯爱三个功能块：第一季蝶形运算，第二级蝶形运算，第三级至log2N级蝶形运算。4.正弦系数表和余弦系数表：正余弦系数表由数据文件coeff.inc给出，主程序通过.copy汇编命令将正余弦系数表与程序代码汇编在一起。数据文件给出了1024复数的正弦余弦系数各512个，利用此系数可完成8~1024点的FFT运算，通过把正余弦系数各扩充至1024点和修改程序中的K_FFTSIZE和K_LOGN值就可以完成2048点的FFT运算。五、实验结果1.先输入一个方波脉冲信号对其进行64点的FFT变换根据N值将程序中的两个常数修改为：K_FFT_SIZE.set64K_LOGN.set6输入的方波信号图形如下所示：图3.输入的方波脉冲信号波形图图4.输入的方波脉冲信号经FFT变换后输出的sinc波形图5.输入方波脉冲信号波形属性图6.经FFT变换后的信号波形属性由以上输入输出波形可知对于输入的方波信号经该程序的FFT变换后得到了输出的sinc信号波形。2.修改程序使之能够进行2048点的FFT运算1.根据N值将FFT.ASM中的两个常数修改为K_FFT_SIZE.set2048K_LOGN.set11程序的的其他地方修改详见修改后的程序FFT.ASM，相应的注释的数据为以前的8-1024点FFT运算的程序数据。2.准备数据文件d_input,输入数据按实部，虚部，实部，虚部.....顺序排放，本实验输入的是一脉冲方波信号该信号由4096个数据点组成，其中有2048点表示实部，2048个点表示虚部。3.将COEFF.INC中sine1中的数据增至1024点，cosine1中的数据增至1024点，sine1的数据值按如下公式得到：y=SIN((i-1)*3.1415926/1024)*32768i=1,2,3..1024cosine1的数值按如下公式得到：y=CON((i-1)*3.1415926/1024)*32768i=1,2,3..1024在EXCEL中运行以上公式即可得到2048点的COEFF系数表，得到的数据取整后放入COEFF.INC文件中中供主程序调用。4.修改FFT.CMD文件将未定义数据段及程序段的空间增大，具体见修改后的FFT.CMD文件5.输入的数据为一正弦波，该正弦波由4096点组成，实部虚部交替即2048个复数点，每个复数点的虚部为0，输入数据由以下公式在EXCEL中生成。Y=.word&ROUND(SIN((i-1)*2*3.1415926*2/20)*32767,0)*MOD(i,2)i=1,2..4096将生成的数据文件放入sineinputdata.dat文件中，供主程序调用。3.程序运行结果如下图所示：输入的正弦波数据d_input数据显示如下图图7.输入的正弦波形图d_input图8.输出的FFT变换后的波形图fft_out图9.输入数据图像属性图10.输出数据图像属性由输入输出结果可知修改后的fft程序能够进行2048点复数的FFT运算4.实现2048点FFT运算的DSP汇编程序1.FFT运算的汇编主程序FFT.ASM文件:.titleFFT.asm.mmregs.copycoeff.inc;.copyFFTIN64P.dat.copysineinputdata.dat.def_c_int00sine:.usectsine,1024;;;;;;;;;;;;;;;;;;;512cosine:.usectcosine,1024;;;;;;;;;;;;;;;;;;;512d_input:.usectd_input,4096;;;;;;;;;;;;;;;;;;;2048fft_data:.usectfft_data,4096;;;;;;;;;;;;;;;;;;;2048fft_out:.usectfft_out,2048;;;;;;;;;;;;;;;;;;;1024STACK:.usectSTACK,10K_DATA_IDX_1.set2K_DATA_IDX_2.set4K_DATA_IDX_3.set8K_TWID_TBL_SIZE.set1024;;;;;;;;;;;;;;;512K_TWID_IDX_3.set128K_FLY_COUNT_3.set4K_FFT_SIZE.set2048;N=8K_LOGN.set11;LOG(N)=LOG(32)=5PA0.set0PA1.set1.bssd_twid_idx,1.bssd_data_idx,1.bssd_grps_cnt,1.sectfft_prg;****位码倒序**************************************.asgAR2,REORDERED.asgAR3,ORIGINAL_INPUT.asgAR7,DATA_PROC_BUF.text_c_int00:SSBXFRCTSTM#STACK+10,SPSTM#d_input,AR1;从PA1口输入2N个数据到d_input中;STM#fft_data,AR1;inputthe2NdataRPT#2*K_FFT_SIZE-1;fromPA1MVPDd_input1,*AR1+;PORTRPA1,*AR1+STM#sine,AR1;movecoeftofsinRPT#1023;;;;;511MVPDsine1,*AR1+;fromprogramtodataSTM#cosine,AR1;movecoeffofcosRPT#1023;;;511MVPDcosine1,*AR1+;fromprogramtodataSTM#d_input,ORIGINAL_INPUTSTM#fft_data,DATA_PROC_BUFMVMMDATA_PROC_BUF,REORDEREDSTM#K_FFT_SIZE-1,BRCRPTBDbit_rev_end-1STM#K_FFT_SIZE,AR0MVDD*ORIGINAL_INPUT+,*REORDERED+MVDD*ORIGINAL_INPUT-,*REORDERED+MAR*ORIGINAL_INPUT+0Bbit_rev_end:;****FFT代码********.asgAR1,GROUP_COUNTER.asgAR2,PX.asgAR3,QX.asgAR4,WR.asgAR5,WI.asgAR6,BUTTERFLY_COUNTER.asgAR7,STAGE_COUNTER;***stage1*********STM#0,BK;LD#0,ASMLD#-1,ASM;protectflowover,alloutput/2STM#fft_data,PXLD*PX,16,A;AH:=Re[x(0)]STM#fft_data+K_DATA_IDX_1,QX;QXpointRe[x(4)]STM#K_FFT_SIZE/2-1,BRC;stage1isN/2-1RPTBDstage1end-1STM#K_DATA_IDX_1+1,AR0SUB*QX,16,A,B;BH=Re[x(0)]-Re[x(4)]ADD*QX,16,A;AH=Re[x(0)]+Re[x(4)]STHA,ASM,*PX+STB,*QX+||LD*PX,ASUB*QX,16,A,B;BH=Im[x(0)]-Im[x(4)]ADD*QX,16,A;AH=Im[x(0)]+Im[x(4)]STHA,ASM,*PX+0STB,*QX+0%;BK=0whyherecircleaccess||LD*PX,Astage1end:;***Stage2***STM#fft_data,PXSTM#fft_data+K_DATA_IDX_2,QXSTM#K_FFT_SIZE/4-1,BRCLD*PX,16,A;AH=Re[x(0)]+Re[x(4)],AH=1RPTBDstage2end-1STM#K_DATA_IDX_2+1,AR0;1stbutterflySUB*QX,16,A,BADD*QX,16,ASTHA,ASM,*PX+STB,*QX+||LD*PX,ASUB*QX,16,A,B;BH={Im[x(0)]+Im[x(4)]}-{Im[x(2)]+Im[x(6)]},BH=0ADD*QX,16,A;AH={Im[x(0)]+Im[x(4)]}+{Im[x(2)]+Im[x(6)]},BH=0STHA,ASM,*PX+STHB,ASM,*QX+;2ndbutterflyMAR*QX+ADD*PX,*QX,A;AH=+,=1SUB*PX,*QX-,B;BH=-,=1STHA,ASM,*PX+SUB*PX,*QX,A;AH=-,=-1,note:subSTB,*QX||