s实验六FFT变换.docx
《s实验六FFT变换.docx》由会员分享，可在线阅读，更多相关《s实验六FFT变换.docx（16页珍藏版）》请在冰豆网上搜索。
s实验六FFT变换
实验四FFT实验
一实验目的
（1）了解FFT的原理；
（2）了解在DSP中FFT的设计及编程方法；
（3）熟悉对FFT的调试方法
二实验内容
本实验要求使用FFT变换求一个时域信号的频域特性，并从这个频域特性求出该信号的频域值。
使用DSP汇编语言实现对FFT的DSP编程。
三实验原理
对于有限长离散数字信号{x[n]},0
，它的频谱离散数学值{X（K）}可由离散傅立叶变换（DFT）求得。
DFT定义为：
也可以方便的把它改写成如下形式：
式中
（有时简写为
）代表
。
不难看出，
是周期性的，且周期为N，即
m,l=0,
的周期性是DFT的关键之一。
为了强调起见，常用表达式
取代
以便明确地给出
地周期为N。
由DFT地定义可以看出，x[n]为复数序列地情况下，完全可以直接运算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点的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算法的两倍。
假定序列{x[n]}的点数N是2的幂，按上述处理方法，定义两个分别为x[n]的偶数项和奇数项的（N/2）点序列x1[n]和x2[n],即：
x1[n]=x[2n]n=0,1,…,（N/2）-1
x2[n]=x[2n+1]n=0,1,…,（N/2）-1
{x[n]}的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大小的数据处理缓冲区的下半部分。
.asgAR1,FFT_TWID_P
LD#FFT_DP,DP
STM#SYSTEM_STACK,SP
rfft_task:
STM#FFT_ORIGIN,AR3;todata
STM#data_input,AR4
RPT#K_FFT_SIZE*2-1
；K_FFT_SIZE=128
MVDD*AR4+,*AR3+
0C00h
0C01h
…
0CFEh
0CFFh
0D00h
a（0）
0D01h
a
（1）
…
…
0DFFEh
a（254）
0DFFFh
a（255）
然后通过子程序调用，有以下几步计算
驱除原始数据的直流分量。
由于采样的需要，信号发生器输出的为叠加了直流分量的正弦信号，所以有此一步。
.asgAR3,FFT_INPUT
.defdata_ave
data_ave:
STM#FFT_ORIGIN,FFT_INPUT;AR31storiginalinput
LD#0,A
RPT#K_FFT_SIZE*2-1,A
ADD*FFT_INPUT+,A;累加原始数据到A
NOP
LDA,#-1-k_LOGN;A右移8位，即A/256，输入数据取平均
STM#K_FFT_SIZE*2-1,BRC
STM#FFT_ORIGIN,FFT_INPUT;AR3---1storiginalinput
RPTBave_end-1
SUBFFT_INPUT,A,B;B=A-（FFT_INPUT）
NEGB；B=（*FFT_INPUT）-A,即将原始数据去除直流分量
STLB,-1,*FFT_INPUT+;除以2以防止溢出后存入原址
ave_end
NOP
RET
将输入序列位倒序存储，以便算法结束后的输出序列为自然顺序。
首先，将2N点的实输入序列拷贝到标记为FFT_ORIGIN的连续内存空间中，理解为一个N点的复数序列d（n）。
时序列中索引为偶数的作为d（n）的实部，，索引为奇数的作为d（n）的虚部。
这一部分就是压缩，然后，将得到的复数序列位倒序存储到数据处理缓冲区fft_data.
0C00H
R（0）=a（0）
0D00H
a（0）
0C01H
I（0）=a
（1）
0D01H
a
（1）
0C02H
I（64）=a（128）
0D02H
a
（2）
0C03H
I（64）=a（129）
0D03H
a（3）
……
……
……
……
0CFCH
R（63）=a（126）
0DFCH
a（252）
0CFDH
R（63）=a（127）
0DFDH
a（253）
0CFEH
R（127）=a（254）
0DFEH
a（254）
0CFFH
I（127）=a（255）
0DFFH
a（255）
;BitReversalRoutine
.asgAR2,REORDERED_DATA
.asgAR3,ORIGINAL_INPUT
.asgAR7,DATA_PROC_BUF
.defbit_rev
;.sect"rfft_prg"
bit_rev:
SSBXFRCT
;ST#FFT_ORIGIN,d_input_addr;fractionalmodeison
STM#FFT_ORIGIN,ORIGINAL_INPUT;AR3->1storiginalinput
STM#fft_data,DATA_PROC_BUF;AR7->dataprocessingbuffer
MVMMDATA_PROC_BUF,REORDERED_DATA;AR2->1stbit-reverseddata
STM#K_FFT_SIZE-1,BRC
RPTBDbit_rev_end-1
STM#K_FFT_SIZE,AR0;AR0=1/2sizeofcircbuffer
MVDD*ORIGINAL_INPUT+,*REORDERED_DATA+
MVDD*ORIGINAL_INPUT-,*REORDERED_DATA+
MAR*ORIGINAL_INPUT+0B
bit_rev_end:
RET;returntoRealFFTmainmodule
一个N点的复数序列存储在数据处理缓冲区
对d（n）作fft变换，结果得到D（k）=f{d（n）}=R（k）+jI（k）。
0C00H
R（0）
0D00H
a（0）
0C01H
I（0）
0D01H
a
（1）
0C02H
I
（1）
0D02H
a
（2）
0C03H
I
（1）
0D03H
a（3）
……
……
……
……
0CFCH
R（126）
0DFCH
a（252）
0CFDH
R（126）
0DFDH
a（253）
0CFEH
R（127）
0DFEH
a（254）
0CFFH
I（127）
0DFFH
a（255）
；=====================================================
;256-PointRealFFTRoutine
.asgAR1,GROUP_COUNTER
.asgAR2,PX
.asgAR3,QX
.asgAR4,WR
.asgAR5,WI
.asgAR6,BUTTERFLY_COUNTER
.asgAR7,DATA_PROC_BUF;forStages1&2
.asgAR7,STAGE_COUNTER;fortheremainingstages
.deffft
;.sect"rfft_prg"
fft:
;Stage1--------------------------------------------
STM#K_ZERO_BK,BK;BK=0sothat*ARn+0%==*ARn+0
LD#-1,ASM;outputsdivby2ateachstage
MVMMDATA_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
RPTBDstage1end-1
STM#K_DATA_IDX_1+1,AR0
SUB*QX,16,A,B;B:
=PR-QR
ADD*QX,16,A;A:
=PR+QR
STHA,ASM,*PX+;PR’:
=（PR+QR）/2
STB,*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
STHA,ASM,*PX+0;PI’:
=（PI+QI）/2
STB,*QX+0%;QI’:
=（PI-QI）/2
||LD*PX,A;A:
=nextPR
stage1end:
;Stage2----------------------------------------------------------------------
MVMMDATA_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
RPTBDstage2end-1
STM#K_DATA_IDX_2+1,AR0
;1stbutterfly
SUB*QX,16,A,B;B:
=PR-QR
ADD*QX,16,A;A:
=PR+QR
STHA,ASM,*PX+;PR’:
=（PR+QR）/2
STB,*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
STHA,ASM,*PX+;PI’:
=（PI+QI）/2
STHB,ASM,*QX+;QI’:
=（PI-QI）/2
;2ndbutterfly
MAR*QX+
ADD*PX,*QX,A;A:
=PR+QI
SUB*PX,*QX-,B;B:
=PR-QI
STHA,ASM,*PX+;PR’:
=（PR+QI）/2
SUB*PX,*QX,A;A:
=PI-QR
STB,*QX;QR’:
=（PR-QI）/2
||LD*QX+,B;B:
=QR
STA,*PX;PI’:
=（PI-QR）/2
||ADD*PX+0%,A;A:
=PI+QR
STA,*QX+0%;QI’:
=（PI+QR）/2
||LD*PX,A;A:
=PR
stage2end:
;Stage3thruStagelogN-1----------------------------------------------------
STM#K_TWID_TBL_SIZE,BK;BK=twiddletablesizealways
ST#K_TWID_IDX_3,d_twid_idx;initindexoftwiddletable
STM#K_TWID_IDX_3,AR0;AR0=indexoftwiddletable
STM#cosine,WR;initWRpointer
STM#sine,WI;initWIpointer
STM#K_LOGN-2-1,STAGE_COUNTER;initstagecounter
ST#K_FFT_SIZE/8-1,d_grps_cnt;initgroupcounter
STM#K_FLY_COUNT_3-1,BUTTERFLY_COUNTER;initbutterflycounter
ST#K_DATA_IDX_3,d_data_idx;initindexforinputdata
stage:
STM#fft_data,PX;PX->PR
LDd_data_idx,A
ADD*（PX）,A
STLMA,QX;QX->QR
MVDKd_grps_cnt,GROUP_COUNTER;AR1containsgroupcounter
;AcircularbufferofsizeRmuststartonaN-bitboundary（thatis,theNLSBs
;ofthebaseaddressofthecircularbuffermustbe0）,whereNisthesmallest
;integerthatsatisfies2N>R.ThevalueRmustbeloadedintoBK.
group:
MVMDBUTTERFLY_COUNTER,BRC;#ofbutterfliesineachgrp
RPTBDbutterflyend-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
STB,*PX;PR’:
=（（QR*WR+QI*WI）+PR）/2
||SUB*PX+,B;B:
=PR-（QR*WR+QI*WI）||PX->PI
STB,*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
STB,*QX+;QI’:
=（（QR*WI-QI*WR）+PI）/2||QX->QR
||SUB*PX,B;B:
=PI-（QR*WI-QI*WR）
LD*WR,T;T:
=WR
STB,*PX+;PI’:
=（PI-（QR*WI-QI*WR））/2||PX->PR
||MPY*QX+,A;A:
=QR*WR||QX->QI
butterflyend:
;Updatepointersfornextgroup
PSHMAR0;preserveAR0
MVDKd_data_idx,AR0
MAR*PX+0;incrementPXfornextgroup
MAR*QX+0;incrementQXfornextgroup
BANZDgroup,*GROUP_COUNTER-
POPMAR0;restoreAR0
MAR*QX-
;Updatecountersandindicesfornextstage
LDd_data_idx,A
SUB#1,A,B;B=A-1
STLMB,BUTTERFLY_COUNTER;BUTTERFLY_COUNTER=#flies-1
STLA,1,d_data_idx;doubletheindexofdata
LDd_grps_cnt,A
STLA,ASM,d_grps_cnt;1/2theoffsettonextgroup
LDd_twid_idx,A
STLA,ASM,d_twid_idx;1/2theindexoftwiddletable
BANZDstage,*STAGE_COUNTER-
MVDKd_twid_idx,AR0;AR0=indexoftwiddletable
fft_end:
RET;returntoRealFFTmainmodule
将fft的计算结果分离为RP（偶实部），RM（奇实部），IP（偶实部），IM（奇实部）
;=================================================
;Unpack256-PointRealFFTOutput
.defunpack
;.sect"rfft_prg"
unpack:
;ComputeintermediatevaluesRP,RM,IP,IM
.asgAR2,XP_k
.asgAR3,XP_Nminusk
.asgAR6,XM_k
.asgAR7,XM_Nminusk
STM#fft_data+2,XP_k;AR2->R[k]（tempRP[k]）
STM#fft_data+2*K_FFT_SIZE-2,XP_Nminusk;AR3->R[N-k]（tempRP[N-k]）
STM#fft_data+2*K_FFT_SIZE+3,XM_Nminusk;AR7->tempRM[N-k]
STM#fft_data+4*K_FFT_SIZE-1,XM_k;AR6->tempRM[k]
STM#-2+K_FFT_SIZE/2,BRC
RPTBDphase3end-1
STM#3,AR0
ADD*XP_k,*XP_Nminusk,A;A:
=R[k]+R[N-k]=2*RP[k]
SUB*XP_k,*XP_Nminusk,B;B:
=R[k]-R[N-k]=2*RM[k]
STHA,ASM,*XP_k+;storeRP[k]atAR[k]
STHA,ASM,*XP_Nminusk+;storeRP[N-k]=RP[k]atAR[N-k]
STHB,ASM,*XM_k-;storeRM[k]atAI[2N-k]
NEGB;B:
=R[N-k]-R[k]=2*RM[N-k]
STHB,ASM,*XM_Nminusk-;storeRM[N-k]atAI[N+k]
ADD*XP_k,*XP_Nminusk,A;A:
=I[k]+I[N-k]=2*IP[k]
SUB*XP_k,*XP_Nminusk,B;B:
=I[k]-I[N-k]=2*IM[k]
STHA,ASM,*XP_k+;storeIP[k]atAI[k]
STHA,ASM,*XP_Nminusk-0;storeIP[N-k]=IP[k]atAI[N-k]
STHB,ASM,*XM_k-;storeIM[k]atAR[2N-k]
NEGB;B:
=I[N-k]-I[k]=2*IM[N-k]
STHB,ASM,*XM_Nminusk+0;storeIM[N-k]atAR[N+k]
phase3end:
ST#0,*XM_k-;RM[N/2]=0
ST#0,*XM_k;IM[N/2]=0
升级会员 