DSP常见算法的实现

3.6 常见的算法实现

在实际应用中虽然信号处理的方式多种多样，但其算法的基本要素却大多相同，在本节中介绍几种较为典型的算法实现，希望通过对这些例子（单精度，16bit ）的分析，能够让大家熟悉DSP 编程中的一些技巧，在以后的工作中可以借鉴，达到举一反三的效果。

1. 函数的产生

在高级语言的编程中，如果要使用诸如正弦、余弦、对数等数学函数，都可以直接调用运行库中的函数来实现，而在DSP 编程中操作就不会这样简单了。虽然TI 公司提供的实时运行库中有一些数学函数，但它们所耗费的时间大多太长，而且对于大多数定点程序使用双精度浮点数的返回结果有点“大材小用”的感觉，因此需要编程人员根据自身的要求“定制”数学函数。实现数学函数的方法主要有查表法、迭代法和级数逼近法等，它们各有特点，适合于不同的应用。

查表法是最直接的一种方法，程序员可以根据运算的需要预先计算好所有可能出现的函数值，将这些结果编排成数据表，在使用时只需要根据输入查出表中对应的函数值即可。它的特点是速度快，但需要占用大量的存储空间，且灵活度低。当然，可以对上述查表法作些变通，仅仅将一些关键的函数值放置在表中，对任意一个输入，可根据和它最接近的数据采用插值方法来求得。这样占用的存储空间有所节约，但数值的准确度有所下降。

迭代法是一种非常有用的方法，在自适应信号处理中发挥着重要的作用。作为函数产生的一种方法，它利用了自变量取值临近的函数值之间存在的关系，如时间序列分析中的AR 、MA 、ARMA 等模型，刻画出了信号内部的特征。因为它只需要存储信号模型的参量和相关的状态变量，所以所占用的存储空间相对较少，运算时间也较短。但它存在一个致命的弱点，由于新的数值的产生利用了之前的函数值，所以它容易产生误差累积，适合精度要求不高的场合。

级数逼近法是用级数的方法在某一自变量取值范围内去逼近数学函数，而将自变量取值在此范围外的函数值利用一些数学关系，用该范围内的数值来表示。这种方法最大的优点是灵活度高，且不存在误差累积，数值精度由程序员完全控制。该方法的关键在于选择一个合适的自变量取值区间和寻找相应的系数。

下面通过正弦函数的实现，具体对上述三种方法作比较。 查表法较简单，只需要自制一张数据表，也可以利用C5400 DSP ROM 内的正弦函数表。 迭代法的关键是寻找函数值间的递推关系。假设函数采样时间间隔为T ，正弦函数的角频率为ω，那么可以如下推导：

令()()()T T ω?β?αω?-+=+sin sin sin 等式的左边展开为

T T side left ω?ω?sin cos cos sin _+=

等式的右边展开为

()T T side right ω?βωα?sin cos cos sin _-+=

对比系数，可以得到1,cos 2-==βωαT 。令nT =?，便可以得到如下的递推式： [][][]21cos 2---=n s n s T n s ω

令[][]T A s s ωsin 2,01-=-=-，逐一迭代就能够获得采样间隔为T 的正弦序列了。从迭代公式可以更清楚地看出，这种方法存在误差累积。 再来看级数逼近法，首先需要寻找一个合适的自变量取值区间和寻找相应的系数。从正弦函数的对称性可知，只需要计算取值在],0[π内的函数就可以推断出所有取值范围内的函数。接下来寻找系数，对于()x sin 可以作如下变换()()sin_new(y)sin sin ==y x π，那么y 的取值范围在]1,0[内，而sin_new( )与sin( )同构，所以在下面的分析中将sin_new( )替代sin( )，提到的正弦函数即指sin_new( )。若汇编编程时采用的数据为Q15格式，那么取值与实际的弧度的对应关系如下图所示。

π

2

π

-

图3- 算法取值与弧度的对应关系

在]1,0[内，正弦函数的修正级数（五次）展开如下式： 5

4

3

2

1.800293x

0.5446778x

5.325196x

-x

0.02026367 3.140625x sin_new(x)+++=

根据上式，可以写出正弦函数的生成程序。 ;compute polynomial stlm A, T ;T=x stm #SinePolyCoeff, AR2 ld *AR2+, 16, A ;AH=c5 ld

*AR2+, 16, B ;BH=c4

poly *AR2+ ;A=c5*x+c4

poly *AR2+ ;A=c5*x^2+c4*x+c3 poly *AR2+ ;A=c5*x^3+c4*x^2+c3*x+c2 poly *AR2+ ;A=c5*x^4+c4*x^3+c3*x^2+c2*x+c1 mpya A sfta

A, 3

;adjust AH to Q15

SinePolyCoeff: ;in Q12 format .word 0x1cce ; 1.800293 (coef for x^5 = c5)

.word 0x08b7 ; 0.5446778

(coef for x^4 = c4)

.word 0xaacc ; -5.325196 (coef for x^3 = c3)

.word 0x0053 ; 0.02026367 (coef for x^2 = c2) .word 0x3240 ; 3.140625 (coef for x^1 = c1)

在编程过程中，使用到了POLY 语句，它能够使多项式的计算简洁快速地完成。该函数的结果可以通过实验X 来验证。

2. FIR 滤波器的实现

FIR 滤波器由于具有线性相位而且延迟能够确定，因而在信号处理中广泛应用。FIR 的基本模型如下图所示。

图3- FIR 模型

其数学表达式为[][][]∑-=-=

10

N i i n x i h n y ，根据该表达式可以给出一种最为直接的实现形

式。直接形式中采用线性地址来存放数据，如图3- 所示。

图3- 直接形式

程序中可以采用MACD 来实现程序如下： ld *(Input), A stl A, *(x_n)

stm #x_n_Nm1, AR2 mpy *AR2-, #h_Nm1, B rpt

#N-2

macd *AR2-, h_Nm2, B

程序首先将新的数据放置在x[n]处，然后对状态缓存由下而上计算，在循环语句中每执行一次状态变量自动向下移动一级，即自动更新。它的计算量基本为N个时钟周期。

当然，数据存放还可以采用循环缓存。另外，由于FIR的系数存在对称性，其结构如图3- 所示。那么可以利用这个特点，实现更为快速的计算。

图3- 对称结构的FIR

为计算方便将状态变量分别存放在两个缓存区间内，一块命名为Buffer_new，存放图3- 上半部分的数据，另一块命名为Buffer_old，存放图3- 下半部分的数据。它们都当循环缓存使用，大小为FIR阶数的一半，即N/2（常数SIZE）。根据上述原理编写的汇编程序如下：stm #Buffer_new, AR2

stm #Buffer_old, AR3

stm #SIZE, BK

stm # -1, AR0

fir_loop:

;read input

ld *(Input), A

stl A, *AR2

;filtering

add *AR2+0%, *AR3+0%, A

rptz B, #SIZE-1

firs *AR2+0%, *AR3+0%, FIR_Coeff

sth B, *(Output) ;store output

mar *+AR2(2)%

mar *AR3+%

mvdd *AR2, *AR3+0% ;update buffer

b fir_loop

为方便理解，在图3- 中，将状态变量更新的过程标明，左边的是Buffer_new，右边的是Buffer_old。开始时，指针AR2和AR3分别指向Buffer_new和Buffer_old中的x[0]与x[-19]，将最“新”的数据存进Buffer_new（步骤1）。利用FIRS实现FIR，结果放在BH中。计算结束后，AR2和AR3分别指向x[-1]和x[-18]（步骤2）。然后调整AR2，让它指向Buffer_new 中最“老”的数据x[-9] （步骤3）；调整AR3，让它指向Buffer_old中最“老”的数据x[-19]

（步骤4）。接下来进行数据更新，将Buffer_new 中最“老”的数据放进Buffer_old 中，成为最“新”的数据。最后AR2指向x[-9]的位置，这个位置将在下次计算时放置新的输入；AR3指向x[-18]，即Buffer_old 中最“老”的数据，便于下次计算（步骤5）。

图3- 对称结构的FIR 实现中状态变量的更新

利用对称结构的实现在计算量上有了减少，大约为N/2个时钟周期。当然，利用上述结构必须注意安排好数据的位置，以保证能进行正常的循环寻址。

3. IIR 滤波器的实现

IIR 滤波器也是数字信号处理的主要工具之一，但由于它不具备线性相位，而且无法准确知道其延迟，所以使用较FIR 少。下面，来观察一下IIR 的结构，其数学表达式如下：

∑∑==-+-=

M

i N i i i i n y a i n x b n y 0

1

][][][

从其数学公式可以看出，我们可以仿照在FIR 处理来直接实现。定义两段缓存分别对应于x[ ]和y[ ]，然后采用类似于FIR 的计算方式，采用MACD 语句，便能很快地完成IIR 滤波。

在直接实现中同时需要保存x[ ]和y[ ]，当其阶数较大时，会占用比较大的数据空间。为此，我们来考察IIR 的另一种结构。如图3- ，在这种结构中存储的状态变量为d[ ]，所以存储空间大大地减少了。

图3- IIR 的另一种结构

通过对FIR 和IIR 算法的分析，一方面向读者介绍这些基本处理中的设计技巧，另一方面也提醒读者在算法实现过程中应充分考察算法本身的特点，以求利用它们获得高效的运算和存储空间的节省。

4. FFT 的实现

在信号处理中，为突出信号的特征，经常在时域和频域之间作变换，而傅立叶变换是这当中最为典型的。基2的快速傅立叶变换有比特翻转排序和蝶形运算组成，前者在3.x 节已经作了介绍，这里着重介绍蝶形运算的实现。 蝶形运算是快速傅立叶变换中设计得极为精巧的部分，它充分揭示了傅立叶变换系数间内在的关系，而且还能实现同址计算，如图所示。完成一次蝶形运算需要一次复数的乘法和两次复数的加法。

N

图3- 蝶形运算的示意图

对于时间抽取的FFT 而言，在其第一级是乘法的系数为0

N W （也就是1），那么这个复数的乘法就名存实亡了，因而在计算FFT 时，第一级可以直接用复数的加法实现。第一级的程序如下： ;*** stage 1 *** stm #0, BK ld

#-1, ASM

stm #FFT_Data, PX stm #FFT_Data+K_DA TA_IDX_1, QX

ld *PX, 16, A ;AH=PX.x

stm #K_FFT_SIZE/2-1, BRC rptbd stage1end-1 stm #K_DA TA_IDX_1+1, AR0

sub *QX, 16, A, B ;BH=PX.x-QX.x add *QX, 16, A

;AH=PX.x+QX.x

sth A, ASM, *PX+ st B, *QX+ ||

ld *PX, A

;AH=PX.y sub *QX, 16, A, B ;BH=PX.y-QX.y add *QX, 16, A

;AH=PX.y+QX.y

sth A, ASM, *PX+0 st

B, *QX+0%

|| ld *PX, A stage1end:

在代码中，为方便与图3- 对应，使用.asg 伪指令将寄存器的名字替换成示意图中的表达方式，PX 和QX 分别指向蝶形运算的数据的地址。可见所有的操作完全是由加减实现的。

在第二级中乘法的系数为0N W 和4

/N

N

W （即+1和-j ），那么乘法操作只是影响参数的符

号，在复数的加减运算时只要弄清与-j 相乘造成的结果，所有的操作仍然可以只用加减法来实现。第二级的实现代码如下： ;*** Stage 2 *** stm #FFT_Data, PX stm #FFT_Data+K_DA TA_IDX_2, QX

ld *PX, 16, A ;AH=PX.x

stm #K_FFT_SIZE/4-1, BRC rptbd stage2end-1

stm #K_DA TA_IDX_2+1, AR0

; 1st butterfly

sub *QX, 16, A, B ;BH=PX.x-QX.x add *QX, 16, A

;AH=PX.x+QX.x

sth A, ASM, *PX+ st

B, *QX+ || ld *PX, A ;AH=PX.y

sub *QX, 16, A, B ;BH=PX.y-QX.y add *QX, 16, A

;AH=PX.y+QX.y

sth

A, ASM, *PX+

sth B, ASM, *QX+ ; 2nd butterfly

mar *QX+

add *PX, *QX, A ;AH=PX.x+QX.y

sub *PX, *QX-, B ;BH=PX.x-QX.y

sth A, ASM, *PX+

sub *PX, *QX, A ;AH=PX.y-QX.x

st B, *QX

|| ld *QX+, B ;BH=QX.x

st A, *PX

|| add *PX+0%, A ;AH=PX.y+QX.x

st A, *QX+0%

|| ld *PX, A

stage2end:

第二级与第一级相同，计算中不包含乘法。从第三级开始，乘法系数的特征就没有第一、第二级那样好了，所以只能直接采用图3- 的方法来计算，代码如下。

;*** Stage 3 through Stage logN ***

stm #K_TWID_TBL_SIZE, BK

st #K_TWID_IDX_3, *(d_twid_idx)

stm #K_TWID_IDX_3, AR0

stm #Cos, WR

stm #Sin, WI

stm #K_LOGN-2-1, STAGE_COUNTER

st #K_FFT_SIZE/8-1, *(d_grps_cnt)

stm #K_FLY_COUNT_3-1, BUTTERFLY_COUNTER

st #K_DA TA_IDX_3, *(d_data_idx)

stage:

stm #FFT_Data, PX

ld *(d_data_idx), A

add *(PX), A

stlm A, QX

mvdk *(d_grps_cnt), GROUP_COUNTER

group:

mvmd BUTTERFLY_COUNTER, BRC

rptbd butterflyend-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

st B, *PX ; PR':=((QR*WR+QI*WI)+PR)/2

|| sub *PX+, B ; B=PR-(QR*WR+QI*WI)

; || PX->PI

st B, *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

st B, *QX+ ; QI':=((QR*WI-QI*WR)+PI)/2 || sub *PX, B ; B=PI-(QR*WI-QI*WR)

ld *WR, T ; T := WR

st B, *PX+ ; PI':= (PI-(QR*WI-QI*WR))/2 || mpy *QX+, A; || PX->PR

; A:= QR*WR || QX->QI butterflyend:

; Update pointers for next group

pshm AR0

mvdk *(d_data_idx), AR0

mar *PX+0

mar *QX+0

banzd group, *GROUP_COUNTER-

popm AR0

mar *QX-

; Update counters and indices for next stage

ld *(d_data_idx), A

sub #1, A, B

stlm B, BUTTERFLY_COUNTER

stl A, 1, *(d_data_idx)

ld *(d_grps_cnt), A

stl A, ASM, *(d_grps_cnt)

ld *(d_twid_idx), A

stl A, ASM, *(d_twid_idx)

mvdk *(d_twid_idx), AR0

banz stage, *STAGE_COUNTER-

fft_end:

在FFT的程序编制中可以发现，在很多高级语言编程中忽略的地方，对于汇编编程时却需要仔细分析，对于能够减少计算量的部分应分割出来，并单独编写处理程序，以获得更好的计算性能。

通过对上面几种常见算法的分析，可以明白在算法实现应注意一下几个要点：一，数据的精度，对于不同的实现方法能够达到的数据精度是不相同的，必须根据实际要求，综合地考虑结果的数据精度、计算速度和占用存储空间等因素后，才能确定最终的实现方法。二，实现的手段，基于同一个目的可以采用不同的手段实现，它们在存储空间分配和计算量上都存在着差别，应根据算法在计算速度上的要求，选择合适的手段。当使用一些经过精心设计的高效结构时，汇编程序的可读性大大下降了，应注意作好技术文档，详细地解释实现的具体细节，便于维护和改进。三，分别处理，对于一些在高级语言编程中不存在差别的处理过程，在汇编编程时应仔细分析，当发现其中存在某些能够更充分运用DSP硬件结构的环节，可以考虑打破原有的软件结构，把这些过程独立出来，采用高效的方法来处理。当然，在实际编程中，使用到的技巧远远不只这些，还需要编程人员在实践中不断的摸索、尝试和总结。