二次IIR均衡器(equalizer)

RY DSP开发 2005/07/11
也许很少有人知道,无论是古老的盒式录音机还是现代的流行数码音响设备,以及众多的音乐播放软件,其中绝大多数的均衡器都只是一系列简单的二次IIR滤波器组合而成。本文将详细介绍二次IIR均衡器的系数设计方法,以及如何用C语言和TI 5000的汇编高效地实现。

关于IIR滤波器的详细知识,应该可以在很多关于信号处理的书籍中找到,本文将不做详细介绍。

二次IIR滤波器的传递函数如下,

             b0 + b1*z^-1 + b2*z^-2
H(z) = -------------------------------------
a0 + a1*z^-1 + a2*z^-2

它看上去有6个参数,其实独立的只有5个,把上式中每个系数都除以a0,可以得到

                (b0/a0) + (b1/a0)*z^-1 + (b2/a0)*z^-2
H(z) = --------------------------------------------------------
1 + (a1/a0)*z^-1 + (a2/a0)*z^-2

传递函数可以直接转换为如下的直接1型计算公式:

y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]
- (a1/a0)*y[n-1] - (a2/a0)*y[n-2]

这里的n表示当前时刻(取样点),x是IIR滤波器的输入,y是输出,这个计算公式表明:n时刻的输出是n,n-1,n-2时刻的输入和n-1,n-2时刻的输出的线性组合。

实现这个公式的C语言程序不会比公式本身长多少,所以在写具体的计算程序之前,我们还是来看看如何设计二次IIR均衡器的系数b0,b1,b2,a0,a1,a2吧。

二次IIR均衡器有中心频率f0,增益dBgain和Q(Q值决定频率响应的宽窄),如下图所示。

[图缺]

数字滤波器还有一个额外的参数:取样频率Fs。
给出上面的4个参数,就可以按照下面的公式求出b0,b1,b2,a0,a1,a2了。

A  = sqrt( 10^(dBgain/20) )
w0 = 2*pi*f0/Fs
alpha = sin(w0)/(2*Q)
b0 = 1 + alpha*A
b1 = -2*cos(w0)
b2 = 1 - alpha*A
a0 = 1 + alpha/A
a1 = -2*cos(w0)
a2 = 1 - alpha/A

这个公式可以根据频率响应的特征列出方程组求得,由于二次IIR的频率响应的通用公式很容易通过其传递函数求得,因此只需要确定频率响应几个特征值,就可以确定传递函数的系数了。使用这个方法可以求出各种类型的二次IIR滤波器的系数公式,详细可以参考http://shoko.calarts.edu/pipermail/music-dsp/2001-March/008604.html

下 面我们来看看如何用C语言实现二次IIR滤波器。上文给出直接1型计算公式就是计算二次IIR滤波器全部算法了。不过在DSP具体实现中,为了能够实时地 处理声音信号,经常使用乒乓缓存和块处理等技术。无论是输入还是输出流,都有一个乒乓缓存与之对应。当系统正在写入或者读出乒缓存的时候,我们的程序就处 理乓缓存中的数据,当系统完成写入或者读出之后,就立刻转到乓缓存继续,而我们则开始对乒缓存处理。通过这样的交替,使得整个输入输出数据流可以不间断地 工作。所以为了实时处理数据的,必须一次能够处理一个乒缓存或者乓缓存,简单地说就是一次能处理一块数据。下面先给出二次IIR滤波器块处理程序,然后再 做分析。

void iir2(short *buf, float *pars, int size, int inc, float *prev_buf){
int i,j;
float out;
float *t;
short *p=buf;
t=prev_buf;
for(i=0;i
t[0]=(float)(*p);
out=0;
for(j=0;j<5;j++)
out+=t[j]*pars[j];
if(out>0x7fff) *p=0x7fff;
else
if(out<0x8000) *p=0x8000;
else
*p=(short)out;
t[2]=t[1];
t[1]=t[0];
t[4]=t[3];
t[3]=out;
p+=inc;
}
}

由 于要计算当前块的第一个输出的时候,必须使用前一个块的信息,所以我们使数组prev_buf来保存前一个块的信息。它的元素分别是:x0,x1,x2, y1,y2。x0表示当前的输入值,x1表示上一个输入值,x2则是上上个输入值,y1是上个输出值,y2是上上个输出值。
Pars是IIR滤波器的系数数组,由程序

for(j=0;j<5;j++)
out+=t[j]*pars[j];

可知,系数数组的值是[(b0/a0), (b1/a0), (b2/a0), -(a1/a0), -(a2/a0)]。

Buf 是滤波器的输入数据,它的长度为size(样本数)。一般的声音都是16bit的整型数据,所以它的类型是short *,而系数和中间的计算结果都用float保存。在实际的DSP程序中,为了获得最佳的计算速度,系数和中间结果也都使用short保存,这样的话就需要 做定点小数运算了。

最后,由于声音数据有stereo(双声道)和mono(单声道)之分,所以增加一个inc参数,这样如果是stereo输入的话,就把inc设置为2,只对一个声道的数据进行处理。

优化

尽 管程序看上去很简单,但是还是有许多优化工作需要做。通常的多段均衡器都是由一系列的二次均衡器串联而成,因此提高程序的效率可以大幅度地提高整个均衡器 的效率。我们上面给出的程序用输出数据覆盖输入数据,因此在块计算中不需要把最近的两个输出数据放入prev_buf,而可以直接使用buf中的已计算的 输出数据。当块计算结束的时候再把最后两个输出数据写入prev_buf。这样一来就不能用循环来计算输出数据了。需要把

out=0;
for(j=0;j<5;j++)
out+=t[j]*pars[j];
改写为
out=pars[0]*t[0];
out+= pars[1]*t[1];
out+= pars[2]*t[2];
out+= pars[3]* *(p-inc);
out+= pars[4]* (p-2*inc);

同 时计算块最开始的两个输出的时候,由于需要用上一个块的保留信息prev_buf,因此需要做特殊的处理,在块计算结束的时候,还需要把最后两个输出保存 入prev_buf。这样修改之后,对于块中的每个输出不需要写入prev_buf,一个输出可以节省两次赋值操作。同理,也可以节省输入数据的两次赋 值,不过这需要把运算结果储存到另外一个数组中,而不是直接覆盖原数据。

到了这一步,程序逻辑上的优化已经基本完成了。剩下的优化工作主要针对DSP芯片进行。

DSP优化
大多数DSP都不直接支持浮点运算,所以把浮点运算改写为定点小数运算,将能大幅度提高效率。而在DSP中通常有能够加快定点小数运算的指令,手写DSP汇编程序可以对程序进行终极优化。

下面针对TI5400系列的芯片介绍一下可以提高运算速度的指令。

累加指令MAC。使用MAC可以在一个指令周期完成加法和乘法,所以DSP很擅长做类似out+= pars[1]*t[1];的计算。

设置状态寄存器的FRCT位。FRCT位是专门为定点小数运算设定的,当FRCT=1的时候,乘法计算单元的输出会左移一位,这样如果我们取计算结果的高16位的话,就相当于把结果右移15位,这正好是Q15定点小数乘法所需要的移位操作。

设 置状态寄存器的OVM位。在一般的计算机芯片中,如果运算超出范围的话,就会从最大值转回到最小值,因此对于一个16位的短整型数来说,0x7fff+ 0x0001=0x8000。但是在DSP中,如果OVM=1的话,那么结果若超过最大(小)值,就保留结果为最大(小)值。这样,判断结果是否溢出的 if语句就不再需要了。

具体的汇编程序在TI5400的算法库中可以找到。注意算法库中的程序使用Q15定点小数来进行运算,这样要求所有 的系数都不大于1,对于系数大于1的情况,则采取取系数的1/2进行乘法运算(系数一般都是小于2的),然后把结果加两次到累加寄存器的办法。因为 TI5400的累加寄存器是40位的,所以不会在累加寄存器中产生溢出。算法库的程序只对a1>1进行了处理,而由前面的均衡器系数公式可知, b0,b1都有可能大于1,因此需要手工添加这部分程序。

另外,16位定点小数的运算精度有限,如果均衡器的中心频率太低,则系数都十分接 近1,2,-1,-2,这时的运算误差将非常大,以至于输出波形变形。根据个人经验,在44100Hz的取样频率时,为了获得较高的计算精度,对于中心频 率小于250Hz的均衡器最好采用32位定点小数进行运算。可惜TI5400不直接支持32位乘法,因此必须把乘法分开做,详细参见定点小数运算一文。虽 然TI5400的算法库也提供了32位二次IIR滤波器的汇编源程序,可是我并没有调试通过,因此就自己动手写了一个,对于单声道44100Hz的输入数 据的处理量大约是5MIPS左右,这几乎是16位程序的5倍。在DSP运用中经常会出现这种效率和效果的矛盾。

性能实测

为了在效率和效果之间做出正确的选择,需要对均衡器的性能进行实测。当实测值与理论值的误差较大的时候就应该考虑使用32位运算。

通 常我们把标准测量信号输入到系统中,录制系统的输出信号并分析输出信号和输入信号的频谱,就可以测得系统的频率响应。输入信号一般采用白噪声或者 Sweep音(频率滑动的正弦波),而用Sweep音测量的结果较为准确。由于频率响应的频率轴通常都使用对数坐标,所以Sweep音的频率变化应该按照 指数增加,这样只要看输出信号的振幅就知道系统的频率响应了。