1、快速傅立叶变换（FFT）程序设计DFT(Discrete Fourier Transformation)是数字信号分析与处理，如图形、语音及图像等领域的重要变换工具，直接计算DFT的计算量与变换区间长度N的平方成正比。当N较大时，因计算量太大，直接用DFT算法进行谱分析和信号的实时处理是不切实际的。快速傅立叶变换(Fast Fourier Transformation，简称FFT)使DFT运算效率提高12个数量级，其原因是当N较大时，对DFT进行了基4和基2分解运算。FFT 算法利用DFT系数的对称性、周期性和可约性等性质将长序列的DFT分解为若干个短序列的DFT计算，然后再按一定规则将其合并
2、，从而得到整个的DFT。按时间抽取的基2 FFT算法是FFT最常用的一种表现形式。设 （ =2）点序列 的离散傅立叶变换为 ，则有：()xn()Xk(1)210NjnkNnXke(2)210()()jnknx式中： =0,1，, ； =0,1, 。nk1将 ( =0,1, )按 的奇偶分成以下两组：)N12(xn2)将上式带入(1)式并化简，就得到(3)和(4) 式。这样N 点序列 就变为两个 点序()Xk2N列 和 离散傅立叶变换的组合：1()xn2(3)12()()KNXkWXk(4)2式中： =0,1, ，k其中 12120()()NnknXkx1220()()nkNnWjKNe(3)式
3、和(4)式可以用下面的碟形节信号流图图1来表示。1()Xk2KNW12()()KNXkWk图1 2点碟形节对于 和 可以继续利用(3)和(4) 式进一步分解，直到分解为对应图1的基本1()Xk2()的两点变换为止。现以8点序列的FFT为例，通过其碟形流图来分析它的数字特点（图 2）。图2 8点按时间抽选法FFT 运算流图通过分析，不失一般性，可以得到 点序列FFT碟形图的特点如下。这里需要说明的N是下面所提到的“ 同种碟形” 即具有相同碟形节系数( )的基本两点碟形。KNW(1) 点碟形图共有 级( )NM2(2)对于第 ( )级m0不同种基本碟形的种数 ；1m相邻不同种基本碟形的碟形 节系数
4、的增量；MmN同种基本碟形的个数 ；2相邻同种基本碟形的间距 ；每个基本碟形中两接点的间距 ；1m(3)输入序列为倒序排列，输出序列为自然序排列。为了说明这个特点，将输入输出序列的序号变为2进制，列表1比较。表1 输入倒序和输出正序对照表x(0) x(4) x(2) x(6) x(1) x(5) x(3) x(7)x(n)000 100 010 110 001 101 011 111X(k) X(0) X(1) X(2) X(3) X(4) X(5) X(6) X(7)000 001 010 011 100 101 110 111从表1中不难看出所谓倒序排列就是指序列的2进制表示的序号采用高位
5、加1并且反向进位方式递增的排列。所以在进行FFT之前要对输入序列进行处理，是输出序列按正序排列。简易程序流程图如图3所示：D S P 初始化调入系数表调入输入数据F F T 蝶形运算归一化N O M = ？N O M = 0N O M = 1图3 FFT程序流程图FFT程序及说明N .set 512 .global _resave .global _input ;存放 FFT 运算中用到的数据,按二进制反序排列.global _indatr ;512 点 FFT 所需数据;-; 反序排列子程序部分;-_resave：MAR * ，AR3LAR AR2，#_inputLAR AR3，#_inda
6、trLAR AR0，#NLAR AR4，#(N-1) RESAV1LACC *+，0 ，AR2SACL *BR0+，AR4BANZ RESAV1，*-，AR3MAR *，AR1SBRK #2LAR AR0，*-PSHD *RET;-; FFT 应用子程序;-N .set 512 ;点数M .set 9 ;N=2M_input .usect “.data0“，2*N ;输入数据 .bss _sintab，N ;SIN 和 COSIN 函数的存储表.bss _nom，1 ;当 _nom=1 时， FFT 需要归一化处理，为 0 时则不需要.global _fft.global _sintab.gl
7、obal _input .global _nom.text_fft： ;-; 与 C 语言兼容的代码部分 ;-POPD *+ ;存储返回地址 ADDRESSSAR AR6，*+ ;存储 AR6SAR AR7，*+ ;存储 AR7SAR AR0，*+ ;存储 AR0SAR AR1，* ;堆栈分布情况：ADDRESS/AR6/AR7/AR0/AR1;ARP=AR1， AR1：AR1LAR AR0，#08h LAR AR3，*0+，AR3 ;AR3：FP， SP=SP+size(size=frame 的长度);ARP=AR3LAR AR2，* ;AR2：AR1LAR AR7，#_nom ;AR7 指
8、向 _nom;-; 初始化一些寄存器;- SPLK #(N-1)，*+ SPLK #(M-1)，*+ ;堆栈分布情况：ADDRESS/AR6/AR7;/AR0/N/M/Y(其中 Y 为没有确定的量 );ARP=AR3， AR2：N， AR3：YSPLK #1，*+ ，AR2 ;ID=1，ARP=AR2;堆栈分布情况：ADDRESS/AR6/AR7/AR0/N/M/ID/Y;ARP=AR2， AR2：N， AR3：;-; FFT 运算处理部分;-SETC OVM ;使能溢出模式SETC SXM ;符号扩展使能SPM 1 ;PREG 寄存器的输出左移一位，;自动将两个 Q15 相乘后化为 Q15
9、的格式LACC *+，AR3ADD #1SACL *+，1，AR2 ;IW=2*(N+1);堆栈分布情况：ADDRESS/AR6/AR7/AR0/N/M/ID/IW/Y;ARP：AR2 ， AR2：M， AR3： Y， LAR AR5，*+ ;AR5=M，ARP ：AR2， AR2：ID， AR3：Y， AR5=MLOOP3 LAR AR6，#_input ;AR6：inputRi ，ARP：AR2， ;AR2：ID， AR3：Y， AR5=M， AR6：input LACC *，1SACL *+ ;ID=ID*2，ARP：AR2，AR2 ：IW，;AR3：Y， AR5=M， AR6：inpu
10、tLACC *，15SACH * ;IW=IW/2LACC *-，15，AR3SACH *+，AR2 ;C2=IW/2;堆栈分布情况：ADDRESS/AR6/AR7/AR0/N/M/ID/IW/C2/Y;ARP：AR2 ， AR2：ID， AR3：Y， AR5=M， AR6：inputLARAR0，* ;AR0=IDLOOP2LAR AR4，#_sintab ;AR4：sintab;ARP：AR2 ，AR0=ID，AR2：ID， AR3：Y，;AR4：sintab， AR5=M， AR6：inputLACC *+，15，AR3SACH *+，AR6 ;C1=ID/2=1;堆栈分布情况：ADDR
11、ESS/AR6/AR7/AR0/N/M/ID/IW/C2/C1/Y;ARP：AR6 ， AR0=ID， AR2：IW ， AR3：Y，;AR4：sintab， AR5=M， AR6：inputMAR *0+，AR4;ARP：AR4 ，AR0=ID，AR2：IW，AR3：Y，;AR4：sintab，AR5=M，AR6=AR6+IDRjLOOP1LACC #0LT *+，AR6 ;TREG=COSlkMPY*+，AR4 ;Rj* COSlk，ARP=AR4，AR4：SIN lk，AR6：IjLT *，AR6MPYA *-，AR3 ;ACC=ACC+Rj*COSlk， PREG=Ij*SINlk;A
12、RP=AR3， AR4：SIN lk， AR6： RjSPAC ;ACC=ACC-Ij*SINlkSACH *+，AR4 ;XT=Rj*COSlk-Ij*SINlk;堆栈分布情况：ADDRESS/AR6/AR7/AR0/N/M/ID/IW;/C2/C1/XT/Y;ARP：AR4 ， AR0=ID， AR2：IW ， AR3：Y，;AR4：SINX ， AR5=M， AR6：Q.XLACC #0LT *-，AR6MPY *+，AR4 ;Rj*SINlk，ARP=AR4，AR4：COS lk，AR6：IjLT *，AR6MPYA *-，AR3 ;ACC=ACC+Rj*SINlk， PREG=Ij*
13、COSlk;ARP=AR3， AR4：COS lk， AR6：RjAPACSACH *-，AR7 ;YT=Rj*SINlk+Ij*COSlk;堆栈分布情况：ADDRESS/AR6/AR7/AR0/N/M/ID/IW;/C2/C1/XT/YT;ARP：AR7 ， AR0=ID， AR2：IW ， AR3：XT ，;AR4：COSX ， AR5=M， AR6：Q.X LACC *，AR6BCND D2，NEQ ;当_nom 不为 0 时，需要在运算过程中进行归一化操作;-; 归一化处理程序部分;-D2MAR *0- ;AR6RiLAC *，15，AR3ADD *，15，AR6SACH *0+，AR
14、3 ;Ri =( Ri +XT)/2;ARP：AR3 ， AR0=ID， AR2：IW ， AR3：XT ，;AR4：COS lk， AR5=M， AR6：RjSUB *+，16，AR6SACH *+ ;Rj=Ri-XT， AR6：Ij， AR3：YT;ARP：AR6 ， AR0=ID， AR2：IW ， AR3：YT ，;AR4：COS lk， AR5=M， AR6：IjMAR *0- ;AR6：IiLACC *，15，AR3ADD *，15，AR6SACH *0+，AR3 ;Ii=(Ii+YT)/2， AR6：Ij;ARP：AR3 ， AR0=ID， AR2：IW ， AR3：YT ，;A
15、R4：COS lk， AR5=M， AR6：IjSUB *-，16，AR6SACH *+，0，AR2 ;Ij=Ii-YT;STACK：ADDRESS/AR6/AR7/AR0/N/M/ID/IW;/C2/C1/XT/YT;ARP：AR2 ，AR0=ID，AR2：IW，AR3：XT ，;AR4：COS lk， AR5=M，AR6：下一个 RjDLAR AR0，*-，AR4 ;AR0=IW;ARP=AR4， AR0=IW， AR2： ID，;AR3：XT， AR4：COS lk， AR5=M，AR6：下一个 RjMAR *0+，AR2 ;AR4=AR4+IW下一个 COSlkLAR AR0，* ;A
16、R0=IDADRK #3 ;AR2：C1;ARP=AR2， AR0=ID， AR2：C1， AR3：XT，;AR4： 下一个 COSlk， AR5=M， AR6：RjLACC *SUB #1SACL*- ;C1=C1-1;ARP=AR2， AR0=ID， AR2：C2， AR3：XT，;AR4：COS lk， AR5=M， AR6：RjBCND LOOP4，LEQ ;跳转至 LOOP4， IF C10MAR *-，AR4 ;AR2：IW;ARP=AR4， AR0=ID， AR2：IW， AR3：XT ，;AR4： 下一个 COSlk， AR5=M， AR6：RjB LOOP1 LOOP4LAC
17、C *SUB #1SACL *-，AR3MAR *-，AR2 ;ARP=AR2， AR0=ID， AR2：IW， AR3：C1，;AR4： 下一个 COSlk， AR5=M， AR6：RiMAR *-BCND LOOP2，GTMAR *，AR3MAR *-，AR5 ;ARP=AR5， AR0=ID， AR2：IW， AR3：C2，;AR4 下一个 COSlk， AR5=M， AR6：RiBANZ LOOP3，*-，AR2;-; 与 C 语言兼容的代码部分;-CLRC OVMSPM 0 MAR *，AR1SBRK #09LAR AR0，*-LAR AR7，*-LAR AR6，*-PSHD *RET程序说明：在进行FFT之前，需要做一些准备工作，对FFT应用到的数据进行反序排列。反序排列子程序部分用到间接寻址选项*BR0+，它提供了了位反转变址寻址的能力，可以对基-2FFT程序中数据点重排序以完成原位运算操作，利用这个特殊的指令可以很容易实现FFT操作。在FFT之前要对一些寄存器进行初始化。
