s实验六FFT变换.docx

1、s实验六FFT变换实验四 FFT实验一 实验目的（1） 了解FFT的原理；（2） 了解在DSP中FFT的设计及编程方法；（3） 熟悉对FFT的调试方法二 实验内容本实验要求使用FFT变换求一个时域信号的频域特性，并从这个频域特性求出该信号的频域值。使用DSP汇编语言实现对FFT的DSP编程。三 实验原理对于有限长离散数字信号xn,0，它的频谱离散数学值X(K)可由离散傅立叶变换（DFT）求得。DFT定义为： 也可以方便的把它改写成如下形式：式中（有时简写为）代表。不难看出，是周期性的，且周期为N，即 m,l=0,的周期性是DFT的关键之一。为了强调起见，常用表达式取代以便明确地给出地周期为N。

2、由DFT地定义可以看出，xn为复数序列地情况下，完全可以直接运算N点DFT需要次复数乘法和N(N-1)次复数加法。因此，对一些相当大的N值（如1024点）来说，直接计算它的DFT所需要的计算量是很大地。一个优化的实数FFT算法是一个组合后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，然后对复序列进行N点的FFT运算，最后再由N点复数输出拆散成2N点的复数序列，这2N点的复数序列与原始的2N点的实数输入序列的DFT一致。FFT的基本思想在于：将原来的N点序列分成两个较短的序列，这些序列的DFT可很简单地组合起来得到原序列的DFT。例如，若N为偶数，将原有的N点系列分成两个（N/2）点

3、序列，那么计算N点的DFT将只需要（N/2）2*2=N2/2次复数乘法，即比直接计算少作一半乘法。因子(N/2)2表示直接计算（N/2）点DFT所需要地乘法次数，而乘数2代表必须完成两个DFT。上述处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定 N/2为偶数），从而又少作一半乘法，使用这种方法，在组合输入和拆散输出地操作中，FFT地运算量减半。这样，利用实数FFT算法来计算实输入序列地DFT的速度几乎是一般复FFT算法的两倍。假定序列xn的点数N是2的幂，按上述处理方法，定义两个分别为xn的偶数项和奇数项的(N/2)点序列x1n和x2n,即：x1n=x

4、2n n=0,1,(N/2)-1x2n=x2n+1 n=0,1,(N/2)-1xn的N点DFT可写成： + 因考虑到可写成：=故X(k)可写为： + =X1(k)+ X2(k)式中X1(k)和X2(k)是x1(n)和x2(n)的（N/2）点DFT。上式表明，N点DFT 可分解为按上式地规则加以组合地两个（N/2）点DFT。四FFT的DSP编程FFT地实现步骤主要有以下四步：（1） 将输入序列压缩和位倒序。（2） N点的复数FFT。（3） 奇数号部分和偶数号部分分离。（4） 产生最后的输出数据。最初，将原始输入的2N点的实序列a(n),存储在4N大小的数据处理缓冲区的下半部分。.asg AR1,

5、FFT_TWID_P LD #FFT_DP,DP STM #SYSTEM_STACK,SPrfft_task: STM #FFT_ORIGIN,AR3 ; to data STM #data_input,AR4 RPT #K_FFT_SIZE*2-1 ；K_FFT_SIZE=128 MVDD *AR4+,*AR3+0C00h0C01h0CFEh0CFFh0D00ha(0)0D01ha(1) 0DFFEh a(254)0DFFFh a(255)然后通过子程序调用，有以下几步计算驱除原始数据的直流分量。由于采样的需要，信号发生器输出的为叠加了直流分量的正弦信号，所以有此一步。 .asg AR3,F

6、FT_INPUT.def data_avedata_ave:STM #FFT_ORIGIN,FFT_INPUT ;AR3 1st originalinputLD #0,ARPT #K_FFT_SIZE*2-1,AADD *FFT_INPUT+,A ;累加原始数据到ANOPLD A,#-1-k_LOGN ;A右移8位，即A/256，输入数据取平均STM #K_FFT_SIZE*2-1,BRCSTM #FFT_ORIGIN,FFT_INPUT ;AR3-1st original inputRPTB ave_end-1SUB FFT_INPUT,A,B ;B=A-(FFT_INPUT)NEG B ；

7、B=(*FFT_INPUT)-A,即将原始数据去除直流分量STL B,-1,*FFT_INPUT+ ;除以2以防止溢出后存入原址ave_endNOPRET将输入序列位倒序存储，以便算法结束后的输出序列为自然顺序。首先，将2N点的实输入序列拷贝到标记为FFT_ORIGIN的连续内存空间中，理解为一个N点的复数序列d(n)。时序列中索引为偶数的作为d(n)的实部，索引为奇数的作为d(n)的虚部。这一部分就是压缩，然后，将得到的复数序列位倒序存储到数据处理缓冲区fft_data.0C00HR(0)=a(0)0D00Ha(0)0C01HI(0)=a(1)0D01Ha(1)0C02HI(64)=a(12

8、8)0D02Ha(2)0C03HI(64)=a(129)0D03Ha(3)0CFCHR(63)=a(126)0DFCHa(252)0CFDHR(63)=a(127)0DFDHa(253)0CFEHR(127)=a(254)0DFEHa(254)0CFFHI(127)=a(255)0DFFHa(255);Bit Reversal Routine .asg AR2,REORDERED_DATA .asg AR3,ORIGINAL_INPUT .asg AR7,DATA_PROC_BUF .def bit_rev; .sect rfft_prgbit_rev: SSBX FRCT ; ST #FFT

9、_ORIGIN,d_input_addr ; fractional mode is on STM #FFT_ORIGIN,ORIGINAL_INPUT ; AR3 - 1 st original input STM #fft_data,DATA_PROC_BUF ; AR7 - data processing buffer MVMM DATA_PROC_BUF,REORDERED_DATA ; AR2 - 1st bit-reversed data STM #K_FFT_SIZE-1,BRC RPTBD bit_rev_end-1 STM #K_FFT_SIZE,AR0 ; AR0 = 1/2

10、 size of circ buffer MVDD *ORIGINAL_INPUT+,*REORDERED_DATA+ MVDD *ORIGINAL_INPUT-,*REORDERED_DATA+ MAR *ORIGINAL_INPUT+0Bbit_rev_end: RET ; return to Real FFT main module一个N点的复数序列存储在数据处理缓冲区对d(n)作fft变换，结果得到D(k)=fd(n)=R(k)+jI(k)。0C00HR(0)0D00Ha(0)0C01HI(0)0D01Ha(1)0C02HI(1)0D02Ha(2)0C03HI(1)0D03Ha(3)0

11、CFCHR(126)0DFCHa(252)0CFDHR(126)0DFDHa(253)0CFEHR(127)0DFEHa(254)0CFFHI(127)0DFFHa(255)；= ;256-Point Real FFT Routine .asg AR1,GROUP_COUNTER .asg AR2,PX .asg AR3,QX .asg AR4,WR .asg AR5,WI .asg AR6,BUTTERFLY_COUNTER .asg AR7,DATA_PROC_BUF ; for Stages 1 & 2 .asg AR7,STAGE_COUNTER ; for the remaining

12、 stages .def fft; .sect rfft_prgfft:; Stage 1 - STM #K_ZERO_BK,BK ; BK=0 so that *ARn+0% = *ARn+0 LD #-1,ASM ; outputs div by 2 at each stage MVMM DATA_PROC_BUF,PX ; PX - PR LD *PX,A ; A := PR STM #fft_data+K_DATA_IDX_1,QX ; QX - QR STM #K_FFT_SIZE/2-1,BRC RPTBD stage1end-1 STM #K_DATA_IDX_1+1,AR0 S

13、UB *QX,16,A,B ; B := PR-QR ADD *QX,16,A ; A := PR+QR STH A,ASM,*PX+ ; PR:= (PR+QR)/2 ST B,*QX+ ; QR:= (PR-QR)/2 |LD *PX,A ; A := PI SUB *QX,16,A,B ; B := PI-QI ADD *QX,16,A ; A := PI+QI STH A,ASM,*PX+0 ; PI:= (PI+QI)/2 ST B,*QX+0% ; QI:= (PI-QI)/2 |LD *PX,A ; A := next PR stage1end:; Stage 2 - MVMM

14、DATA_PROC_BUF,PX ; PX - PR STM #fft_data+K_DATA_IDX_2,QX ; QX - QR STM #K_FFT_SIZE/4-1,BRC LD *PX,A ; A := PR RPTBD stage2end-1 STM #K_DATA_IDX_2+1,AR0; 1st butterfly SUB *QX,16,A,B ; B := PR-QR ADD *QX,16,A ; A := PR+QR STH A,ASM,*PX+ ; PR:= (PR+QR)/2 ST B,*QX+ ; QR:= (PR-QR)/2 |LD *PX,A ; A := PI

15、SUB *QX,16,A,B ; B := PI-QI ADD *QX,16,A ; A := PI+QI STH A,ASM,*PX+ ; PI:= (PI+QI)/2 STH B,ASM,*QX+ ; QI:= (PI-QI)/2; 2nd butterfly MAR *QX+ ADD *PX,*QX,A ; A := PR+QI SUB *PX,*QX-,B ; B := PR-QI STH A,ASM,*PX+ ; PR:= (PR+QI)/2 SUB *PX,*QX,A ; A := PI-QR ST B,*QX ; QR:= (PR-QI)/2 |LD *QX+,B ; B :=

16、QR ST A, *PX ; PI:= (PI-QR)/2 |ADD *PX+0%,A ; A := PI+QR ST A,*QX+0% ; QI:= (PI+QR)/2 |LD *PX,A ; A := PRstage2end:; Stage 3 thru Stage logN-1 - STM #K_TWID_TBL_SIZE,BK ; BK = twiddle table size always ST #K_TWID_IDX_3,d_twid_idx ; init index of twiddle table STM #K_TWID_IDX_3,AR0 ; AR0 = index of t

17、widdle table STM #cosine,WR ; init WR pointer STM #sine,WI ; init WI pointer STM #K_LOGN-2-1,STAGE_COUNTER ; init stage counter ST #K_FFT_SIZE/8-1,d_grps_cnt ; init group counter STM #K_FLY_COUNT_3-1,BUTTERFLY_COUNTER ; init butterfly counter ST #K_DATA_IDX_3,d_data_idx ; init index for input datast

18、age: STM #fft_data,PX ; PX - PR LD d_data_idx, A ADD *(PX),A STLM A,QX ; QX - QR MVDK d_grps_cnt,GROUP_COUNTER ; AR1 contains group counter;A circular buffer of size R must start on a N-bit boundary (that is, the N LSBs;of the base address of the circular buffer must be 0), where N is the smallest;i

19、nteger that satisfies 2 N R. The value R must be loaded into BK.group: MVMD BUTTERFLY_COUNTER,BRC ; # of butterflies in each grp RPTBD butterflyend-1 LD *WR,T ; T := WR MPY *QX+,A ; A := QR*WR | QX-QI MACR *WI+0%,*QX-,A ; A := QR*WR+QI*WI | QX-QR ADD *PX,16,A,B ; B := (QR*WR+QI*WI)+PR ST B,*PX ; PR:

20、=(QR*WR+QI*WI)+PR)/2 |SUB *PX+,B ; B := PR-(QR*WR+QI*WI) | PX-PI ST B,*QX ; QR:= (PR-(QR*WR+QI*WI)/2 |MPY *QX+,A ; A := QR*WI T=WI | QX-QI MASR *QX,*WR+0%,A ; A := QR*WI-QI*WR ADD *PX,16,A,B ; B := (QR*WI-QI*WR)+PI ST B,*QX+ ; QI:=(QR*WI-QI*WR)+PI)/2 | QX-QR |SUB *PX,B ; B := PI-(QR*WI-QI*WR) LD *WR

21、,T ; T := WR ST B,*PX+ ; PI:= (PI-(QR*WI-QI*WR)/2 | PX-PR |MPY *QX+,A ; A := QR*WR | QX-QIbutterflyend:; Update pointers for next group PSHM AR0 ; preserve AR0 MVDK d_data_idx,AR0 MAR *PX+0 ; increment PX for next group MAR *QX+0 ; increment QX for next group BANZD group,*GROUP_COUNTER- POPM AR0 ; r

22、estore AR0 MAR *QX-; Update counters and indices for next stage LD d_data_idx,A SUB #1,A,B ; B = A-1 STLM B,BUTTERFLY_COUNTER ; BUTTERFLY_COUNTER = #flies-1 STL A,1,d_data_idx ; double the index of data LD d_grps_cnt,A STL A,ASM,d_grps_cnt ; 1/2 the offset to next group LD d_twid_idx,A STL A,ASM,d_t

23、wid_idx ; 1/2 the index of twiddle table BANZD stage,*STAGE_COUNTER- MVDK d_twid_idx,AR0 ; AR0 = index of twiddle tablefft_end: RET ; return to Real FFT main module将fft的计算结果分离为RP（偶实部），RM（奇实部），IP（偶实部），IM（奇实部）;=;Unpack 256-Point Real FFT Output .def unpack; .sect rfft_prgunpack:; Compute intermediate

24、values RP, RM, IP, IM .asg AR2,XP_k .asg AR3,XP_Nminusk .asg AR6,XM_k .asg AR7,XM_Nminusk STM #fft_data+2,XP_k ; AR2 - Rk (temp RPk) STM #fft_data+2*K_FFT_SIZE-2,XP_Nminusk ; AR3 - RN-k (tempRPN-k) STM #fft_data+2*K_FFT_SIZE+3,XM_Nminusk ; AR7 - temp RMN-k STM #fft_data+4*K_FFT_SIZE-1,XM_k ; AR6 - t

25、emp RMk STM #-2+K_FFT_SIZE/2,BRC RPTBD phase3end-1 STM #3,AR0 ADD *XP_k,*XP_Nminusk,A ; A := Rk+RN-k =2*RPk SUB *XP_k,*XP_Nminusk,B ; B := Rk-RN-k =2*RMk STH A,ASM,*XP_k+ ; store RPk at ARk STH A,ASM,*XP_Nminusk+ ; store RPN-k=RPk at ARN-k STH B,ASM,*XM_k- ; store RMk at AI2N-k NEG B ; B := RN-k-Rk

26、=2*RMN-k STH B,ASM,*XM_Nminusk- ; store RMN-k at AIN+k ADD *XP_k,*XP_Nminusk,A ; A := Ik+IN-k =2*IPk SUB *XP_k,*XP_Nminusk,B ; B := Ik-IN-k =2*IMk STH A,ASM,*XP_k+ ; store IPk at AIk STH A,ASM,*XP_Nminusk-0 ; store IPN-k=IPk at AIN-k STH B,ASM,*XM_k- ; store IMk at AR2N-k NEG B ; B := IN-k-Ik =2*IMN-k STH B,ASM,*XM_Nminusk+0 ; store IMN-k at ARN+kphase3end: ST #0,*XM_k- ; RMN/2=0 ST #0,*XM_k ; IMN/2=0