s实验六_FFT变换

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

2、看出，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）点序列，那么计算N点

3、的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=x2n n=0,1,

4、(N/2)-1x2n=x2n+1 n=0,1,(N/2)-1xn的N点DFT可写成： + (n为奇数)(n为偶数) 因考虑到可写成：=故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

5、AR1,FFT_TWID_P LD #FFT_DP,DPSTM #SYSTEM_STACK,SPrfft_task:STM #FFT_ORIGIN,AR3 ; to dataSTM #data_input,AR4 RPT #K_FFT_SIZE*2-1 ；K_FFT_SIZE=128MVDD *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_prg"bit_rev:SSBX FRCT ;S

9、T #FFT_ORIGIN,d_input_addr; fractional mode is onSTM #FFT_ORIGIN,ORIGINAL_INPUT ; AR3 -> 1 st original inputSTM #fft_data,DATA_PROC_BUF ; AR7 -> data processing bufferMVMM DATA_PROC_BUF,REORDERED_DATA ; AR2 -> 1st bit-reversed dataSTM #K_FFT_SIZE-1,BRCRPTBD bit_rev_end-1STM #K_FFT_SIZE,AR0

10、; AR0 = 1/2 size of circ bufferMVDD *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)0D0

11、3Ha(3)0CFCHR(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 remain

12、ing stages.def fft;.sect "rfft_prg"fft:; Stage 1 -STM #K_ZERO_BK,BK ; BK=0 so that *ARn+0% = *ARn+0LD #-1,ASM ; outputs div by 2 at each stageMVMM DATA_PROC_BUF,PX ; PX -> PRLD *PX,A ; A := PRSTM #fft_data+K_DATA_IDX_1,QX ; QX -> QRSTM #K_FFT_SIZE/2-1,BRCRPTBD stage1end-1STM #K_DATA_

13、IDX_1+1,AR0SUB *QX,16,A,B ; B := PR-QRADD *QX,16,A ; A := PR+QRSTH A,ASM,*PX+ ; PR:= (PR+QR)/2ST B,*QX+ ; QR:= (PR-QR)/2|LD *PX,A ; A := PISUB *QX,16,A,B ; B := PI-QIADD *QX,16,A ; A := PI+QISTH A,ASM,*PX+0 ; PI:= (PI+QI)/2ST B,*QX+0% ; QI:= (PI-QI)/2|LD *PX,A ; A := next PR stage1end:; Stage 2 -MVM

14、M DATA_PROC_BUF,PX ; PX -> PRSTM #fft_data+K_DATA_IDX_2,QX ; QX -> QRSTM #K_FFT_SIZE/4-1,BRCLD *PX,A ; A := PRRPTBD stage2end-1STM #K_DATA_IDX_2+1,AR0; 1st butterflySUB *QX,16,A,B ; B := PR-QRADD *QX,16,A ; A := PR+QRSTH A,ASM,*PX+ ; PR:= (PR+QR)/2ST B,*QX+ ; QR:= (PR-QR)/2|LD *PX,A ; A := PIS

15、UB *QX,16,A,B ; B := PI-QIADD *QX,16,A ; A := PI+QISTH A,ASM,*PX+ ; PI:= (PI+QI)/2STH B,ASM,*QX+ ; QI:= (PI-QI)/2; 2nd butterflyMAR *QX+ADD *PX,*QX,A ; A := PR+QISUB *PX,*QX-,B ; B := PR-QISTH A,ASM,*PX+ ; PR:= (PR+QI)/2SUB *PX,*QX,A ; A := PI-QRST B,*QX ; QR:= (PR-QI)/2|LD *QX+,B ; B := QRST A, *PX

16、 ; PI:= (PI-QR)/2|ADD *PX+0%,A ; A := PI+QRST 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 alwaysST #K_TWID_IDX_3,d_twid_idx ; init index of twiddle tableSTM #K_TWID_IDX_3,AR0 ; AR0 = index of twiddle tableSTM #c

17、osine,WR ; init WR pointerSTM #sine,WI ; init WI pointerSTM #K_LOGN-2-1,STAGE_COUNTER ; init stage counterST #K_FFT_SIZE/8-1,d_grps_cnt ; init group counterSTM #K_FLY_COUNT_3-1,BUTTERFLY_COUNTER ; init butterfly counterST #K_DATA_IDX_3,d_data_idx ; init index for input datastage:STM #fft_data,PX ; P

18、X -> PRLD d_data_idx, AADD *(PX),ASTLM A,QX ; QX -> QRMVDK 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;integer that satisfies

19、 2 N > R. The value R must be loaded into BK.group:MVMD BUTTERFLY_COUNTER,BRC ; # of butterflies in each grpRPTBD butterflyend-1LD *WR,T ; T := WRMPY *QX+,A ; A := QR*WR | QX->QIMACR *WI+0%,*QX-,A ; A := QR*WR+QI*WI | QX->QRADD *PX,16,A,B ; B := (QR*WR+QI*WI)+PRST B,*PX ; PR:=(QR*WR+QI*WI)+

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

21、ST B,*PX+ ; PI:= (PI-(QR*WI-QI*WR)/2 | PX->PR|MPY *QX+,A ; A := QR*WR | QX->QIbutterflyend:; Update pointers for next groupPSHM AR0 ; preserve AR0MVDK d_data_idx,AR0MAR *PX+0 ; increment PX for next groupMAR *QX+0 ; increment QX for next groupBANZD group,*GROUP_COUNTER-POPM AR0 ; restore AR0MA

22、R *QX-; Update counters and indices for next stageLD d_data_idx,ASUB #1,A,B ; B = A-1STLM B,BUTTERFLY_COUNTER ; BUTTERFLY_COUNTER = #flies-1STL A,1,d_data_idx ; double the index of dataLD d_grps_cnt,ASTL A,ASM,d_grps_cnt ; 1/2 the offset to next groupLD d_twid_idx,ASTL A,ASM,d_twid_idx ; 1/2 the ind

23、ex of twiddle tableBANZD 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_prg"unpack:; Compute intermediate values RP, RM,

24、 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-kSTM #fft_data+4*K_FFT_SIZE-1,XM_k ; AR6 -> temp R

25、MkSTM #-2+K_FFT_SIZE/2,BRCRPTBD phase3end-1STM #3,AR0ADD *XP_k,*XP_Nminusk,A ; A := Rk+RN-k =2*RPkSUB *XP_k,*XP_Nminusk,B ; B := Rk-RN-k =2*RMkSTH A,ASM,*XP_k+ ; store RPk at ARkSTH A,ASM,*XP_Nminusk+ ; store RPN-k=RPk at ARN-kSTH B,ASM,*XM_k- ; store RMk at AI2N-kNEG B ; B := RN-k-Rk =2*RMN-kSTH B,

26、ASM,*XM_Nminusk- ; store RMN-k at AIN+kADD *XP_k,*XP_Nminusk,A ; A := Ik+IN-k =2*IPkSUB *XP_k,*XP_Nminusk,B ; B := Ik-IN-k =2*IMkSTH A,ASM,*XP_k+ ; store IPk at AIkSTH A,ASM,*XP_Nminusk-0 ; store IPN-k=IPk at AIN-kSTH B,ASM,*XM_k- ; store IMk at AR2N-kNEG B ; B := IN-k-Ik =2*IMN-kSTH B,ASM,*XM_Nminu

27、sk+0 ; store IMN-k at ARN+kphase3end:ST #0,*XM_k- ; RMN/2=0ST #0,*XM_k ; IMN/2=0; Compute AR0,AI0, ARN, AIN.asg AR2,AX_k.asg AR4,IP_0.asg AR5,AX_NSTM #fft_data,AX_k ; AR2 -> AR0 (tempRP0)STM #fft_data+1,IP_0 ; AR4 -> AI0 (tempIP0)STM #fft_data+2*K_FFT_SIZE+1,AX_N ; AR5 -> AINADD *AX_k,*IP_0

28、,A ; A := RP0+IP0SUB *AX_k,*IP_0,B ; B := RP0-IP0STH A,ASM,*AX_k+ ; AR0 = (RP0+IP0)/2ST #0,*AX_k ; AI0 = 0MVDD *AX_k+,*AX_N- ; AIN = 0STH B,ASM,*AX_N ; ARN = (RP0-IP0)/2; Compute final output values ARk, AIk.asg AR3,AX_2Nminusk.asg AR4,COS.asg AR5,SIN STM #fft_data+4*K_FFT_SIZE-1,AX_2Nminusk ; AR3 -

29、> AI2N-1(temp RM1)STM #cosine+K_TWID_TBL_SIZE/K_FFT_SIZE,COS ; AR4 -> cos(k*pi/N)STM #sine+K_TWID_TBL_SIZE/K_FFT_SIZE,SIN ; AR5 -> sin(k*pi/N)STM #K_FFT_SIZE-2,BRCRPTBD phase4end-1STM #K_TWID_TBL_SIZE/K_FFT_SIZE,AR0 ; index of twiddle tablesLD *AX_k+,16,A ; A := RPk |AR2->IPkMACR *COS,*AX_k,A ; A :=A+cos(k*pi/N)*IPkMASR *SIN,*AX_2Nminusk-,A ; A := A-sin(k*pi/N)*RMk | AR3->IMkLD *AX_2Nminusk+,16,B ; B := IMk |AR3->RMkMASR *SIN+0%,*AX_k-,B ; B :=