FFT算法的DSP实现.docx

1、FFT算法的DSP实现FFT 算法的 DSP 实现对于离散傅里叶变换（DFT）的数字计算，FFT是一种有效的方法。一般假定输入序列是 复数。当实际输入是实数时，利用对称性质可以使计算 DFT 非常有效。一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成 一个 N 点的复序列， 之后对复序列进行 N 点的 FFT 运算， 最后再由 N 点的复数输出拆散成 2N点的复数序列，这 2 N点的复数序列与原始的 2N点的实数输入序列的 DFT输出一致。使用这种方法， 在组合输入和拆散输出的操作中， FFT 运算量减半。 这样利用实数 FFT 算法来计算实输入序列的 DFT的速

2、度几乎是一般 FFT算法的两倍。下面用这种方法来实现 一个 256 点实数 FFT（2N=256 ）运算。1. 实数 FFT 运算序列的存储分配如何利用有限的 DSP 系统资源，合理的安排好算法使用的存储器是一个比较重要的问 题。本文中，程序代码安排在 0x3000 开始的存储器中，其中 0x30000x3080 存放中断向量 表，FFT程序使用的正弦表和余弦表数据 （.data段）安排在OxcOO开始的地方，变量（.bss段定义 ） 存放在 0x80 开始的地址中。另外，本文中 256 点实数 FFT 程序的数据缓冲位 Ox23OOOx23ff , FFT 变换后功率谱的计算结果存放在 Ox

3、22OOOx22ff 中。连续定位 .cmd 文件 程序如下：MEMORY PAGE O:IPROG: origin = Ox3O8O,len=Ox1F8OVECT: lorigin=Ox3OOO,len=Ox8OEPROG: origin=Ox38OOO,len=Ox8OOOPAGE 1:USERREGS: origin=Ox6O,len=Ox1cBIOSREGS: origin=Ox7c,len=Ox4IDATA: origin=Ox8O,len=OxB8OSECTIONSEDATA: origin=OxCOO,len=Ox14OO.vectors: VECT PAGE O.sysregs

4、:.trcinit:.gblinit: BIOSREGS PAGE 1 IPROG PAGE O IPROG PAGE O.bios:frt: IPROG PAGE O IPROG PAGE O.text: IPROG PAGE O.cinit: IPROG PAGE O.pinit: IPROG PAGE O.sysinit: IPROG PAGE O.data: .bss:.far:.const: EDATA PAGE 1 IDATA PAGE 1 IDATA PAGE 1 IDATA PAGE 1.switch: IDATA PAGE 1.sysmem: IDATA PAGE 1cio:

5、 IDATA PAGE 1.MEM$obj: IDATA PAGE 1.sysheap: IDATA PAGE 12.基2实数FFT运算的算法该算法主要分为以下四步进行：1）输入数据的组合和位排序首先，原始输入的2N=256个点的实数序列复制放到标记有 “ d_input_addr 的相邻单元，当成N=128点的复数序列dn，其中奇数地址是 dn实部，偶数地址是 dn的虚部，这个 过程叫做组合（n为序列变量，N为常量）。然后，把输入序列作位倒序，是为了在整个运 算最后的输出中得到的序列是自然顺序，复数序列经过位倒序，存储在数据处理缓冲其中， 标记为”_data ”如图一所示，输入实数序列为 a

6、n, n=0,1,2,3, , ,255。分离an成两个序列，如图二所 示，原始的输入序列是从地址 0x2300到0x23FF,其余德从0x2200到0x22FF的是经过为倒 序之后的组合序列：n=0，1, 2，3，,， 127。dn表示复合FFT的输入，rn表示实部，in表示虚部：dn=rn+j in按位倒序的方式存储 d n到数据处理缓冲中，如图二所示。2200h2200hR0=a02201h2201hI0=a12202h2202hR64=a1282203h2203hl64=a1292204h2204hR32=a642205h2205hI32=a652206h2206hR96=a19222

7、07h2207hI96=a1932208h2208hR16=a322209h2209hI16=a33220ah220ahR80=a160220bh.220bh.I80=a16122fehR127=a25422ffh22ffhl127=a2552300hA02300hA02301hA12301hA12302hA22302hA22303hA32303hA32304hA42304hA42305hA52305hA52306hA62306hA62307hA72307hA72308h:A82308hA82309hA92309hA9230ahA10230ahA10230bhA11230bhA11A25423

8、ffhA25523ffhA255程序设计中，在用 C54x进行位倒序组合时，使用位倒序寻址方式可以大大提高程序执 行的速度和使用存储器的效率。 在这种寻址方式中，AR0存放的整数N是FFT点数的一半，一个辅助寄存器指向一个数据存放的单元。 当使用位倒序寻址把 AR0加到辅助寄存器中时，地址以位倒序的方式产生，即进位是从左到右，而不是从右到左。例如，当 AR0=0x0060 ，AR2=0x0040时，通过指令：MAR AR2+0B就可以得到 AR2位倒序寻址后的值为 0x0010。下面是0x0060(1100000)与0x0040(1000000)以位倒序方式相加的过程：1100000+ 100

9、00000010000实现256点数据位倒序存储的具体程序段如下：:AR2( REORDERED_DATA):装入第一个位倒序数据指针注意，在上面的程序中，输入缓冲指针 AR3 (即ORIGINAL_INPUT )在操作时先加1再减1，是因为我们把输入数据相邻的两个字看成一个复数，在用寄存器间接寻址移动了一 个复数(两个字的数据)之后，对AR3进行位倒序寻址之前要把 AR3的值恢复到这个复数 的首字的地址，这样才能保证位倒序寻址的正确。2)N点复数FFT运算在数据处理器里进行 N点复数FFT运算。由于在FFT运算中要用到旋转因子 WN，它是一个复数。我们把它分为正弦和余弦部分，用 Q15格式将

10、它们存储在两个分离的表中。每个表中有128项，对应从01800。因为采用循环寻址地址来对表寻址， 128=27 28，因此每张表排队的开始地址就必须是 8个LSB位为0的地址。按照系统的存储器配置，把正弦表第一项sine_table”放在 0x0D00的位置，余弦表第一项 cos-table”放在 OxOEOO 的位置。根据公式N-1D(k)八 dnW k=0,1, ,N-1n=0利用蝶形结对dn进行N=128点复数FFT运算，其中nk (2-/N)nk 2WN =e cos( n k) js in( nk)N N所需的正弦值和余弦值分别以 Q15的格式存储于内存区以 0x0E00开始的余弦表

11、中。把128点的复数FFT分为七级来算， 形结，第二级是计算 4点的FFT蝶形结，然后是8点、16点、 结计算。最后所有的结果表示为0x0D00开始的正弦表和以第一级是计算2点的FFT蝶32点、64点、128点的蝶形DK=Fd n=Rk+jlk其中，Rk、 Ik分别是DK的实部和虚部。FFT完成以后，结果序列 DK就存储到数据处理缓冲器的上半部分，如图三所示，下半部分仍然保留原始的输入序列 an,这半部分将在第三步中被改写。这样原始的 a n序列的所有DFT的信息都在D(k)中了，第三步中需要做的就是把 D(k)变为最终的2N=256点复合序列，Ak=Fa(n)。2200hR02201hI02

12、202hR12203hI12204hR22205hI22206hR3 :2207hI32208hR42209hI4220ahR5220bhI522fehR12722ffhl1272300hA02301hA12302hA22303hA32304hA42305hA52306hA62307hA72308hA82309hA9 230ahA10230bhA1123fehA25423ffhA255注意，在实际的 DSP的编程中为了节约程序运行时间，提高代码的效率，往往要用到 并行程序指令。比如：ST B, *AR3II LD *AR3+, B并行指令的执行效果是，使原本分开要两个指令周期才能执行完的两条指

13、令在一个指 令周期中就能执行完。上述指令时将 B移位(ASM-16)所决定的位数，存于 AR3所指定的存储单元中，同时并行执行，将 AR3所指的单元中的值装入到累加器 B的高位中。由于指令的src和dst都是Acc、B，所以存入*AR3中的值是这条指令执行以前的值。这一步中，实现FFT计算的具体程序如下：Fft:;计算FFT的第一步，两点的 FFTasgAR1,GROUP_COUNTER;定义FFT计算的组指针asgAR2,PX；AR2为指向参加蝶形运算第一个；数据的指针asgAR3,QX；AR2为指向参加蝶形运算第二个；数据的指针asgAR4,WR;定义AR4为指向余弦表的指针asgAR5,

14、WT;定义AR5为指向正弦表的指针asgAR6,BUTTERFL Y_COUNTER;定义AR6为指向蝶形结的指针asgAR7,DATA_PROC_BUF；定义在第一步中的数据处理缓冲指针asgAR7,STAGE_COUNTER;定义剩下几步中的数据处理缓冲指针pshmst0pshmar0pshmbk;保存坏境变量SSBXSXM；开启符合扩展模式STM#K_ZERO_BK, BK;让 BK=0 使 *ARn+O%=*ARn+oLD#-1, ASM；为避免溢出在每一步输出时右移一位MVMM DATA_PROC_BUF, PX；PX 指向参加蝶形结运算的第一个数；的实部（ PR）LD*PX, 16

15、, A； AH: =PRSTM#fft_data+K_DATA_IDX_1,QX ；指向参加蝶形运算的第二个数的实;部(QR)STM#K_FFT_SIZE/2-1, BRC；设置块循环计数器RPTBD STAGELEND-1;语句重复执行的范围到地址; STAGELEND-1 处STM#K_DATA_IDX_1+1, AR0;延迟执行的两字节的指令（该指令;不重复执行）SUB*QX, 16, A, B; BH: =PR-QRADD*QX, 16, A; AH: =PR+QRSTHA, ASM, *PX+; PR: =(PR+QR)/2STB,*QX+; QR: =(PR-QR)/2II LD*

16、PX, A; AH: =PISUB*QX, 16, A, B; BH： =PI-QIADD* QX, 16, A; AH： =PI+QISTHA, ASM, *PX+0%; PI： =(PI+QI)/2ST B, *QX+0%; QI： =( PI-QI )/2I LD*PX, A; AH: =next PRstagelend:;计算 FFT 的第二步，四点的 FFTMVMM DATA_PROC_BUF, PX; PX 指向参加蝶形运算第一个数据的;实部 (PR)STM#fft_data+K_DATA_IDX_2, QX; QX 指向参加蝶形运算第二个数据的;实部 (QR)STM#K_FFT_

17、SIZE/4-1, BRC;设置块循环计数器LD*PX, 16, A; AH: =PRRPTBD STAGEL2END-1;语句重复执行的范围到地址; STAGEL1END-1 处STM#K_DA TA_IDX_2+1, AR0;初始化 AR0 以被循环寻址以下是第;一个蝶形;运算过程SUB*QX, 16, A, B; BH: =PR-QRADD*QX, 16, A; AH:= PR+QRSTHA, ASM, *PX+; PR:= (PR+QR)/2ST B, *QX+; QR:= (PR-QR)/2I LD*PX, A; AH:=PISUB*QX, 16, A, B; BH:=PI-QIAD

18、D*QX, 16, A; AH:=PI+QISTHA, ASM, *PX+; PI:= (PI+QI)/2STH B,ASM, *QX+; QI:= (PI-QI)/2;以下是第二个蝶形运算过程MAR*QX+;QX中的地址加一ADD*PX, *QX, A； AH:= PR+QISUB*PX, *QX-, B； BH: =PR-QISTHA, ASM, *PX+； QR := (PR+QI)/2SUB*PX, *QX, A； AH:=PI+QRSTB, *QX； QR := (PR-QI)/2II LD *QX, B； BH:=ORST A, *PX； PI:= (PI-QR)/2I ADD *

19、PX+0%, A； AH:= PI+QRSTA, *QX+0%； QIR := (PI+QR)/2I LD *PX, A； AH:=PRStage2end:；stage 3 thru stage logN-1: 从第三个蝶；形到第六个蝶形的过程如下STM#K_TWID_TBL_SIZE, BK； BK= 旋转因子表格的大小值ST#K_TWID_IDX_3, d_twid_idx；初始化旋转表格索引值STW#K_TWID_IDX_3, AR0； AR0= 旋转表格初始索引值STW#COS_TABLE, WR；初始化 WR 指针STW#SINE_TABLE, WI；初始化 WI 指针STW#K_L

20、OGN-2-1, STAGE_COUNTER；初始化步骤指针ST#K_FFT_SIZE/8-1, d_grps_cnt；初始化组指针STM#K_FL Y_COUNT_3-1, BUTTERFL Y_COUNTER；初始化蝶形指针ST#K_DA TA_IDX_3, d_data_idx；初始化输入数据的索引Stage:；以下是每一步的运算过程STM#fft_data, PX； PX 指向参加蝶形运算第一个数据；的实部（ PR ）LDd_data_idx, AADD*(PX), ASTLMA, QX； QX 指向参加蝶形运算第二个数据 ；的实部（ QR）MVDKd_grps_cnt, GROUP_

21、COUNTER； AR1 是组个数计数器group:；以下是每一组的运算过程MVMDBUTTERFLY_COUNTER, BRC；将每一组中的蝶形的个数装入 BRCRPTBDbutterflyend-1；重复执行至 butterflyend-1 处LD*WR, TMPY*QX+, A； A:=QR*W RI QX*QIMAC*WI+0%, *QX-, A； A:=QR*WR+QI*WIADD*PX, 16, A, B； B:=( QR*WR+QI*WI)+PR；I QX 指向 QRSTB, *PX； PR:=( (QR*WR+QI*WI)+PR)/2I SUB*PX+ B； B:=PR-( Q

22、R*WR+QI*WI)ST B, *QX； QR :=(PR-( QR*WR+QI*WI)/2I MPY*QX+, A； A:=QR*WI T=WTQX指向QIMAS *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)/2II SUB *PX, BLD *WR, TST B, *PX+I MPY *QX+, Abutterflyend:PSHM AR0MVDK d_data_idx, AR0MAR *PX+0MAR *QX+0BA

23、NZS group, *GROUP_COUNTER-POPM AR0MAR *QX-;I QX 指向 QR; B:=PI- (QR*WI-QI*WR); T:=WR; PI:=(PI- (QR*WI-QI*WR)/2;I PX 指向 PR; A:=QR*WR I QX 指向 QI;更新指针以准备下一组蝶形的运算;保存 AR0;ARO中装入在该步运算中每一组所；用的蝶形的数目;增加PX准备进行下一组的运算;增加QX准备进行下一组的运算；当组计数器减一后不等于零时，延; 迟跳转至 group 处;恢复ARO (一字节)；修改QX以适应下一组的运算更新指；针和其他索引数据以便进入下一个；步骤的运算L

24、D d_data_idx, ASUB #1, A, B; B=A-1STLM B, BUTTERFLY_COUNTERSTL A, 1, d_data_idxLD d_grps_cnt, ASTL A, -1, d_grps_cntLD d_twid_idx, ASTL A, -1, d_twid_idxBANZD stage, *STAGE_COUNTER_MVDK d_twid_idx, AR0POPM bkPOPM ar0POPM st0Fft_end:RET3)；修改蝶形个数计算器；下一步计算的数据索引翻倍；下一步计算的组数目减少一半；下一步计算的旋转因索引减少一半;人只0=旋转因子索

25、引(两字节)；恢复环境变量分离复数FFT的输出为奇部分和偶部分分离FFT输出为相关的四个序列： RP RM IP和IM,即偶实数、奇实数、偶虚数和奇虚数四部分，以便第四部形成最终结果。利用信号分析的理论我们把 D(K)通过下面的公式分为偶实数 RPK、奇实数RMK、偶虚数 IPK 和奇虚数 IMK:RPK=RPN-K=0.5*(RK+RN-K)RMK=-RMN-K=0.5*( RK-RN-K)IPK=IPN-K=0.5*(lK+lN-K)IMK=-IMN-K=0.5*( IK-IN-K)RPO=ROIP0=I0RM0=IM0=RMN/2=IMN/2=0RPN/2=RN/2IPN/2=IN/2R

26、MK图四显示了第三步完成以后存储器中的数据情况， RPK和IPK存储在上半部分,和IMK存储在下半部分。2200hRpO=RO2201hlpO=lO2202hRP12203hIP12204hRP22205hIP22206hRP32207hIP32208hRP42209hIP4220ahRP5220bhIP522fehRP12722ffhIP1272300hA02301hA12302hIM1272303hRM1272304hIM1262305hRM1262306hIM1252307hRM1252308hIM1242309hRM124230ahIM123230bhRM12323fehIM123ffhRM1这一过程的程序代码如下删除两个字unpackasg AR2, XP_K.asg AR3, XP_Nmi nusk.asg AR6, XM_K.asg AR7, XM_Nmi nuskSTM #fft_data+2, XP_K；AR2 指向 RK （ temp RPK）ST