IIR数字滤波器在定点DSP上的实现

2006年第 2 期 声学与电子工程 总第 82 期IIR 数字滤波器在定点DSP 上的实现马 洁, 董金明(北京航空航天大学电子信息工程学院, 北京，100083)摘要: 介绍了IIR(Infinite Impulse Response)数字滤波器在Motorola56300系列的定点DSP(Digital Signal Processor)上的实现方法，给出了双二阶滤波器的单精度和双精度两种设计，并比较了它们的精度。

同时还给出了源程序，对在Motorola 的DSP 上实现双二阶滤波器具有一定的参考价值。

关键词: 数字滤波器；定点DSP ；IIR ；双二阶滤波器随着数字化时代的到来，数字信号处理已经发展成为一门极其重要的学科和技术领域。

数字信号处理技术在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。

在数字信号处理应用中，数字滤波器是十分重要的部分。

本文主要介绍了数字滤波器在Motorola 的56300系列的定点DSP 上的实现，这里采用了单精度和双精度两种方法，并比较了这两种方法在实现精度上的差别。

并且介绍了数字滤波器设计的模块化以及可重复利用的工程化设计思路，使数字滤波器具备可重复利用的价值。

1 IIR 数字滤波器原理和DSP56300介绍1.1 IIR 数字滤波器的实现一个IIR 数字滤波器可用系统函数表示为：)()(1)(1z X z Y z a z bz H N k kk Mk kk=−=∑∑=−=− 由此式可得出表示输入输出关系的常系数线性差分方程为：∑∑==−+−=Mk k N k k k n x b k n y a n y 01)()()(数字滤波器可分为IIR 和FIR 两种类型。

不同类型的滤波器在实现的时候，可以采用多种形式的网络结构来实现。

IIR 滤波器的结构包括直接型，二阶节级联和并联结构。

FIR 包括有横截型，级联型结构，频率抽样型，快速卷积型结构。

一个好的滤波器结构应该具有易于控制滤波器性能，适合于模块化实现的特点。

对于客户来说，一般要求在四阶以内，并且要求阶数可控制。

根据这些要求，只要设计一个具有良好特性的双二阶滤波器，然后通过控制级联的个数，就可以实现任意阶数的滤波器。

为下文叙述方便，将单精度双二阶节滤波器记为Biquad ， 双精度双二阶节滤波器记为Biquad2。

1.2 DSP56300系列芯片介绍DSP56300的ALU （数据逻辑运算单元）包括完全流水线式的并行乘法器-累加器（MAC, Multiplier-Accumulator ）；四个24位独立的通用输入寄存器：X1、 X0、 Y1 和 Y0， 或者由这四个寄存器连接成两个48位输入寄存器：X 和 Y ；六个数据寄存器：A2、 A1、 A0 和B2, B1、 B0，可连接成为两个通用56位累加器A 和B ；在进行乘法、加减法和乘累加运算的同时，如果满足某些特定的条件，DSP 最多允许两条并行的数据读取，以提高运算速度。

另外，DSP56300支持线性寻址、循环寻址和逆向寻址三种寻址模式。

因为数字滤波器的实现结构实际上都是由乘加运算组成的。

因此在运算时可以充分运用乘累加指令，这样在一条语句里同时完成乘法和加法，另外可以并行地读取下次乘累加的操作数。

通常把滤波器的系数和输入输出的过去值分别放在X 存储空间和Y 存储空间，这样可以大大提供滤波器的工作效率。

总之，用DSP 来实现数字滤波器具有很多得天独厚的优势，可以实时地实现数字滤波器。

2 Biquad 和 Biquad2 在DSP 上的实现为了实现数字滤波器的模块化，便于在以后的项目中重复利用。

这里将Biquad 和Biquad2都写成了函数的形式，具有明确的与外界的接口，只要滤波器的系数和滤波中的过去值按照规定的顺序和位置放置，就可以十分方便的实现滤波。

注意，在这里系数要量化为24bits ，用Q22定标，待滤波的采样值也量化为24bits ，这里用Q23定标。

由于系数采用了Q22格式，所以在滤波之前，需要设置SR (状态寄存器) 的 [S1 S0] 为 [1 0] ，增大MAC 中舍入位置的标度，也使得从累加器读出到X 数据总线和Y 数据总线的数据左移一位。

双二阶滤波器要实现的差分方程为：[]1[1]2[2]0[]1[1]2[2]y n A y n A y n B x n B x n B x n +×−+×−=×+×−+×−2.1 Biquad 的实现Biquad 的结构图如图1所示，系数表按照 -A2、 –A1、 B2、 B1、B0的顺序在X 存储空间放置， 状态表按照x[n-1]、 y[n-1]、 y[n-2]、 x[n-2]的顺序在Y 存储空间放置。

图1 Biquad 的 结构图其汇编程序如下。

输入参数： r1指向状态表头，r4指向系数表头，n0为待处理的数据个数，r0指向待处理的数据块。

输出参数：用滤波后的数据更新数据块中的原始数据 ；加载初始化条件move #3,n1move y:(r1)+,y1 ；取出x[n -1]存入y1 move x:(r4)+,x0 y:(r1)+,a ；取出-A2存入x0，y[n -1]存入aasr a #-4,n4 move a,x1 y:(r1),y0 ；将y[n -1]存入x1, y[n -2]存入y0；对数据块中的数据进行滤波的主循环体部分dor n0,_loop ；重复从下一指令到loop 间的所有指令n0次mpy x0,y0,a x:(r4)+,x0 a,y:(r1)+；乘-A2和y[n -2]，并取出-A1,更新y[n -2] mac x0,x1,a x:(r4)+,x0 y:(r1),y0 ；乘累加并取出B2和x[n -2] mac x0,y0,a x:(r4)+,x0 y1,y:(r1)- ；乘累加并取出B1，并更新x[n-2] mac x0,y1,a x:(r4)+n4,x0 y:(r0),y1 ；乘累加并取出B0和x[n]macr x0,y1,a x:(r4)+,x0 y:(r1),y0 ；乘累加并经过舍入处理后得到输出采样值y[n],； 并取出-A2和y[n-2] move a,x1 a,y:(r0)+；输出计算出的输出采样值 y[n] _loop；为下一个数据块的滤波存储必要的状态值 move (r1)- move (r1)-move x:(r4)-,x0 y1,y:(r1)+ ；更新x[n-1] move x1,y:(r1)+n1 ；更新y[n-1] 2.2 Biquad2的实现因为用定点DSP 实现时，所有的数据都是定长的，运算也都是定点运算，因而会产生有限字长效应。

所产生的误差有系数量化引起的误差以及运算过程中的舍入误差。

产生误差在用定点 DSP 时是不能避免的，但可以通过一些方法减小误差。

对于IIR 滤波器而言，可以将反馈的滞后输出用双字表示。

这里给出了实现这种双精度IIR 数字滤波器的方法。

Biquad2和Biquad 的基本算法是一样的。

主要的区别在于滞后输出采用了双字表示。

这里还需要说明数的双字表示概念。

在DSP 中，数值的表示都是基于2的补码表示形式。

设N 是48 bits 的二进制数，N 表示的整数为N_D ，则N_D=N<2^47?N:N-2^48，N 表示的小数为N_F ，则N_F=N_D*2^-47。

设 Nls 是N 的最低有效24位，Nls_D 是Nls 表示的整数。

Nrnd 是N 按照2的补码舍入后的48位结果 (2的补码舍入法是指一切大于或等于一半的值都向上舍入，小于一半的所有值则舍去)。

则有： N_F=Nrnd_D*2^-47+(Nls_D*2^-23)*2^-24 更进一步来说，Nmr 是Nrnd 的最高有效24位，Nmr_D 是Nmr 表示的整数。

所以N 可以用两个24位的字来存储（Nmr,Nls ）。

因此N_F=Nmr_D*2^-23+(Nls_D*2^-23)*2^-24 设C 是24位字，C_D 和C_F 分别是它表示的十进制整数和十进制小数。

C_F*N_F=(2*C_D*Nmr_D+(ROUND(2*C_D*Nls_D))>>24)*2^-47。

这种算法对2的补码舍入法是正确的，但对于其他的舍入方法有可能导致错误。

系数的放置和Biquad 一致。

状态表按照yls[n -1], ymr[n -1], yls[n -2], ymr[n -2], x[n -2],x[n -1]的顺序，在Y 存储空间放置。

Biquad2的结构图如图2所示。

图2 Biquad2 的结构其汇编程序如下(这里给出了主循环体部分)，这个函数的入口参数和出口参数和Biquad 一致。

在主循环体部分之前，需要说明初始化条件：x1= ymr[n-1]，x0=-A2，y0=yls[n-2]，y1=yls[n-1]。

dor n0,_loop；重复从下一指令到loop 间的所有指令n0次 mpy x0,y0,a y:(r1)-,y0 ；乘-A2和yls[n-2]，并取出ymr[n-2]存入y0 mpy x0,y0,b x:(r4)+,x0 y1,y:(r1)+ ；乘-A2和ymr[n-2]，并取出-A1,更新yls[n-2]macr x0,y1,a x1,y:(r1)+；乘累加得到：-A2*yls[n-2]+(-A1)*yls[n-1]，并舍入处理，更新ymr[n-2]asr #24,a,a ；舍入后的结果右移24位 mac x0,x1,b y:(r1)+,y1 ；乘累加得到：-A2*ymr[n-2]+(-A1)* ymr[n-1]，并取出x[n-2]存入y1add b,a x:(r4)+,x0 y:(r1)-,y0 ；相加得到V[n]，并取出B2，x[n-1]存入y0mac x0,y1,a x:(r4)+,x0 y0,y:(r1)+ ；乘累加并取出B1，更新x[n-2] mac x0,y0,a x:(r4)+n4,x0 y:(r0),y0 ；乘累加并取出B0，x[n]存入y0 mac x0,y0,a y0,y:(r1)-n1 ；乘累加得到w[n],更新x[n-1]asl #1,a,b ；得到最终输出结果b=2*w[n] rnd a x:(r4)+,x0 y:(r1)+,y0；舍入处理后得到wmr[n]，并取出-A2和yls[n-2]move b0,y1 ；更新yls[n-1]存入y1 move a,x1 a,y:(r0)+；更新ymr[n-1]=2*wmr[n]存入x1，输出计算出的输出采样值y[n]_loop3 Biquad 和 Biquad2 比较3.1 比较存储空间和指令数表1 存储空间和指令数比较存储空间 (word) 指令数 （Mips Million instruction per second ）Biquad 20 0.393 Biquad230 0.635注: Mips 的计算是基于采样率为48000Hz ，每128个采样值调用一次Biquad 或者Biquad23.2 比较滤波精度 3.2.1测试方法(1) 用MATLAB 设计出采样率SR 为48000Hz, 截止频率Fc=80Hz 的巴特沃斯高通滤波器。