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

专题八

通用数字信号处理方法的DSP实现DSP常见的几种信号处理算法：DSP基本算术运算指令除法运算平方根运算级数展开产生正弦波FIR滤波器的实现FFT的实现自适应滤波（LMS）的实现实现整数（定点）加法指令C5400，C5500中提供了多条用于加法的指令，如ADD。加法指令中不同的情况，如用于无符号数的加法运算，用于带进位的加法运算（如32位扩展精度加法），以及用于立即数的加法。它们的具体加法指令有所不同。要实现饱和运算加法，需要设置标志位SATD。C5500下还提供了同时完成加/减运算的指令ADDSUB等，以及条件执行加减，辅助寄存器加法等增强加法指令。实现C6000加法指令C6000也提供了ADD加法指令。无符号加法（ADDU），饱和运算（SADD）加法，立即数加法（ADDK）等都有专用指令完成，等等也提供同时完成加/减运算的ADDSUB、ADDSUB2指令。也提供ADD2，ADD4的SIMD操作。C674x以及C6700浮点加法指令使用ADDSP，ADDDP完成循环寻址的加法C5500的循环寻址通过ARx辅助寄存器在指令中实现，可在任意间接寻址模式中使用。使用时可设置对应的ARnLC比特位（ST2_55）。也可在指令后加.CRC6000的循环寻址可以通过指令ADDAB，ADDAD，ADDAH，ADDAW，并使用B4-B7/A4-A7来实现实现减法操作减法运算指令与加法类似，毕竟加一个负数就等效于做一个减法操作。减法指令中，SUBC为移位减，DSP中的除法就是用该指令来实现的。C5400，C5500，C6000都提供SUBC指令。参考教材的5.1.1小节（Page247），如何实现整数除法以及求模运算。教材Page252程序示例C6X整数除法C程序。实现16位定点整数除法在C54X中没有提供专门的除法指令，一般有两种方法来完成除法。一种是用乘法来代替，除以某个数相当于乘以其倒数，所以先求出其倒数，然后相乘。这种方法对于除以常数特别适用。另一种方法是使用SUBC指令，重复16次减法完成除法运算。temp1/temp2为例，说明如何使用SUBC指令实现整数除法:

ldtemp1,T

；将被除数装入T寄存器

mpytemp2,A

；除数与被除数相乘，结果放入A寄存器

ldtemp2,B

；将除数temp2装入B寄存器的低16位

absB

；求绝对值

stlB,temp2

；将B寄存器的低16位存回temp2

ldtemp1,B

；将被除数temp1装入B寄存器的低16位

absB

；求绝对值

rpt#15

；重复SUBC指令16次

subc temp2,b

；使用SUBC指令完成除法运算

bcd div_end,agt

；延时跳转，先执行下面两条指令，

;然后判断A，若A>0，则跳转到标号

;div_end，结束除法运算stlB,quot_i ;将商（B寄存器的低16位）存入变量quot_isthB,remain_i

；将余数（B寄存器的高16位）存入变量remain_ixorB

；若两数相乘的结果为负，则商也应为负。Subquot_i,B

；将商反号stlB,quot_i

；存回变量quot_i中div_end:

实现乘法运算乘法指令在C5500，C6000中都有大量的指令来实现。在C5000中提供了大量的乘法运算指令，16位相乘，其结果都是32位，放在累加器中。乘数在C5000的乘法指令很灵活，可以是T寄存器、立即数、存贮单元和累加器的高16位。有符号与无符号数乘法是不同的指令。实现乘法运算-C6000C6000也有MPY乘法指令，但其乘法的乘数有16位，32位，其结果也有32位，64位的不同指令。也提供有符号与无符号的区别指令。C674x以及C67+提供MPYSP/MPYDP浮点数乘法指令。C6000的乘法指令需要2个或更多的时钟周期，如MPYSP（4），MPYDP（10），等。具体请参考指令手册。小数乘法小数乘法与整数乘法主要差别在乘积多出一个符号比特，需要左移一位来消除C5000下使用状态寄存器的FRCT位来控制，当FRCT=1时进行小数乘法。C6000使用独立指令来实现乘积左移，实现小数乘法，如SMPY，SMPY32，SMPY2，SMPYHL，等待16位定点整数乘法例子在C54X中，小数的乘法与整数乘法基本一致，只是由于两个有符号的小数相乘，其结果的小数点的位置在次高的后面，所以必须左移一位，才能得到正确的结果。C54X中提供了一个状态位FRCT，将其设置为1时，系统自动将乘积结果左移移位。rsbx FRCT

；清FRCT标志，准备整数乘ld temp1,T

；将变量temp1装入T寄存器mpy temp2,a

；完成temp2*temp1，结果放；入A寄存器（32位）16位定点小数乘法例子0ccdH（十进制的0.1）x0599aH(十进制的0.7），两数相乘后B寄存器的内容为08f5f0a4H（十进制的0.07000549323857）。可以使用RND或使用MPYR指令对低16位做四舍五入的处理。SsbxFRCT

；FRCT=1，准备小数乘法ld temp1,16,a

；将变量temp1装入寄存器A的高16位mpya temp2

；完成temp2乘寄存器A的高16位，结；果在B中，同时将temp2装入T寄存器sth b,mpy_f

；将乘积结果的高16位存入变量mpy_f实现16定点小数除法在C54X中实现16位的小数除法与前面的整数除法基本一致，也是使用循环的SUBC指令来完成。但有两点需要注意：第一，小数除法的结果一定是小数（小于1），所以被除数一定小于除数。这与整数除法正好相反。所以在执行SUBC指令前，应将被除数装入A或B寄存器的高16位，而不是低16位。其结果的格式与整数除法一样，A或B寄存器的高16位为余数，低16位为商。第二，与小数乘法一样，应考虑符号位对结果小数点的影响。所以应对商右移一位。乘累加运算-C5500MAC/MAS是DSP的核心运算，C5500提供直接的乘累加或乘累减指令。C5500在完成MAC/MAS同时，还可以直接完成数据的搬移，如MACMZC5500提供MAC/MPY/MOV的并行执行能力，如MAC::MAC，MAC::MAS，MAC::MPY，MACMZ::MOV，等等乘累加运算-C6000C62x没有直接的MAC指令，主要通过功能单元的并行能力来实现这个DSP核心操作！但从C64开始，增加了DOTP2，DOTPN2等乘累加操作指令。从C64+开始，增加了DDOTP2等多次乘累加操作指令（需要pipeline配合）。在C66x，有DOTP4H等4个乘累加操作。特殊乘法指令C5500的FIRSADD/FIRSSUB，LMS，SQA，SQDST，SQR等等C64+的复数乘法：CMPYC64+的Galoisfield乘法（伽罗瓦域）：GMPY，XORMPY,等。可用于加密处理的多项式运算C6700的平方根指令：RSQRSP，等C6700的倒数指令：RCPSP，RCPDP，等除法的实现5.1.4小节（Page257）除法的实现DSP中没有直接除法指令5.1.1小节和5.1.2小节可以使用SUBC指令实现整数和小数除法更通常的做法是将除法改为乘以倒数。故除法问题转为倒数计算。5.1.4小节给出了使用牛顿迭代法计算倒数。Newton-RhapsonC6700提供倒数计算指令，但精度有限。牛顿迭代法计算倒数先定义函数，a为输入参数则f(x)的根就是a的倒数。利用牛顿迭代法求去f(x)的根。假设f’(x)为f(x)的一阶导数，则如果x[n]为f(x)的根，则迭代方程为：Matlab代码%计算输入数据的倒数，其迭代的初始值应该保证x2始终为正！a=input('Plsinputthenumber=');x2=0.3;error=1;number=0;whileerror>0.00001x1=x2;x2=x1*(2-a*x1);number=number+1;error=abs(x2-x1);disp([num2str(number),'',num2str(x2)]);enddisp(['Resultis',num2str(x2)]);C6700的除法实现教材的Page259提供了单精度浮点除法运算的汇编程序实例代码。同时给出了扩展精度的倒数计算方法。教材Page260给出了C代码函数实现除法例子。如果是其它定点DSP，可以参考前面的定点数运算，用牛顿迭代法计算倒数，从而实现定点除法计算！平方根的实现5.2小节（Page267）平方根计算求平方根也是数字信号处理中常见的算术运算之一。利用上一节介绍的Newton-Rhapson迭代算法，可以很方便地完成平方根或者平方根倒数的运计算。将平方根运算作如下变换：由上式可以看出，求平方根和求平方根倒数的运算实际上是等价的。平方根计算因此，定义函数只要求得方程g(y)=0的根，进行x=ay的运算可得到a的平方根。为了使用Newton-Rhapson迭代算法，计算g(y)的一阶导数并带入迭代方程：平方根计算-C6700若使用C67xx，可将y[0]=RSQRSP(a)为初始值，即使用C67xx的求倒数的指令计算值作为初值。X0=_rsqrsp(arg1) ;getinitialTLUvalueto8bitsV=0.5*arg1 ;generatetermjustonceX1=X0*(1.5-V*X02) ;1stiterationto16-bitsX2=X1*(1.5–V*X12) ;2nditerationto32-bitsAns=arg1*X2 ;generatesquareroot*_rsqrsp为C6700的平方根计算指令，这里用作初值，提高精度。级数求和的实现5.3小节（PLOY指令的应用）利用POLY指令计算泰勒级数三角函数、对数、指数等超越函数都可以用级数展开，如常用的泰勒级数。不同的级数展开有不同的收敛速度。以指数的展开为例：（实际一般取9项）利用POLY指令计算泰勒级数利用POLY指令计算泰勒级数POLY指令的含义:polySmemRound(A(32-16))xT+B->ASmem<16->B指数展开

x->TSmem->1/n!taylor:STMa9,AR3 ;AR3pointstocoefficientin ;Taylor’sequ.LDX,T ;setuprunningenvironmentusingLD*AR3+,16,A;powerfulpolyinstructionon54xDSPLD*AR3+,16,B

RPT#7 ;loop8timesenoughforaudioapp.POLY*AR3+ ;AH=fractionalpartinQ15 ;format利用POLY指令计算泰勒级数1,7,46,273,1365,5461,16384,32767,0,0C5500计算级数C5400提供了POLY指令实现级数迭代计算，但C5500，C6000都没有类似指令。但C5500可以通过并行指令来实现：上面的并行指令完全可以替代POLY的功能C6000则通过功能单元并行来实现级数运算MACAC1,T0,AC2,AC1::MOV*AR1+<<#16,AC2正弦波的产生5.4小节（Page278）查表法产生正弦波可以构造一个查找表，直接给出正弦函数的值。需要建立一个±/2范围输入所对应的正弦(余弦)函数的值。不妨将表的大小初设为257项,表中的第一个值对应于0°，最后一个值对应于180°，或者说。这样，表中相邻两点之间的间隔为:180/256=0.7031250°正弦函数和余弦函数之间只有一个90°的相移。查表法产生正弦波表格中的第1项是sin(0°)的值，第2项是sin(0.7031250°)的值，第3项是sin(1.406250°)的值，依此类推，最后一项是sin(180°)的值。利用该表所能构造的波形的相位步进，通常是0.703125°的整数倍。记录波形的初始相位和当前相位，就可以根据相位值去查找正弦表中相应的位置，得到当前输出点的幅度值。查表法产生正弦波例如，在语音处理中，如果需要产生的正弦波的频率为10kHz，采样频率为44.1kHz，对应的相位步进为：根据输出信号的初始相位确定第一个点，则下一个输出点的相位为81.6°，再下一个输出点的相位为163.2°，依此类推。得到这些相位值后，就可以从查找表中的对应位置，读取当前输出点的幅度。需要产生余弦波时，将相位值加上/2（90°）数字振荡器产生正弦数字振荡器的本质是，使用一个IIR（InfiniteImpulseResponse）滤波器，将一对共轭极点放在单位圆上，来产生正弦或余弦函数。利用正弦波sinx的指数形式（即系统的冲击相应）：可以得到正弦序列x[n]的z变换其中，A=2cos(T)，B=-1，C=sin(T)数字振荡器产生正弦转换为二阶差分方程该系统的结构框图：通过递推可以逐点计算系统输出，即正弦信号数字振荡器产生正弦例如，设该振荡器的频率为

，采

样频率为

，则

系数计算为：教材的5.4.2小节给出DSP实现的代码例子

FIR滤波器

教材6章（6.3节）

FIR滤波器的结构y(n)=h(0)x(n)+h(1)x(n-1）

+...+h(N-1)x[n-(N-1)]DSP实现FIR滤波器可以使用数字滤波器辅助设计软件包或自行计算FIR滤波器的系数。本实验例子中使用的是一个34阶的对称结构的FIR低通滤波器，其采样频率Fs为25KHZ，通带截止频率

1.5KHZ，阻带截止频率为3KHZ，阻带衰减为-40dB。(系数使用DFDP4自动生成）FIR滤波器的数据存储方式FIR滤波器的DSP实现为了能正确使用循环寻址，数据缓冲区和系数的开始地址必须正确指定。滤波系数指针初始化时指向h(N-1)，经过一次FIR滤波计算后，在循环寻址的作用下，仍然指向h(N-1)。而数据缓冲区指针指向的是需要更新的数据，如x(n)。在写入新数据并完成FIR运算后，该指针指向x(n-(N-1))。使用循环寻址可以方便地完成滤波窗口数据的自动更新.FIR滤波器的DSP实现使用带MAC指令的循环寻址模式实现FIR滤波器，程序片段如下：(输入数据在AL中,滤波结果在AH中)STM #1,AR0 ;AR0=1STM #N,BK ;BK=N,循环寻址BUFFER大小为NSTL A,*FIR_DATA_P+%;更新滤波窗口中的采样数据RPTZA,#(N-1)

；重复MAC指令N次,先将A清零MAC*FIR_DATA_P+0%,*FIR_COFF_P+0%,A

;完成滤波计算。注意FIR滤波系数存

；放在数据存储区

.title"exampleforFIRfilter!".mmregs.globalmainstartOFF_INTR_3.set0Chk_win_size.set34;34tapsFIRLPFilterwin_data.set2000h;FIRdatawindow.datafilter_coff.word8ch;low_pass,1.5kHz(3kHz),34h(N-1).word59h;fs=25kHz,fc=1.5kHz.word0FF78h

.word0Fe56h.word0FF78h.word59hend_coff.word8ch;h(0)

.text_c_int00:

ssbxintm;disableallinterrupt!

stm#VECTOR,ar0;vector'startaddress->ar0

st#INSTR_B,*ar0(#OFF_INTR_3)

st#fir,*ar0(#OFF_INTR_3+1) ;initA/Dintvector!

ld#temp,dp;thefollowingcodesforMAC

st#win_data,t_ar3;ar3->2000hdatawindows

st#filter_coff,t_ar2;ar2->filtercoff

xora,a;cleara

xorb,b;clearb

rsbxintm;enableallint!g:idle1

bg;;Thesecodesmaybecalledbyserialrevint!allregistersdon't;change!;Whenenterthissubroutine,theADdatahasbeenputinA!;fir:

ld#temp,dp

stla,temp;a->ADdata!

calllow_pass_mac

lda,-16,b

ld#0,dp

stlb,2,TDXR;sentresulttoDA!

rete;********************************************************;LOWPASSFILTER(useMAC);Inputdata->A,outputdata->A;usedAR2->coff,AR3->data_buffer!;********************************************************low_pass_mac:

pshmst1

pshmst0

pshmar0

pshmbk

mvdm#t_ar2,ar2;restorear2

mvdm#t_ar3,ar3;restorear3

stm#1,ar0

stm#N,bk;setcircularaddressingsize

stla,*ar3+%

rptza,#(N-1);0->a,thenrepeat34times

mac*ar2+0%,*ar3+0%,a;doneFIRfilter,resultina

mvmdar3,#t_ar3;savear3

mvmdar2,#t_ar2;savear2

popmbk

popmar0

popmst0

popmst1

ret程序、数据的存储器安排程序功能框图相关外设的准备：DSP，AD/DA，……相关外设的软件设置：McBSP串口初始化、AC01的初始化、……硬件电路设计、调试软件设计、调试FIR滤波器的工程实现初始化串口初始化AC01

等待新数据？调用滤波程序串口发送中断服务程序串口接受中断服务程序是llCCS图形工具显示FIR滤波效果图形显示FIR输入/输出频谱FIRS指令来实现FIR滤波器一种有限单位冲激响应呈现对中心点对称的FIR滤波器。长度为N的线性相位FIR滤波器的输出表达式为：FIRS指令来实现FIR滤波器FIRSXmem,Ymem,pmad含义：FIRS指令来实现FIR滤波器

B+（A(32-16))xPmad->B(Xmem+Ymem)<<16->APAR++Pmad寻址FIR滤波器系数，Xmem和Ymem分别指向窗口数据的上下两部分。16点FIRS滤波数据存放*ar2*ar3new*ar2++16点FIRS滤波数据存放FIR*AR2+0%,*AR3+0%#FIR_COFAR2--,AR3-=2*ar2*ar3new*ar2++使用FIRS指令完成滤波利用FIRS指令,需要将输入数据缓冲分成两个,大小为N/2。初始状态将AR2指到缓冲区1的顶部，AR3指到缓冲区2的底部。每次滤波之前,应先将缓冲区1顶部的数据移到缓冲区2的底部,新来的一个样本存储到缓冲区1中时,并对缓冲区1指针AR2加1(使用循环寻址)。使用FIRS指令完成滤波处理器然后使用FIRS指令进行乘加运算。当然,在使用FIRS指令前,需要预先计算一次求和,以初始化A。在RPTZ重复指令和循环寻址的配合下,完成FIR滤波。滤波完成后，需要对两个数据缓冲的指针进行修正，以便对下一个点进行处理。将Buffer1的指针减1和Buffer2的指针减2，使他们指向各自缓冲的数据队列的最后。STM #1,AR0

；AR0=1STM #（N/2）,BK

；BK=N/2,循环寻址BUFFER大小为NMVDD*ar2,*ar3

；更新Buffer2STL A,*ar2+%

；更新滤波窗口中的采样数据ADD*ar2+0%,*ar3+0%

；初始化ARPTZB,#(N/2-1)

；重复FIRS指令N/2次,先将B清零FIRS*ar2+0%,*ar3+0%,filter_coff+N/2

；完成滤波计算。注意FIR滤波系数存放；在程序存贮区，filter_coff为系数起始地址MAR*ar2-%

；修改Buffer1指针MAR*+ar3(-2)%

；修改Buffer2指针使用带FIRS指令的循环寻址模式实现FIR滤波器，程序片段如下：(输入数据在AL中,滤波结果在B中)FIRSADDAcx,ACy,Cmem,Xmem,YmemFIRSADD指令-C5500

ACy=ACy+(ACx*Cmem)::ACx=(Xmem<<16)+(Ymem<<16)FIRADD.CR可使得后面的Cmem，Xmem，Ymem间接寻址都是用循环寻址。在满足C5500的并行条件时，可以与其它指令并行执行。C6000实现滤波器-乘累加loop:

ldh.d1 *A8++,A2

||ldh.d2 *B9++,B3

nop 4 mpy.m1x A2,B3,A4

nop

add.l1 A4,A6,A6 sub.l2 B0,1,B0

[b0]b.s1 loop

nop 5Whatcanyou

putinparallel?LoadInstructionsMPY2x=+a1a0A0LDW.D1*A4++,A0x1x0B0||LDW.D2*B4++,B0A2A3a1*x1a0*x0MPY2A0,B0,A3:A2||MPYH.M2A0,B0,B5+a1x1+a3x3...a0x0+a2x2...ADD.L1A2,A6,A6||ADD.L2A3,B6,B6A6B6finalsumADD.L1A6,B6,A4A4+DOTP2withLDDW=+a2a0A1:A0LDDW.D1*A4++,A1:A0||LDDW.D2*B4++,B1:B0A2B2a3*x3+a2*x2a1*x1+a0*x0DOTP2A0,B0,A2||DOTP2A1,B1,B2+intermediatesumADDA2,A3,A3a1a3:A5x2x0B1:B0x1x3:finalsumADDA3,B3,A4A4+||ADDB2,B3,B3intermediatesumA3B3InCh8,we'llgetalltheseinstructionsworkinginparalleld0*c0d1*c1d2*c2d3*c3+yd4*c4d5*c5d6*c6d7*c7coefdataBlockRealFIRfor(i=0;I<ndata;i++){ sum=0; for(j=0;j<ncoef;j++){ sum=sum+(d[i+j]*c[j]); }

y[i]=sum;}loopIteration[i,j][0,0][0,1]d0c0d1c1d1c0d2c2d2c1d3c3d3c2...BlockRealFIRExample(DDTOPL2)for(i=0;I<ndata;i++){ sum=0; for(j=0;j<ncoef;j++){ sum=sum+(d[i+j]*c[j]); }

y[i]=sum;}loopIteration[i,j][0,0][0,1]d0c0

+d1c1d1c0+d2c1d2c2d3c3d3c2...Four16x16multipliesIneach.Muniteverycycle

addsupto8MACs/cycle,or 8000MMACSBottomLine:TwoloopiterationsforthepriceofoneDDOTPL2d3d2:d1d0,

c1c0,

sum1:sum0loopIteration[i,j][0,0][0,1][0,2][0,3][0,4][0,5]d0c0

+d1c1d1c0+

d2c1d2c2d2c0d3c3d3c2d3c1d3c0d4c4d4c3d4c2+d5c3d4c1d4c0d5c5d5c4d5c2+d6c3d5c1d6c0d6c6d6c5d6c4d6c2d6c1d7c7d7c6d7c5d7c4d7c3d6c2d8c7d8c6d8c5d8c4d6c3

DDOTPL2.M1d3d2:d1d0,c1c0,sum1:sum0||DDOTPL2.M2d7d6:d5d4,c3c2,sum3:sum2ParallelDDOTPL2’sRealBlockFIR[A_j]SPLOOPD4||MVC.S2B_i0,ILC;setILC||ADD.L2B_i0,1,B_i0;T/4||ADDAH.D2B_DLYaddr,nCoefs-1+4,B_DLYOUTaddr||MVK.S1nCoefs/4-3,A_T||ADDAB.D1A_DLYaddr,8,A_DLYaddr*-stageA*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_DLYaddr[1],A_dbda:A_d9d8;[3,1]||^LDDW.D2T2*B_INaddr++,B_TEMP1:B_TEMP0;ld1stinput||^MVC.S2B_i0,RILC;setRILCSPMASK||^LDDW.D2T1*B_INaddr++,A_TEMP1:A_TEMP0;ld2ndinput||^ZERO.L1A_st;clearstflag||^ADDAB.D1DP,outputs+8,A_OUTaddr||^MVK.S2nCoefs/4-1,B_TC||^MVK.S1nCoefs/4-1,A_TC*-stageB*SPMASK||^SUB.L2XA_OUTaddr,8,B_OUTaddrNOP1

DMV.S2XB_d5d4,A_d3d2,B_d5d4_:B_d3d2_;[7,1]||DDOTPL2.M1A_d3d2:A_d1d0,A_c1c0,A_pb:A_pf;[7,1]||DDOTPL2.M2B_d7d6:B_d5d4,B_c3c2,B_p2:B_p6;[7,1]SPMASK||DMV.S1XA_d9d8,B_d7d6,A_d9d8_:A_d7d6_;[8,1]||DDOTPL2.M2B_d5d4_:B_d3d2_,B_c3c2,B_pa:B_pe;[8,1]||DDOTPL2.M1A_dbda:A_d9d8,A_c3c2,A_p0:A_p4;[8,1]||^STNDW.D2T2B_TEMP1:B_TEMP0,*B_DLYOUTaddr++;st1stinputNewdouble-throughput16-bitmpyinstructionsABCABCABCDDDABCD4ABCABCABCDDDABCDO1O2O3O4O5e1e2,p1e3,p2e4,p3p4DDOTP4DDOTP4(.unit)src1,src2,dst_o:dst_eBlockRealFIRC64xImplementationDOTP2

d1d0,c1c0,s0;d[1]*c[1]+d[0]*c[0]

DOTP2

d2d1,c1c0,s1;d[2]*c[1]+d[1]*c[0]C64xPlusImplementationDDOTPL2d3d2:d1d0,c1c0,s1:s0Reduces.MrequirementbyhalfC64x: 194cycles,624bytes(N=40,T=16)

C64x+: 126cycles,496bytes(N=40,T=16)N=Lengthofexampleblock,T=datasizeCMPY...ComplexMultiply(CMPY)A0r1i1xxA1r2i2==CMPYA0,A1,A3:A2r1*r2-i1*i2:i1*r2+r1*i2

32-bits

32-bits

Four16x16multipliesper.MunitUsingtwoCMPYs,atotalofeight16x16multipliespercycleFloating-pointversion(CMPYSP)uses:64-bitinputs(registerpair)128-bitpackedproducts(registerquad)Youthenneedtoadd/subtracttheproductstogetthefinalresultsingle.MunitComplexMPYCMPY(.M)src1,src2,dst_o:dst_eFour16-bitinputs(real0,imag0,real1,imag1)Two32-bitoutputs(real,imag)CMPYR(.M)src1,src2,dstFour16-bitinputs(real0,imag0,real1,imag1)Onepacked32-bitoutput(16bitreal,16bitimag)Roundedbyadding2**15CMPYR1(.M)src1,src2,dstFour16-bitinputs(real0,imag0,real1,imag1)Onepacked32-bitoutput(real,imag)Roundedbyadding2**14plusleftshiftBlockComplexFIRC64xImplementationDOTP2 dre_dim,cim_cre,pim

DOTPN2 dre_dim,cre_cim,preC64xPlusImplementationCMPY dre_dim,cre_cim,pre:pimReduces.MrequirementbyhalfC64x: 674cycles,572bytes(N=40,T=16)

C64x+: 344cycles,460bytes(N=40,T=16)N=Lengthofexampleblock,T=datasizeDSPLIBOptimizedDSPFunctionLibraryforCprogrammersusingC62x/C67xandC64xdevicesTheseroutinesaretypicallyusedincomputationallyintensivereal-timeapplicationswhereoptimalexecutionspeediscritical.Byusingtheseroutines,youcanachieveexecutionspeedsconsiderablyfasterthanequivalentcodewritteninstandardANSIClanguage.Andtheseready-to-usefunctionscansignificantlyshortenyourdevelopmenttime.TheDSPlibraryfeatures:C-callableHand-codedassembly-optimizedTestedagainstCmodelandexistingrun-time-supportfunctionsAdaptive

filteringMathDSP_firlms2DSP_dotp_sqrCorrelationDSP_dotprodDSP_autocorDSP_maxvalFFTDSP_maxidxDSP_bitrev_cplxDSP_minvalDSP_radix2DSP_mul32DSP_r4fftDSP_neg32DSP_fftDSP_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_r8DSP_fltoq15DSP_fir_symDSP_minerrorDSP_iirDSP_q15toflTechnicalTrainingOrganizationT

TOFFT实现教材7章(7.3小节及7.4小节)FFT是数字信号处理中重要的工具FFT是一种高效实现离散付氏变换的算法。(8.2.2小节介绍了Goertzel算法）

离散付氏变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由付氏逆变换到时域。DFT的定义为：DFT的定义可以方便的把它改写为如下形式：WN（旋转因子）的周期性是DFT的关键性质之一。为了强调起见，常用表达式WN取代W以便明确其周期是N。FFT是DFT的快速算法由DFT的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N点DFT需要（N-1）2次复数乘法和N（N-1）次加法。FFT的基本思想在于，将原有的N点序列序列分成两个较短的序列，这些序列的DFT可以很简单的组合起来得到原序列的DFT。FFT是DFT的快速算法例如，若N为偶数，将原有的N点序列分成两个（N/2）点序列，那么计算N点DFT将只需要约[(N/2)2·2]=N2/2次复数乘法。即比直接计算少作一半乘法。该处理方法可以反复使用，即（N/2）点的DFT计算也可以化成两个（N/4）点的DFT（假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT运算的情况。FFT是DFT的快速算法比如，一个N=8点的FFT运算按照这种方法来计算FFT可以用下面的流程图来表示：实数FFT运算对于离散傅立叶变换（DFT）的数字计算，FFT是一种有效的方法。一般假定输入序列是复数。当实际输入是实数时，利用对称性质可以使计算DFT非常有效。一个优化的实数FFT算法是一个组合以后的算法。原始的2N个点的实输入序列组合成一个N点的复序列，之后对复序列进行N点的FFT运算，最后再由N点的复数输出拆散成2N点的复数序列。实数FFT运算使用这种方法，在组合输入和拆散输出的操作中，FFT运算量减半。这样利用实数FFT算法来计算实输入序列的DFT的速度几乎是一般复FFT算法的两倍。本实验就用这种方法实现了一个256点实数FFT（2N=256）运算。实数FFT运算序列的存储分配基二实数FFT运算的算法第一步，输入数据的组合和位倒序：（1）把输入序列作位倒序，是为了在整个运算最后的输出中得到的序列是自然顺序。（2）原始的输入的2N=256个点的实数序列复制放到标记有“d_input_addr”的相邻单元，当成N=128点的复数序列d[n]。其中，奇数地址是d[n]的实部，偶数地址是d[n]的虚部。复数序列经过位倒序，存储在数据处理缓冲器中“fft-data”。bit_rev:

STM#d_input_addr,ar3

；在AR3中放入输入地址

STM#fft_data,ar7

；在AR7中放入处理后输出的地址

MVMMDATA_PROC_BUF,ar2;AR2中装入第一个位倒序

;数据指针

STM#K_FFT_SIZE-1,BRC

STM#K_FFT_SIZE,ar0;AR0=输入数据数目的一半(128)

RPTBbit_rev_end

MVDD*ar3+,*ar7+ ;将原始输入缓冲中的数据放入到

;位倒序缓冲中去之后输入缓冲

;(AR3)指针加1,位倒序缓冲(AR2)指

;针也加1

MVDD*ar3-,*ar7+ ;将原始输入缓冲中的数据放入到

;位倒序缓冲中去之后输入缓冲(AR3);指针减1以保证位倒序寻址正确

MAR*ar3+0B

；按位倒序寻址方式修改AR3bit_rev_end:基二实数FFT运算的算法第二步，N点复数FFT：（1）在数据处理缓冲器里进行N点复数FFT运算。由于在FFT运算中要用到旋转因子WN，它是一个复数。我们把它分为正弦和余弦部分，用Q15格式将它们存储在两个分离的表中。每个表中有128项，对应从0度到180度。因为采用循环寻址来对表寻址，128=27<28，因此每张表排队的开始地址就必须是8个LSB位为0的地址。基二实数FFT运算的算法我们把128点的复数FFT分为七级来算，第一级是计算两点的FFT蝶形结，第二级是计算四点的FFT蝶形结，然后是八点、十六点、三十二点六十四点、一百二十八点的蝶形结计算。最后所得的结果表示为：其中，R[k]、I[k]分别是D[k]的实部和虚部。D[k]=F{d[n]}=R[k]+jI[k]这一步中，实现FFT计算的具体程序如下：;stage1：计算FFT的第一步，两点的FFT.asgAR1,GROUP_COUNTER ；定义FFT计算的组指针

.asgAR2,PX ；定义AR2为指向参加蝶形运算第一个数据的指针

.asgAR3,QX ；定义AR3为指向参加蝶形运算第二个数据的指针

.asgAR4,WR ；定义AR4为指向余弦表的指针

.asgAR5,WI ；定义AR5为指向正弦表的指针

.asgAR6,BUTTERFLY_COUNTER

；定义AR6为指向蝶形结的指针

.asgAR7,DATA_PROC_BUF

；定义在第一步中的数据处理缓冲指针

.asgAR7,STAGE_COUNTER

；定义剩下几步中的数据处理缓冲指针

pshmst0

pshmar0

pshmbk

；保存环境变量

SSBXSXM

；开启符号扩展模式

STM#K_ZERO_BK,BK

；让BK=0使*ARn+0%=*ARn+0

LD#-1,ASM

；为避免溢出在每一步输出时右移一位

MVMMDATA_PROC_BUF,PX

；PX指向参加蝶形结运算的

;第一个数的实部（PR）

LD*PX,16,A

；AH:=PR

STM#fft_data+K_DATA_IDX_1,QX

；QX指向参加蝶形

;结运算的第二个数的实部（QR）

STM#K_FFT_SIZE/2-1,BRC

；设置块循环计数器

RPTBDstage1end-1

；语句重复执行的范围到地址

;stage1end-1处

STM#K_DATA_IDX_1+1,AR0

；延迟执行的两字节的指令

;（该指令不重复执行）

SUB*QX,16,A,B

；BH:=PR-QR

ADD*QX,16,A

；AH:=PR+QR

STHA,ASM,*PX+

；PR':=(PR+QR)/2

STB,*QX+

；QR':=(PR-QR)/2

||LD*PX,A

；AH:=PI

SUB*QX,16,A,B

；BH:=PI-QI

ADD*QX,16,A

；AH:=PI+QI

STHA,ASM,*PX+0%

；PI':=(PI+QI)/2

STB,*QX+0%

；QI':=(PI-QI)/2

||LD*PX,A

；AH:=nextPRstage1end:(R[0]+jI[0])+(R[1]+jI[1])*W0

;Stage2：计算FFT的第二步，四点的FFT

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,16,A

；AH:=PR

RPTBDstage2end-1

；语句重复执行的范围到地址

;stage1end-1处

STM#K_DATA_IDX_2+1,AR0；初始化AR0以被循环寻址;以下是第二步运算的第一个蝶形结运算过程

SUB*QX,16,A,B

；BH:=PR-QR

ADD*QX,16,A

；AH:=PR+QR

STHA,ASM,*PX+

；PR':=(PR+QR)/2

STB,*QX+

；QR':=(PR-QR)/2

||LD*PX,A

；AH:=PI

SUB*QX,16,A,B

；BH:=PI-QI

ADD*QX,16,A

；AH:=PI+QI

STHA,ASM,*PX+

；PI':=(PI+QI)/2

STHB,ASM,*QX+

；QI':=(PI-QI)/2;以下是第二步运算的第二个蝶形结运算过程

MAR*QX+

；QX中的地址加一

ADD*PX,*QX,A

；AH:=PR+QI

SUB*PX,*QX-,B

；BH:=PR-QI

STHA,ASM,*PX+

；PR':=(PR+QI)/2

SUB*PX,*QX,A

；AH:=PI-QR

STB,*QX

；QR':=(PR-QI)/2

||LD*QX+,B

；BH:=QR

STA,*PX ;PI':=(PI-QR)/2

||ADD*PX+0%,A

；AH:=PI+QR

STA,*QX+0%

；QI':=(PI+QR)/2

||LD*PX,A

；AH:=PRstage2end:;Stage3thruStagelogN-1：从第三步到第六步的过程如下(略）基二实数FFT运算的算法第三步，分离复数FFT的输出为奇部分和偶部分：分离FFT输出为相关的四个序列：RP、RM、IP和IM，即偶实数，奇实数、偶虚数和奇虚数四部分，以便第四步形成最终结果。利用信号分析的理论我们把D[k]通过下面的公式分为偶实数RP[k]、奇实数RM[k]、偶虚数IP[k]和奇虚数IM[k]：RP[k]=RP[N-k]=0.5*(R[k]+R[N-k])RM[k]=-RM[N-k]=0.5*(R[k]-R[N-k])IP[k]=IP[N-k]=0.5*(I[k]+I[N-k])IM[k]=-IM[N-k]=0.5*(I[k]-I[N-k])RP[0]=R[0]IP[0]=I[0]RM[0]=IM[0]=RM[N/2]=IM[N/2]=0RP[N/2]=R[N/2]IP[N/2]=I[N/2]

这一过程的程序代码(略）基二实数FFT运算的算法第四步，产生最后的N=256点的复数FFT结果:产生2N=256个点的复数输出，它与原始的256个点的实输入序列的DFT一致。输出驻留在数据缓冲器中。通过下面的公式由RP[k]、RM[n]、IP[n]和IM[n]四个序列可以计算出a[n]的DFT：AR[k]=AR[2N-k]=RP[k]+cos(kπ/N)*IP[k]-sin(kπ/N)*RM[k]AI[k]=-AI[2N-k]=IM[k]-cos(kπ/N)*RM[k]-sin(kπ/N)*IP[k]AR[0]=RP[0]+IP[0]AI[0]=IM[0]-RM[0]AR[N]=R[0]-I[0]AI[N]=0其中：

A[k]=A[2N-k]=AR[k]+jAI[k]=F{a(n)}

这一过程的程序代码(略）基二实数FFT运算的算法计算所求信号的功率：由于最后所得的FFT数据是一个复数，为了能够方便的在虚拟频谱仪上观察该信号的特征，我们通常对所得的FFT数据进行处理——取其实部和虚部的平方和，即求得该信号的功率。CCS中的图形工具显示FFT结果DSPLIBOptimizedDSPFunctionLibraryforCprogrammersusingC62x/C67xandC64xdevicesTheseroutinesaretypicallyusedincomputationallyintensivereal-timeapplicationswhereoptimalexecutionspeediscritical.Byusingtheseroutines,youcanachieveexecutionspeedsconsiderablyfasterthanequivalentcodewritteninstandardANSIClanguage.Andtheseready-to-usefunctionscansignificantlyshortenyourdevelopmenttime.TheDSPlibraryfeatures:C-callableHand-codedassembly-optimizedTestedagainstCmodelandexistingrun-time-supportfunctionsAdaptive

filteringMathDSP_firlms2DSP_dotp_sqrCorrelationDSP_dotprodDSP_autocorDSP_maxvalFFTDSP_maxidxDSP_bitrev_cplxDSP_minvalDSP_radix2DSP_mul32DSP_r4fftDSP_neg32DSP_fftDSP_recip16DSP_fft16x16rDSP_