专题八之三-通用数字信号处理方法DSP实现v2013

1、专题八专题八通用数字信号处理方法的通用数字信号处理方法的DSP实现实现DSP常见的几种信号处理算法：常见的几种信号处理算法：uDSP基本算术运算指令基本算术运算指令u除法运算除法运算u平方根运算平方根运算u级数展开级数展开u产生正弦波产生正弦波uFIR滤波器的实现滤波器的实现uFFT的实现的实现u自适应滤波（自适应滤波（LMS）的实现）的实现实现整数（定点）加法指令实现整数（定点）加法指令uC5400，C5500中提供了多条用于加法的指令，中提供了多条用于加法的指令，如如ADD。u加法指令中不同的情况，如用于加法指令中不同的情况，如用于无符号数的加无符号数的加法法运算，用于运算，用于带进位的加

2、法运算带进位的加法运算（如（如32位扩展位扩展精度加法），以及用于精度加法），以及用于立即数的加法立即数的加法。它们的。它们的具体加法指令有所不同。具体加法指令有所不同。u要实现要实现饱和运算加法饱和运算加法，需要设置标志位，需要设置标志位SATD。uC5500下还提供了同时完成加下还提供了同时完成加/减运算的指令减运算的指令ADDSUB等，以及条件执行加减，辅助寄存器等，以及条件执行加减，辅助寄存器加法等增强加法指令。加法等增强加法指令。实现实现C6000加法指令加法指令uC6000也提供了ADD加法指令。u无符号加法（ADDU），饱和运算（SADD）加法，立即数加法（ADDK）等都有专用指

3、令完成，等等u也提供同时完成加/减运算的ADDSUB、ADDSUB2指令。u也提供ADD2，ADD4的SIMD操作。uC674x以及C6700浮点加法指令使用ADDSP，ADDDP完成循环寻址的加法循环寻址的加法uC5500的循环寻址通过ARx辅助寄存器在指令中实现，可在任意间接寻址模式中使用。使用时可设置对应的ARnLC比特位（ST2_55）。也可在指令后加.CRuC6000的循环寻址可以通过指令ADDAB，ADDAD，ADDAH，ADDAW，并使用B4-B7/A4-A7来实现实现减法操作实现减法操作u减法运算指令与加法类似，毕竟加一个负数就减法运算指令与加法类似，毕竟加一个负数就等效于做一

4、个减法操作。等效于做一个减法操作。u减法指令中，减法指令中，SUBC为移位减，为移位减，DSP中的除法中的除法就是用该指令来实现的。就是用该指令来实现的。uC5400，C5500，C6000都提供都提供SUBC指令。指令。u参考教材的参考教材的5.1.1小节（小节（Page247），如何实现），如何实现整数除法以及求模运算。整数除法以及求模运算。u教材教材Page252程序示例程序示例C6X整数除法整数除法C程序。程序。实现实现16位定点整数除法位定点整数除法u在在C54X中没有提供专门的除法指令，一中没有提供专门的除法指令，一般有两种方法来完成除法。一种是用乘般有两种方法来完成除法。一种是用

5、乘法来代替，除以某个数相当于乘以其倒法来代替，除以某个数相当于乘以其倒数，所以先求出其倒数，然后相乘。这数，所以先求出其倒数，然后相乘。这种方法对于除以常数特别适用。另一种种方法对于除以常数特别适用。另一种方法是使用方法是使用SUBC指令，重复指令，重复16次减法完次减法完成除法运算。成除法运算。utemp1/temp2为例，说明如何使用为例，说明如何使用SUBC指令实现整数除法指令实现整数除法: ld temp1,T ；将被除数装入；将被除数装入T寄存器寄存器 mpy temp2,A ；除数与被除数相乘，结果放入；除数与被除数相乘，结果放入A寄存器寄存器 ld temp2,B ；将除数；将除

6、数temp2装入装入B寄存器的低寄存器的低16位位 abs B ；求绝对值；求绝对值 stl B,temp2 ；将；将B寄存器的低寄存器的低16位存回位存回temp2 ld temp1,B ；将被除数；将被除数temp1装入装入B寄存器的低寄存器的低16位位 abs B ；求绝对值；求绝对值 rpt #15 ；重复；重复SUBC指令指令16次次 subctemp2,b ；使用；使用SUBC指令完成除法运算指令完成除法运算 bcddiv_end,agt；延时跳转，先执行下面两条指令，；延时跳转，先执行下面两条指令， ; 然后判断然后判断A，若，若A0，则跳转到标号，则跳转到标号 ; div_en

7、d，结束除法运算，结束除法运算stl B,quot_i ; 将商（将商（B寄存器的低寄存器的低16位）存入变量位）存入变量quot_isth B,remain_i ；将余数（将余数（B寄存器的高寄存器的高16位）存入变量位）存入变量remain_ixor B ；若两数相乘的结果为负，则商也应为负。；若两数相乘的结果为负，则商也应为负。Sub quot_i,B ；将商反号；将商反号stl B,quot_i ；存回变量；存回变量quot_i中中div_end:实现乘法运算实现乘法运算u乘法指令在C5500，C6000中都有大量的指令来实现。u在C5000中提供了大量的乘法运算指令，16位相乘，其结

8、果都是32位，放在累加器中。u乘数在C5000的乘法指令很灵活，可以是T寄存器、立即数、存贮单元和累加器的高16位。u有符号与无符号数乘法是不同的指令。实现乘法运算实现乘法运算-C6000uC6000也有MPY乘法指令，但其乘法的乘数有16位，32位，其结果也有32位，64位的不同指令。u也提供有符号与无符号的区别指令。uC674x以及C67+提供MPYSP/MPYDP浮点数乘法指令。uC6000的乘法指令需要2个或更多的时钟周期，如MPYSP（4），MPYDP（10），等。具体请参考指令手册。小数乘法小数乘法u小数乘法与整数乘法主要差别在乘积多出一个符号比特，需要左移一位来消除uC5000下

9、使用状态寄存器的FRCT位来控制，当FRCT=1时进行小数乘法。uC6000使用独立指令来实现乘积左移，实现小数乘法，如SMPY，SMPY32，SMPY2，SMPYHL，等待16位定点整数乘法例子位定点整数乘法例子u在在C54X中，小数的乘法与整数乘法基本一致，中，小数的乘法与整数乘法基本一致，只是由于两个有符号的小数相乘，其结果的小只是由于两个有符号的小数相乘，其结果的小数点的位置在次高的后面，所以必须左移一位，数点的位置在次高的后面，所以必须左移一位，才能得到正确的结果。才能得到正确的结果。C54X中提供了一个状中提供了一个状态位态位FRCT，将其设置为，将其设置为1时，系统自动将乘积时，

10、系统自动将乘积结果左移移位。结果左移移位。rsbxFRCT；清；清FRCT标志，准备整数乘标志，准备整数乘ldtemp1,T；将变量；将变量temp1装入装入T寄存器寄存器mpytemp2,a；完成；完成temp2*temp1，结果放，结果放 ；入；入A寄存器（寄存器（32位）位）16位定点小数乘法例子位定点小数乘法例子u0ccdH（十进制的（十进制的0.1）x 0599aH(十进制的十进制的0.7），），两数相乘后两数相乘后B寄存器的内容为寄存器的内容为08f5f0a4H（十进（十进制的制的0.07000549323857）。）。u可以使用可以使用RND或使用或使用MPYR指令对低指令对低1

11、6位做四位做四舍五入的处理。舍五入的处理。Ssbx FRCT ；FRCT=1，准备小数乘法，准备小数乘法ldtemp1,16,a ；将变量；将变量temp1装入寄存器装入寄存器A的高的高16位位mpya temp2 ；完成；完成temp2乘寄存器乘寄存器A的高的高16位，结位，结 ；果在；果在B中，同时将中，同时将temp2装入装入T寄存器寄存器sthb,mpy_f ；将乘积结果的高；将乘积结果的高16位存入变量位存入变量mpy_f实现实现16定点小数除法定点小数除法u在在C54X中实现中实现16位的小数除法与前面的整数位的小数除法与前面的整数除法基本一致，也是使用循环的除法基本一致，也是使用

12、循环的SUBC指令来指令来完成。但有两点需要注意：完成。但有两点需要注意：u第一，小数除法的结果一定是小数（小于第一，小数除法的结果一定是小数（小于1），），所以被除数一定小于除数。这与整数除法正好所以被除数一定小于除数。这与整数除法正好相反。所以在执行相反。所以在执行SUBC指令前，应将被除数指令前，应将被除数装入装入A或或B寄存器的高寄存器的高16位，而不是低位，而不是低16位。位。其结果的格式与整数除法一样，其结果的格式与整数除法一样，A或或B寄存器寄存器的高的高16位为余数，低位为余数，低16位为商。位为商。u第二，与小数乘法一样，应考虑符号位对结果第二，与小数乘法一样，应考虑符号位对

13、结果小数点的影响。所以应对商右移一位。小数点的影响。所以应对商右移一位。乘累加运算乘累加运算-C5500uMAC/MAS是DSP的核心运算，C5500提供直接的乘累加或乘累减指令。uC5500在完成MAC/MAS同时，还可以直接完成数据的搬移，如MACMZuC5500提供MAC/MPY/MOV的并行执行能力，如MAC:MAC，MAC:MAS，MAC:MPY，MACMZ:MOV，等等乘累加运算乘累加运算-C6000uC62x没有直接的MAC指令，主要通过功能单元的并行能力来实现这个DSP核心操作！u但从C64开始，增加了DOTP2，DOTPN2等乘累加操作指令。u从C64+开始，增加了DDOTP

14、2等多次乘累加操作指令（需要pipeline配合）。u在C66x ，有DOTP4H等4个乘累加操作。特殊乘法指令特殊乘法指令uC5500的FIRSADD/FIRSSUB，LMS，SQA，SQDST，SQR等等uC64+的复数乘法：CMPYuC64+的Galois field乘法（伽罗瓦域伽罗瓦域）：GMPY，XORMPY,等。可用于加密处理的多项式运算uC6700的平方根指令：RSQRSP，等uC6700的倒数指令：RCPSP，RCPDP，等除法的实现除法的实现5.1.4小节（小节（Page257）除法的实现除法的实现uDSP中没有直接除法指令u5.1.1小节和5.1.2小节可以使用SUBC指

15、令实现整数和小数除法u更通常的做法是将除法改为乘以倒数。故除法问题转为倒数计算。u5.1.4小节给出了使用牛顿迭代法计算倒数。Newton-RhapsonuC6700提供倒数计算指令，但精度有限。牛顿迭代法计算倒数牛顿迭代法计算倒数u先定义函数，a为输入参数u则f(x)的根就是a的倒数。利用牛顿迭代法求去f(x)的根。u假设f (x)为f(x)的一阶导数，则u如果xn为f(x)的根，则迭代方程为：axxf1)()0(a2)(xxf.3 ,2 ,1 ,0,)()(1nnxfnxfnxnx.3 ,2 ,1 ,0),2(1nnxanxnxMatlab代码% 计算输入数据的倒数，其迭代的初始值应该保证

16、x2始终为正！a=input(Pls input the number=);x2=0.3;error=1;number=0;while error 0.00001 x1=x2; x2=x1*(2-a*x1); number=number+1; error=abs(x2-x1); disp(num2str(number), ,num2str(x2);enddisp(Result is ,num2str(x2);C6700的除法实现的除法实现u教材的Page259提供了单精度浮点除法运算的汇编程序实例代码。u同时给出了扩展精度的倒数计算方法。u教材Page260给出了C代码函数实现除法例子。u如果

17、是其它定点如果是其它定点DSP，可以参考前面的，可以参考前面的定点数运算，用牛顿迭代法计算倒数，定点数运算，用牛顿迭代法计算倒数，从而实现定点除法计算！从而实现定点除法计算！平方根的实现平方根的实现5.2小节（小节（Page267）平方根计算平方根计算u求平方根也是数字信号处理中常见的算术运算之一。利用上一节介绍的Newton-Rhapson迭代算法，可以很方便地完成平方根或者平方根倒数的运计算。u将平方根运算作如下变换：u由上式可以看出，求平方根和求平方根倒数的运算实际上是等价的。1(0)xaaaa平方根计算平方根计算u因此，定义函数u只要求得方程g(y) = 0的根，进行x = ay的运算

18、可得到a的平方根。u为了使用Newton-Rhapson迭代算法，计算g(y)的一阶导数并带入迭代方程：21( )(0)g yaay21 1.50.5 (1, 2, 3)y ny na y nn平方根计算平方根计算-C6700u若使用C67xx，可将y0 = RSQRSP(a)为初始值，即使用C67xx的求倒数的指令计算值作为初值。X0 = _rsqrsp(arg1) ; get initial TLU value to 8 bitsV = 0.5*arg1 ; generate term just onceX1 = X0*(1.5-V*X02); 1st iteration to 16-bi

19、tsX2 = X1*(1.5V*X12); 2nd iteration to 32-bitsAns = arg1*X2; generate square root* _rsqrsp为C6700的平方根计算指令，这里用作初值，提高精度。级数求和的实现级数求和的实现5.3小节小节（PLOY指令的应用）指令的应用）利用利用POLY指令计算泰勒级数指令计算泰勒级数u三角函数、对数、指数等超越函数都可以用级数展开，如常用的泰勒级数。u不同的级数展开有不同的收敛速度。u以指数 的展开为例：xe!.! 3! 2132nxxxxenx（实际一般取9项）利用利用POLY指令计算泰勒级数指令计算泰勒级数!.! 3

20、! 2132nxxxxenx)!.! 3! 21 (1121nxxxxn)! 8! 71(.! 41(! 31(! 211 (1xxxxx利用利用POLY指令计算泰勒级数指令计算泰勒级数uPOLY指令的含义: poly SmemRound(A(32-16) x T + B - ASmem Bu指数展开)! 8! 71(.! 41(! 31(! 211 (1xxxxx x - TSmem - 1/n!taylor: STM a9, AR3 ;AR3 points to coefficient in;Taylors equ. LD X, T;set up running environment u

21、sing LD *AR3+, 16, A ;powerful poly instruction on 54x DSP LD *AR3+, 16, B RPT #7;loop 8 times enough for audio app. POLY *AR3+ ;AH=fractional part in Q15;format利用利用POLY指令计算泰勒级数指令计算泰勒级数a9 .int 1, 7, 46, 273, 1365, 5461, 16384, 32767, 0, 0C5500计算级数计算级数uC5400提供了POLY指令实现级数迭代计算，但C5500，C6000都没有类似指令。u但C55

22、00可以通过并行指令来实现：u上面的并行指令完全可以替代POLY的功能uC6000则通过功能单元并行来实现级数运算则通过功能单元并行来实现级数运算MAC AC1,T0,AC2,AC1: MOV *AR1+ ar0 st #INSTR_B,*ar0(#OFF_INTR_3) st #fir,*ar0(#OFF_INTR_3+1) ; init A/D int vector ! ld #temp,dp ;- the following codes for MAC st #win_data,t_ar3 ; ar3 - 2000h data windows st #filter_coff,t_ar2

23、; ar2 - filter coff xor a,a ; clear a xor b,b ; clear b rsbx intm ; enable all int !g: idle 1 b g;-; These codes may be called by serial rev int ! all registers dont ; change !; When enter this subroutine, the AD data has been put in A !;- fir: ld #temp,dp stl a,temp ; a - AD data ! call low_pass_ma

24、c ld a,-16,b ld #0,dp stl b,2,TDXR ; sent result to DA ! rete;*; LOWPASS FILTER (use MAC); Input data - A, output data - A; used AR2 - coff,AR3 - data_buffer !;*low_pass_mac: pshm st1 pshm st0 pshm ar0 pshm bk mvdm #t_ar2,ar2 ; restore ar2 mvdm #t_ar3,ar3 ; restore ar3 stm #1,ar0 stm #N,bk ; set cir

25、cular addressing size stl a,*ar3+% rptz a,#(N-1) ; 0 - a, then repeat 34 times mac *ar2+0%,*ar3+0%,a ; done FIR filter, result in a mvmd ar3,#t_ar3 ; save ar3 mvmd ar2,#t_ar2 ; save ar2 popm bk popm ar0 popm st0 popm st1 retu程序、数据的存储器安排程序、数据的存储器安排u程序功能框图程序功能框图u相关外设的准备：相关外设的准备：DSP，AD/DA，u相关外设的软件设置：相关

26、外设的软件设置：McBSP串口初始串口初始化、化、AC01的初始化、的初始化、u硬件电路设计、调试硬件电路设计、调试u软件设计、调试软件设计、调试FIRFIR滤波器的工程实现滤波器的工程实现初始化串口初始化串口初始化初始化AC01 等待新数据？等待新数据？调用滤波程序调用滤波程序串口串口发送发送中中断服务程序断服务程序串口串口接受接受中中断服务程序断服务程序是llCCS图形工具显示图形工具显示FIR滤波效果滤波效果图形显示图形显示FIR输入输入/输出频谱输出频谱FIRS指令来实现指令来实现FIR滤波器滤波器u一种有限单位冲激响应呈现对中心点对称的一种有限单位冲激响应呈现对中心点对称的FIR滤波

27、器。长度为滤波器。长度为N 的线性相位的线性相位FIR滤波器的滤波器的输出表达式为：输出表达式为：12/0)1()()()(NkkNnxknxkhnyFIRS指令来实现指令来实现FIR滤波器滤波器uFIRS Xmem,Ymem,pmad含义：FIRS指令来实现指令来实现FIR滤波器滤波器 B +（A(32-16) x Pmad - B(Xmem + Ymem) APAR+uPmad寻址FIR滤波器系数，Xmem和Ymem分别指向窗口数据的上下两部分。16点点FIRS滤波数据存放滤波数据存放*ar2 *ar3new *ar2+16点点FIRS滤波数据存放滤波数据存放FIR *AR2+0%,*AR

28、3+0%#FIR_COFAR2-, AR3-=2*ar2 *ar3new *ar2+使用使用FIRS指令完成滤波指令完成滤波u利用利用FIRS指令指令, ,需要将输入数据缓冲分需要将输入数据缓冲分成两个成两个, ,大小为大小为N/2。u初始状态将初始状态将AR2指到缓冲区指到缓冲区1的顶部，的顶部，AR3指到缓冲区指到缓冲区2的底部。每次滤波之前的底部。每次滤波之前, ,应先将缓冲区应先将缓冲区1顶部的数据移到缓冲区顶部的数据移到缓冲区2的底部的底部, ,新来的一个样本存储到缓冲区新来的一个样本存储到缓冲区1中时中时, ,并对缓冲区并对缓冲区1 1指针指针AR2加加1(1(使用循使用循环寻址环

29、寻址) )。使用使用FIRS指令完成滤波指令完成滤波u处理器然后使用处理器然后使用FIRS指令进行乘加运算。指令进行乘加运算。当然当然, ,在使用在使用FIRS指令前指令前, ,需要预先计算一需要预先计算一次求和次求和, ,以初始化以初始化A 。u在在RPTZ重复指令和循环寻址的配合下重复指令和循环寻址的配合下, ,完完成成FIR滤波滤波。u滤波完成后，需要对两个数据缓冲的指针进滤波完成后，需要对两个数据缓冲的指针进行修正，以便对下一个点进行处理。将行修正，以便对下一个点进行处理。将Buffer1的指针减的指针减1和和Buffer2的指针减的指针减2，使，使他们指向各自缓冲的数据队列的最后。他

30、们指向各自缓冲的数据队列的最后。STM #1,AR0 ； AR0=1STM #（N/2）,BK ；BK=N/2,循环寻址循环寻址BUFFER大小为大小为NMVDD *ar2, *ar3 ；更新；更新Buffer2STLA, *ar2+% ； 更新滤波窗口中的采样数据更新滤波窗口中的采样数据ADD ADD * *ar2+0% , ar2+0% , * *ar3+0%ar3+0% ； 初始化初始化A ARPTZ B, #(N/2-1) ) ； 重复重复FIRSFIRS指令指令N/2次次, ,先将先将B B清零清零FIRS *ar2+0%, *ar3+0%,filter_coff+N/2 ；完成滤波

31、计算。注意完成滤波计算。注意FIR滤波系数存放滤波系数存放 ；在程序存贮区，；在程序存贮区，filter_cofffilter_coff为系数起始地址为系数起始地址MAR MAR * *ar2-%ar2-% ； 修改修改Buffer1Buffer1指针指针MAR MAR * *+ar3(-2)%+ar3(-2)% ； 修改修改Buffer2Buffer2指针指针使用带使用带FIRSFIRS指令的循环寻址模式实现指令的循环寻址模式实现FIR滤波器，程序片段滤波器，程序片段如下：如下：( (输入数据在输入数据在AL中中, ,滤波结果在滤波结果在B中中)uFIRSADD Acx,ACy,Cmem,X

32、mem,Ymem FIRSADD指令指令-C5500-C5500 ACy = ACy + (ACx * Cmem): ACx = (Xmem16) + (Ymem16)uFIRADD.CR可使得后面的Cmem，Xmem，Ymem间接寻址都是用循环寻址。u在满足C5500的并行条件时，可以与其它指令并行执行。C6000实现滤波器-乘累加loop: ldh.d1*A8+,A2 | ldh.d2*B9+,B3 nop4 mpy.m1x A2,B3,A4 nop add.l1A4,A6,A6 sub.l2B0,1,B0 b0 b.s1loop nop5What can you put in paral

33、lel?Load InstructionsMPY2x=+a1a0A0LDW .D1 *A4+,A0 x1x0B0| LDW .D2 *B4+,B0A2A3a1 * x1a0 * x0MPY2 A0,B0,A3:A2| MPYH.M2 A0,B0,B5+a1x1 + a3x3 .a0 x0 + a2x2 .ADD .L1 A2,A6,A6| ADD .L2 A3,B6,B6A6B6final sumADD .L1 A6,B6,A4A4+DOTP2 with LDDW=+a2a0A1:A0LDDW .D1 *A4+,A1:A0| LDDW .D2 *B4+,B1:B0A2B2a3*x3 + a2*

34、x2a1*x1 + a0*x0DOTP2 A0,B0,A2| DOTP2 A1,B1,B2+intermediate sumADD A2,A3,A3a1a3:A5x2x0B1:B0 x1x3:final sumADD A3,B3,A4A4+| ADD B2,B3,B3intermediate sumA3B3In Ch 8, well get all these instructions working in paralleld0*c0d1*c1d2*c2d3*c3+ yd4*c4d5*c5d6*c6d7*c7coefdataBlock Real FIRfor (i = 0; I ndata;

35、i+) sum = 0;for (j = 0; j ncoef; j+) sum = sum + (di+j * cj);yi = sum;loop Iterationi,j0,00,1d0c0d1c1d1c0d2c2d2c1d3c3d3c2.Block Real FIR Example (DDTOPL2 )for (i = 0; I ndata; i+) sum = 0;for (j = 0; j ncoef; j+) sum = sum + (di+j * cj);yi = sum;loop Iterationi,j0,00,1d0c0+d1c1d1c0+d2c1d2c2d3c3d3c2.

36、uFour 16x16 multipliesuIn each .M unit every cycle-adds up to 8 MACs/cycle, or8000 MMACSuBottom Line: Two loop iterations for the price of oneDDOTPL2 d3d2:d1d0, c1c0, sum1:sum0loop Iterationi,j0,00,10,20,30,40,5d0c0+d1c1d1c0+d2c1d2c2d2c0d3c3d3c2d3c1d3c0d4c4d4c3d4c2+d5c3d4c1d4c0d5c5d5c4d5c2+d6c3d5c1d

37、6c0d6c6d6c5d6c4d6c2d6c1d7c7d7c6d7c5d7c4d7c3d6c2d8c7d8c6d8c5d8c4d6c3 DDOTPL2.M1 d3d2:d1d0, c1c0, sum1:sum0| DDOTPL2.M2 d7d6:d5d4, c3c2, sum3:sum2Parallel DDOTPL2sReal Block FIR A_j SPLOOPD 4| MVC .S2 B_i0, ILC ;set ILC | ADD .L2 B_i0, 1, B_i0 ;T/4 | ADDAH .D2 B_DLYaddr, nCoefs-1+4, B_DLYOUTaddr | MVK

38、 .S1 nCoefs/4-3, A_T| ADDAB .D1 A_DLYaddr, 8, A_DLYaddr*- stage A -* LDDW .D1T2 *+A_DLYaddr, B_d7d6:B_d5d4 ; 1,1 | LDDW .D2T1 *+B_DLYaddr, A_d3d2:A_d1d0 ; 1,1 LDDW .D2T2 *+B_COEFaddr, B_c3c2:B_c1c0 ; 2,1 | LDDW .D1T1 *+A_COEFaddr, A_c3c2:A_c1c0 ; 2,1 SPMASK| LDDW .D1T1 *A_DLYaddr1, A_dbda:A_d9d8 ; 3

39、,1 | LDDW .D2T2 *B_INaddr+, B_TEMP1:B_TEMP0 ;ld 1st input | MVC .S2 B_i0, RILC ;set RILC SPMASK| LDDW .D2T1 *B_INaddr+, A_TEMP1:A_TEMP0 ;ld 2nd input | ZERO .L1 A_st ;clear st flag | ADDAB .D1 DP, outputs+8, A_OUTaddr| MVK .S2 nCoefs/4-1, B_TC| MVK .S1 nCoefs/4-1, A_TC*- stage B -* SPMASK| SUB .L2X

40、A_OUTaddr, 8, B_OUTaddr NOP 1 DMV .S2X B_d5d4, A_d3d2, B_d5d4_:B_d3d2_ ; 7,1 | DDOTPL2 .M1 A_d3d2:A_d1d0,A_c1c0, A_pb:A_pf ; 7,1 | DDOTPL2 .M2 B_d7d6:B_d5d4,B_c3c2, B_p2:B_p6 ; 7,1 SPMASK| DMV .S1X A_d9d8, B_d7d6, A_d9d8_:A_d7d6_ ; 8,1 | DDOTPL2 .M2 B_d5d4_:B_d3d2_,B_c3c2, B_pa:B_pe ; 8,1 | DDOTPL2

41、.M1 A_dbda:A_d9d8,A_c3c2, A_p0:A_p4 ; 8,1 | STNDW .D2T2 B_TEMP1:B_TEMP0, *B_DLYOUTaddr+ ;st 1st inputNew double-throughput 16-bit mpy instructionsABCABCABCDDDABCD4ABCABCABCDDDABCDO1O2O3O4O5e1e2,p1e3,p2e4,p3p4DDOTP4DDOTP4 (.unit) src1, src2, dst_o:dst_eBlock Real FIRuC64x ImplementationDOTP2 d1d0, c1

42、c0, s0 ;d1*c1 + d0*c0DOTP2 d2d1, c1c0, s1 ;d2*c1 + d1*c0uC64x Plus ImplementationDDOTPL2 d3d2:d1d0, c1c0, s1:s0uReduces .M requirement by halfC64x: 194 cycles, 624 bytes (N=40, T=16)C64x+:126 cycles, 496 bytes (N=40, T=16) N = Length of example block, T = data sizeCMPY.Complex Multiply (CMPY)A0r1i1x

43、xA1r2i2=CMPY A0, A1, A3:A2r1*r2 - i1*i2:i1*r2 + r1*i2 32-bits 32-bits uFour 16x16 multiplies per .M unituUsing two CMPYs, a total of eight 16x16 multiplies per cycleuFloating-point version (CMPYSP) uses: 64-bit inputs (register pair)128-bit packed products (register quad)You then need to add/subtrac

44、t the products to get the final result single .M unitComplex MPYuCMPY (.M) src1, src2, dst_o:dst_eFour 16-bit inputs (real0,imag0, real1,imag1)Two 32-bit outputs (real, imag)uCMPYR (.M) src1, src2, dstFour 16-bit inputs (real0,imag0, real1,imag1)One packed 32-bit output (16 bit real, 16 bit imag)Rou

45、nded by adding 2*15 uCMPYR1 (.M) src1, src2, dstFour 16-bit inputs (real0,imag0, real1,imag1)One packed 32-bit output (real, imag)Rounded by adding 2*14 plus left shiftBlock Complex FIRuC64x ImplementationDOTP2dre_dim, cim_cre, pimDOTPN2 dre_dim, cre_cim, preuC64x Plus ImplementationCMPY dre_dim, cr

46、e_cim, pre:pimuReduces .M requirement by halfC64x: 674 cycles, 572 bytes (N=40, T=16)C64x+: 344 cycles, 460 bytes (N=40, T=16)N = Length of example block, T = data sizeDSPLIBuOptimized DSP Function Library for C programmers using C62x/C67x and C64x devicesuThese routines are typically used in comput

47、ationally intensive real-time applications where optimal execution speed is critical. uBy using these routines, you can achieve execution speeds considerably faster than equivalent code written in standard ANSI C language. And these ready-to-use functions can significantly shorten your development t

48、ime.uThe DSP library features:C-callableHand-coded assembly-optimizedTested against C model and existing run-time-support functionsAdaptive filteringMathDSP_firlms2DSP_dotp_sqrCorrelationDSP_dotprodDSP_autocorDSP_maxvalFFTDSP_maxidxDSP_bitrev_cplxDSP_minvalDSP_radix 2DSP_mul32DSP_r4fftDSP_neg32DSP_f

49、ftDSP_recip16DSP_fft16x16rDSP_vecsumsqDSP_fft16x16tDSP_w_vecDSP_fft16x32MatrixDSP_fft32x32DSP_mat_mulDSP_fft32x32sDSP_mat_transDSP_ifft16x32MiscellaneousDSP_ifft32x32DSP_bexpFilters & convolutionDSP_blk_eswap16DSP_fir_cplxDSP_blk_eswap32DSP_fir_genDSP_blk_eswap64DSP_fir_r4DSP_blk_moveDSP_fir_r8D

50、SP_fltoq15DSP_fir_symDSP_minerrorDSP_iirDSP_q15toflTechnical TrainingOrganizationT TOFFT实现实现教材教材7章章(7.3小节及小节及7.4小节小节)FFT是数字信号处理中重要的工是数字信号处理中重要的工具具uFFT是一种高效实现离散付氏变换的算是一种高效实现离散付氏变换的算法。法。(8.2.2小节介绍了小节介绍了 Goertzel算法）算法） u离散付氏变换的目的是把信号由时域变离散付氏变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信换到频域，从而可以在频域分析处理信息，得到的结果再由付氏逆变换到

51、时域。息，得到的结果再由付氏逆变换到时域。uDFT的定义为：的定义为： 1,.,1 , 0)2(10NkenxkXnkNjNnDFT的定义的定义u可以方便的把它改写为如下形式：可以方便的把它改写为如下形式：uWN（旋转因子）的周期性是（旋转因子）的周期性是DFT的关键的关键性质之一。为了强调起见，常用表达式性质之一。为了强调起见，常用表达式WN取代取代W以便明确其周期是以便明确其周期是N。 1,.,1 , 010NkWnxkXnkNNn)2sin()2cos()2(nkNjnkNeWnkNjnkNFFT是是DFT的快速算法的快速算法u由由DFT的定义可以看出，在的定义可以看出，在xn为复数序为

52、复数序列的情况下，完全直接运算列的情况下，完全直接运算N点点DFT需要需要（N-1）2次复数乘法和次复数乘法和N（N-1）次加法。）次加法。uFFT的基本思想在于，将原有的的基本思想在于，将原有的N点序列点序列序列分成两个较短的序列，这些序列的序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列可以很简单的组合起来得到原序列的的DFT。FFT是是DFT的快速算法的快速算法u例如，若例如，若N为偶数，将原有的为偶数，将原有的N点序列分成点序列分成两个（两个（N/2）点序列，那么计算）点序列，那么计算N点点DFT将将只需要约只需要约(N/2)2 2 2=N2/2次复数乘法。即比次

53、复数乘法。即比直接计算少作一半乘法。直接计算少作一半乘法。u该处理方法可以反复使用，即（该处理方法可以反复使用，即（N/2）点的）点的DFT计算也可以化成两个（计算也可以化成两个（N/4）点的）点的DFT（假定（假定N/2为偶数），从而又少作一半的乘为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后法。这样一级一级的划分下去一直到最后就划分成两点的就划分成两点的FFT运算的情况。运算的情况。FFT是是DFT的快速算法的快速算法u比如，一个比如，一个N = 8点的点的FFT运算按照这种运算按照这种方法来计算方法来计算FFT可以用下面的流程图来可以用下面的流程图来表示：表示：W0W0

54、W2W0W2W0W1W2W3x(0)x(4)x(2)x(6)x(1)x(5)x(3)x(7)X(7)X(6)X(5)X(4)X(3)X(2)X(1)X(0)W0W0W0实数实数FFT运算运算u对于离散傅立叶变换（对于离散傅立叶变换（DFT）的数字计）的数字计算，算，FFT是一种有效的方法。一般假定是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，输入序列是复数。当实际输入是实数时，利用对称性质可以使计算利用对称性质可以使计算DFT非常有效。非常有效。u一个优化的实数一个优化的实数FFT算法是一个组合以算法是一个组合以后的算法。原始的后的算法。原始的2N个点的实输入序列个点的实输入序

55、列组合成一个组合成一个N点的复序列，之后对复序列点的复序列，之后对复序列进行进行N点的点的FFT运算，最后再由运算，最后再由N点的复点的复数输出拆散成数输出拆散成2N点的复数序列。点的复数序列。实数实数FFT运算运算u使用这种方法，在组合输入和拆散输出使用这种方法，在组合输入和拆散输出的操作中，的操作中，FFT运算量减半。运算量减半。u这样利用实数这样利用实数FFT算法来计算实输入序算法来计算实输入序列的列的DFT的速度几乎是一般复的速度几乎是一般复FFT算法算法的两倍。的两倍。u本实验就用这种方法实现了一个本实验就用这种方法实现了一个256点实点实数数FFT（2N = 256）运算。）运算。

56、实数实数FFT运算运算序列序列的存的存储分储分配配程序存储区实数FFT程序1800h1FFFhFFT程序存储空间数据存储区0F00h0FFFh堆栈test_val0C00h对输入数据进行符号扩展的测试值temp_ar50C01h输出缓冲区的指针temp_ar30C02h输入缓冲区的指针input_cnt0C03h输入数据计数器output_cnt0C04h输出数据计数器d_grps_cnt0C05h组指针d_twid_idx0C06h旋转因子指针d_data_idx0C07h数据处理缓冲指针sine_table0D00h0DFFh正弦表cos_table0E00h0EFFh余弦表fft_da

57、ta2200h22FFh数据处理缓冲同时又是功率谱输出缓冲d_input_addr2300h23FFh数据输入缓冲基二实数基二实数FFT运算的算法运算的算法u第一步，输入数据的组合和位倒序：（第一步，输入数据的组合和位倒序：（1）把输入序列作位倒序，是为了在整个运把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序。算最后的输出中得到的序列是自然顺序。（2）原始的输入的）原始的输入的2N = 256个点的实数个点的实数序列复制放到标记有序列复制放到标记有“d_input_addr”的的相邻单元，当成相邻单元，当成N = 128点的复数序列点的复数序列dn。其中，奇数地址是。其中

58、，奇数地址是dn的实部，偶的实部，偶数地址是数地址是dn的虚部。复数序列经过位的虚部。复数序列经过位倒序，存储在数据处理缓冲器中倒序，存储在数据处理缓冲器中 “fft-data”。2206h2207h2208h2209h220Ah220Bh22FFh2300ha02301ha12302ha22303ha32304ha42305ha523FFha2552200h2201h2202h2203h2204h2205h2306ha62307ha72308ha82309ha9230Aha10230Bha112206hr96=a1922207hi96=a1932208hr16=a322209hi16=a33

59、220Ahr80=a160220Bhi80=a16122FFhi127=a2552300ha02301ha12302ha22303ha32304ha42305ha523FFha2552200hr0=a02201hi0=a12202hr64=a1282203hi64=a1292204hr32=a642205hi32=a652306ha62307ha72308ha82309ha9230Aha10230Bha1122FEhr127=a25423FEha254bit_rev: STM #d_input_addr,ar3 ；在；在AR3中放入输入地址中放入输入地址 STM #fft_data,ar7 ；

60、在；在AR7中放入处理后输出的地址中放入处理后输出的地址 MVMM DATA_PROC_BUF,ar2 ; AR2中装入第一个位倒序中装入第一个位倒序 ; 数据指针数据指针 STM #K_FFT_SIZE-1,BRC STM #K_FFT_SIZE,ar0 ; AR0=输入数据数目的一半输入数据数目的一半(128) RPTB bit_rev_end MVDD *ar3+,*ar7+; 将原始输入缓冲中的数据放入到将原始输入缓冲中的数据放入到 ; 位倒序缓冲中去之后输入缓冲位倒序缓冲中去之后输入缓冲 ; (AR3)指针加指针加1, 位倒序缓冲位倒序缓冲(AR2)指指 ; 针也加针也加1 MVDD *ar