1．概述

随着数字化技术的快速、深入发展,人们对数字化电子产品所产生的图像、图形以及声音等质量的要求越来越高。在实时数字处理过程中，与D/A和A/D转换相关的模拟信号重构过程是决定数字系统输出质量的关键^[1]。在声音的拾取过程及通过音响设备的传送过程中，由于设备或器件的原因，其幅度对频率的响应往往不一致，这样就达不到原来的听觉效果，而均衡器就是一个改变放大器频率响应的设备。一个好的均衡器能起到以下几个方面的作用：

(1) 能校正音频设备所产生的频率失真；

(2) 校正由于建筑学共振性的不均匀所带来的传输增益的频率失真；

(3) 调整某些易反馈的频率成分，抑制声反馈，提高会场增益；

(4) 在艺术创作中，利用它来刻画演员和乐器的音色修改，提高艺术效果。

现有的许多播放器都具有均衡器的功能，如Winamp(如图 1所示)等。

图 1 Winamp的均衡器

MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的滤波器设计工具箱为实现声音信号的数字滤波提供了十分方便的函数和命令，但MATLAB的计算速度慢。Visual C++是Windows平台下主要的应用程序开发环境之一，它能方便实现软件开发，开发的系统具有执行速度快等优点，故MATLAB与VC的混合编程正好结合了MATLAB强大的工具箱与VC很快的执行速度。本文将给出基于MATLAB与VC混合编程的数字均衡器的设计方法，功能包括：

(1)．WAV声音信号获取，即从文件(*.wav，*.au等)读取音频信息；

(2)．八段数字均衡器设定，具体按照表 1中的频段，其中前两个频段合并在一起；

(3)．滤波，包括生成低通、高通、带通、带阻的巴特沃斯、切比雪夫滤波器，并显示滤波器的频率特征曲线；

(4)．保存，保存滤波后的音频信息到文件；

(5)．播放，播放滤后音频信号并实时显示波形及频谱特性。

表 1频率的音感特征

30~60Hz	沉闷	如没有相当大的响度，人耳很难感觉。
60~100Hz	沉重	80Hz附近能产生极强的“重感”效果，响度很高也不会给人舒服的感觉，可给人以强烈的刺激作用。
100~200Hz	丰满
200~500Hz	力度	易引起嗡嗡声的烦闷心理。
500~1KHz	明朗	800Hz附近如提升10dB，会明显产生一种嘈杂感，狭窄感。
1K~2KHz	透亮	2800Kz附近明亮感关系最大。
2K~4Kz	尖锐	6800Hz形成尖啸，锐利的感觉。
4K~8Kz	清脆	3400Hz易引起听觉疲劳。
8K~16Kz	纤细	＞7.5KHz音感清彻纤细。
30~60Hz	沉闷	如没有相当大的响度，人耳很难感觉。

2．设计原理

2.1数字滤波器

2.1.1 数字滤波器的原理简介

数字滤波器的功能是把输入序列通过一定的运算，变换成输出序列。数字滤波器一般可用两种方法实现：一种是根据描述数字滤波器的数学模型或信号流程图，用数字硬件构成专用的数字信号处理机，即硬件方式；另一种是编写滤波器运算程序，在计算机上运行，即软件方式。考虑到软件的灵活性及易于实现，本文采用软件方式实现数字滤波器。

数字滤波器有无限冲激响应(IIR)和有限冲激响应(FIR)两种。下面分别介绍：

(1) IIR滤波器

IIR滤波器的特点是：单位冲激响应h(n)是无限长的；系统函数H(z)在有限长Z平面（0<|Z|<∞)有极点存在；结构上存在输出到输入的反馈，也即结构上是递归型的；因果稳定的IIR滤波器其全部极点一定在单位圆内。其系统函数为

                                                    (1)

计算机上实现时则需要用到差分方程的形式，如下

                                        (2)

IIR滤波器有四种基本的网络结构(具体参看文献[3])，直接I型、直接II型、级联型与并联型。其中直接I型需要2N个延迟单元，而直接II型只需要N个延迟单元。因此，用软件实现时，直接II型少占用存储单元。级联型则是将N阶IIR系统函数分解成二阶因式连乘积，并联型则是将系统函数化成部分分式之和，则可得到IIR数字滤波器的并联结构。

(2) FIR滤波器

IIR滤波器的特点是：系统的单位冲激响应h(n)是个有限长序列；系统函数|H(z)|在|z|>0处收敛，极点全部在z=0处(即FIR一定为稳定系统)；结构上主要是非递归结构，没有输出到输入反馈。但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。其系统函数的一般形式为

                                                                  (3)

对应的差分方程为

                                                            (4)

FIR系统的基本结构有直接型、级联型、快速卷积型、频率取样型等。

2.1.2 FIR与IIR滤波器的比较与选择

IIR滤波器可以用比FIR滤波器少的阶数来满足相同的技术指标，这样，IIR滤波器所用的存储单元和所用的运算次数都比FIR滤波器少。FIR滤波器可得到严格的相位，而IIR滤波器不能得到。事实上，IIR滤波器的选频特性越好，它的相位的非线性就越严重。因此在需要严格线性相位的情况下应该选择FIR滤波器。IIR滤波器可利用模拟滤波器现成的设计公式、数据和表格，因而计算工作量较小，对计算工具要求不高。FIR滤波器没有现在的设计公式，对计算工具要求较高，需要借助计算机来设计。另外，IIR滤波器主要是设计规格化的、频率特性为分段常数的标准低通、高通、带通、带阻和全通滤波器，而FIR滤波器可设计出理想正交变换器、理想微分器、线性调频器等各种网络，适应性较广。

总之，IIR和FIR这两种滤波器各有特点，在实际应用中空间选择中哪种滤波器，就从多方面的因素来考虑。例如用于语音通信的滤波器，对相位要求不是主要的，因此选用IIR滤波器较为合适，可以充分发挥其经济和高效的特点。而图像信号处理和数据传输等以波形携带信息的系统，对相位的线性要求较高，因此采用FIR滤波器较好。

对于数字均衡器，一方面是用于处理语言信号，另一方面需要用到频率特性分段的带通滤波器，因此应该IIR滤波器。下面介绍IIR滤波器的设计方法。

2.1.3 IIR数字滤波器的设计方法

实际中的数字滤波器都是用有限精度算法实现的线性非移变离散系统，设计IIR数字滤波器的方法主要有两种：一种是利用模拟滤波器的理论来设计，另一种是计算机辅助设计，即使用最优化技术设计。利用模拟滤波器的设计理论来设计IIR滤波器，就是首先根据实际要求设计一个模拟滤波器，然后再将这个模拟滤波器转换成数字滤波器。由于模拟网络综合理论已经发展得很成熟，故许多常用的模拟滤波器不仅有了简单而严格的设计分式，而且设计参数已经表格化，所以设计起来很方便。因此本文采用第一种方法。

设计步骤大致分以下三步：

(1) 设计模拟滤波器。根据实际需要确定滤波器的参数，利用的滤波器的设计公式设计出模拟滤波器并得到其传递函数H(s)，常用的滤波器有巴特沃斯滤波器、椭圆滤波器和切比雪夫滤波器；

(2) 将模拟滤波器转换成数字滤波器。利用冲激响应不变法或双线性变法将H(s)转换成H(z)，不同的设计方法对应于不同的s平面到z平面的映射公式；

(3) 频率变换。上述方法得到的是低通滤波器，为了得到高通、带通、带阻滤波器，还需要用利用变换公式作频率变换。

以上各步骤都有成熟的理论与公式，具体可参看文献[3]。

2.2均衡器的原理

2.2.1 均衡器总体设计

均衡器的基本功能是调节各频段的信号强弱，为了满足该功能，本文采用如下的方法：

Step1：设计出对应八个频段的八个带通滤波器；

Step2：对原始信号分八路用八个带通滤波器进行滤波；

Step3：将八个滤波器的滤波结果加权求和，权值的设计与均衡器的调节要求一致。

这样最终得到的结果便是所需要的均衡结果。其中第2步中各带通滤波器的输入信号均为原始信号，而不是“串联”地滤波。设原始输入信号为x(n)，第i路的输出信号为，第i路的权值为，均衡器的输出信号为y(n)，则有

                                                                     (5)

         (6)

式中，、为滤波器的参数，可由2.1.3节的方法得到，N为滤波器的阶数。

具体应用中，式(6)有两个问题，下节将介绍这两个问题并给出解决方法。

2.2.2 滤波计算的两个问题

(1) 数组越界问题

即当n<N时，如5阶滤波器处理第一个x时(n=1，N=5)，则式(6)中的与会出现数组下标小于0。

对于声音信号，人的听觉不会感觉出很短时间的异常。因此，可以简单地将数组越界的情况用0代替，即n<N时，规定。这样只有前N个数字信号不准确，从而前N/Fs秒声音信号不准确，N一般不超过10，而Fs一般为44100，从而N/Fs很小，人的听觉根本不会感觉到。

(2) 溢出问题

定点加法运算有可能发生溢出。但是，在采用补码进行运算时，即使中间计算结果发生了溢出，但只要最终累加值的绝对值小于1，就能保证最后得到的总和是正确数值。因此，为防止数字滤波器定点运算产生溢出， y(n)绝对值小于1就够了。

一种解决方式是针对输入x。对于IIR滤波器，有

                                (7)

由不溢出的条件得

                                                                       (8)

　

　　　

即，当(8)式不满足时，需要对输入信号x(n)乘以一个小于1的比例因子，使得(8)式恰好满足。比例因子的计算式易由(8)式推出。

为简单起见，本文采用针对输出的解决方式，即如果y(n)绝对值超过1，则对其赋值为1或，符号选择取决于原值的符号。

2.3 软件设计

2.3.1 数据流图

数据流图(DFD)是一种图形化技术，它描绘信息流和数据从输入移动到输出的过程中所经受的变换。在数据流图中没有任何具体的物理元素，它只是描绘信息在软件中流动和被处理的情况。设计数据流图时只需考虑系统必须完成的基本逻辑功能，而不用考虑具体实现，因而它是进行软件设计很好的出发点。

面向数据流的设计方法的目标是给出设计软件结构的一个系统化的途径。结合上一节的内容，可以得出软件的数据流图如图 2所示。

图 2 均衡器的数据流图

2.3.2 模块划分

模块化就是把程序划分成独立命名且可独立访问的模块，每个模块完成一个子功能，把这些模块集成起来构成一个整体，可以完成指定的功能满足用户需求。根据人类解决一般问题的经验，如果一个问题由两个问题组合而成，那么它的复杂程度大于分别考虑每个问题时的复杂程度之和，也就是说把复杂的问题分解成许多容易解决的小问题，原来的问题也就容易解决了。这就是模块化的根据。

在模块划分时应遵循如下规则^[4]：改进软件结构提高模块独立性；模块规模应该适中；深度、宽度、扇出和扇入都应适当；模块的作用域应该在控制域之内；力争降低模块接口的复杂程度；设计单入口单出口的模块；模块功能应该可以预测。

本着上述的启发式规则，对软件进行如图 3所示的模块划分。

图 3数字均衡器的模块划分

3．软件实现

3.1界面设计

MATLAB是Mathworks公司推出的数学软件，它将数值分析、矩阵计算、信号处理和图形显示结合在一起，为众多学科领域提供了一种简洁、高效的编程工具。它提供的GUIDE工具为可视化编程工具，使得软件的界面设计像VB一样方便。故本文采用MATLAB作为主要编程语言实现数字均衡器。

为了实现预期的功能，设计如图 4所示的界面。

图 4　频谱分析仪的界面设计

左边最上面的部分为标题区，用于显示软件标题等信息，不具人机交互功能。

再往下是信号输入/输出区，包含3种输入方式。界面应该具有：只有当每个单选框被选中时才允许使用对应的输入框、按钮等。支持WAV和AU两种格式的音频文件。“原声播放”用于播放输入的原信频信号，区别于滤波后的信号，用于对比滤波前后的信号。

再往下是均衡器区。八个滚动条对应八个带通滤波器的权值，在中间时为1，最下为0，最上为10。“类型”单选框用于选择巴特沃斯、切比雪夫等各种类型IIR滤波器。“重置”用于将各输入恢复默认。

“滤波”后才能“播放”，这是考虑到MATLAB对多线程机制支持得不好，不宜使用一边滤波一边播放的方法。由于滤波时间较长，故添加了进度显示的功能。条形图是许多音频软件使用的频谱输出方式。

右边的上面是滤波器，可以指定各种类型与参数。“生成滤波器”可以产生一个符合条件的滤波器，并在下面的图中显示滤波器的频率特征曲线。“滤波”用于对输入的音频信号用生成的滤波器滤波。

再往下则分别是时域波形与频谱分析的图形显示，均为实时显示。

在界面设计时为避免错误的操作，在某些按钮不能使用时必须将其设置为不可用，而它们能够使用时立即设置为可用。例如，没有打开音频文件时“滤波”按钮不可用，打开后“滤波”按钮可用；点击“滤波”前，“播放”不可用，点击“滤波”并完成滤波后，“播放”则变为可用。

3.2均衡器模块的实现

均衡器模块功能是生成带通滤波器并对原信号滤波。为了便于程序的扩充与修改，将表 1中的频段数据独立出来，如下

handles.fband=[20 100 200 500 1000 2000 4000 8000 16000];

均衡器的功能主要在“滤波”按钮的回调函数中实现，下面具体介绍

开始需要得到滤波器的参数：频带与阶数，如下

order=str2double(get(handles.order,'String'));

num=8;%8 filters totally

fband=handles.fband;

接着就是产生滤波器并滤波了，每次滤波后产生的y_i(n)如果都存到一个变量中，则需要8个这样的变量，对于输入信号量很大时将需要大量的内存空间，为此，本文采用累加的方式，这样只用一个变量去存储滤波结果。

nn=length(handles.y);

handles.yy=zeros(size(handles.y));

for i=1:num

if get(handles.butterworth,'Value')==1

[b a]=butter(order,2*fband(i:i+1)/handles.Fs);

elseif get(handles.cheby1,'Value')==1

[b a]=cheby1(order,0.5,2*fband(i:i+1)/handles.Fs);

elseif get(handles.cheby2,'Value')==1

[b a]=cheby2(order,20,2*fband(i:i+1)/handles.Fs);

elseif get(handles.ellip,'Value')==1

[b a]=ellip(order,0.5,20,2*fband(i:i+1)/handles.Fs);

else

errordlg('No filter type chosen or filter type error!');

end

eval(sprintf('k=get(handles.band%d,''Value'');',i));

y=(9^k-1)/8*qfilter(b,a,handles.y);

handles.yy=handles.yy+y;

end

guidata(hObject,handles);

在程序中有一句y=(9^k-1)/8*qfilter(b,a,handles.y)，作用是加权以便下一句的求和，其中k是从滚动条中得到的参数，从0到2，但需要的权值是从0到10，令w为需要的权，即，由于实际需要w为指数变化，故设

                                                                           (9)

代入w=0, k=0; w=1, k=1; w=10, k=2解得a=1/8，b=9，，从而有上述代码。

3.3 混合编程实现filter函数

MATLAB中提供的filter函数，但有两个问题，一个是运行速度较慢，另一个是不能解决2.2.2节的溢出问题。由于filter函数的执行速度是整个软件速度的瓶颈，故需要提高其速度。故这里采用VC实现filter函数，具体通过MATLAB与VC的混合编程来实现。这样能够较好解决那两个问题。

具体步骤为：

(1) MATLAB中运行mex –setup将VC选作编译器

(2) 在VC中新建一个DLL工程，添加如下代码

#include "mex.h"

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])

{   double *a,*b,*x,*y;

int M,n,i,j;

b=mxGetPr(prhs[0]); //获得指向矩阵的指针

a=mxGetPr(prhs[1]);

x=mxGetPr(prhs[2]);

M=mxGetM(prhs[2]); //获得矩阵的行数

plhs[0]=mxCreateDoubleMatrix(M,1,mxREAL);

y=mxGetPr(plhs[0]);

n=mxGetN(prhs[0]);

for(i=0;i<n;i++)

{   y[i]=b[0]*x[i];

for(j=1;j<i+1;j++)

y[i]+=b[j]*x[i-j]-a[j]*y[i-j];

if(y[i]>1)y[i]=1;else if(y[i]<-1)y[i]=-1;

}

for(;i<M;i++)

{   y[i]=b[0]*x[i];

for(j=1;j<n;j++)

y[i]+=b[j]*x[i-j]-a[j]*y[i-j];

if(y[i]>1)y[i]=1;else if(y[i]<-1)y[i]=-1;

}

}

(3) 保存为qfilter.cpp，放在MATLAB工作目录，在MATLAB中运行mex qfilter.cpp生成qfilter.dll。

这样，只需调用qfilter函数即能实现滤波。

3.4 实时显示的实现

软件的功能要求实时显示波形与频谱分析结果，这里采用分批处理的方法。即先取一批信号，将其播放出来，同时进行FFT，并绘出时域波形与频谱分析结果，然后取下一批信号……由于各批信号间的时间间隔很小，从而感觉像连续播放，而时域波形与频谱则是实时地变化。

该功能主要在“播放”按钮的回调函数中实现

N=length(handles.yy);

n=4096;

step=n/handles.Fs;

from=1;

to=n;

ymax=max(handles.y);

ymin=min(handles.y);

f=linspace(0,handles.Fs/2,n/2);

fband=handles.fband;

num=length(fband)-1;

tic

while to<=N;

y=handles.yy(from:to,:);

pause(step-toc*0.99);

wavplay(y,handles.Fs,'async');

tic

P=fft(y,n);

Pyy =2*abs(P)/n;

loglog(handles.axes1,f,Pyy(1:n/2));

axis(handles.axes1,[fband(1) fband(num+1) 1e-3 ymax]);

plot(handles.axes2,y);

ylim(handles.axes2,[ymin ymax]);

from=from+n;

to=to+n;

end

采用异步播放方式可以使得播放不占用主线程，从而使得播放效果不会很卡。为了使得播放流畅，需要与pause的使用结合，具体方法为用tic和toc计算出其它代码运行时间，然后用播放时间减去其它代码运行时间，从而得到需要暂停运行的时间。

4．运行实例与误差分析

4.1运行实例

图 5是程序的运行的一个实例，可以看到，基本上可以实现各个频段的均衡与调结。但是并不能完全滤掉某个频段的信号。

另外，为了比较MATLAB的filter函数与VC写的qfilter的执行速度，软件中提供了选择滤波程序的界面，图 5中的运行结果为qfilter的运行时间，为0.031秒，而相同情况下，MATLAB的filter的运行时间为0.157秒，如图 6所示。从而混合编程的方式大大提高了软件的运行效率。所以在均衡器功能上直接采用qfilter。

4.2 误差分析

理想的滤波器是不存在的，实际中只能尽量地接近理想滤波器，所以只能尽量地去接近理想的滤波器。对于IIR数字滤波器，阶数越高，则滤波器的精度越高，即越接近于理想的情况，但计算机运行速度就越慢了，所以实际中应该权衡运行速度与滤波器的精度。另外，计算机本身也有误差(如截断误差)，但相对于滤波器的设计误差，计算机的计算误差可以怱略不计。故提高精度关键在于提高滤波器的设计精度。

图 5 运行结果实例

图 6 MATLAB的filter的运行时间

另外，每个带通滤波器之间的过渡也有设计误差，理论上带通滤波器间的过滤应该满足各滤波器的频率响应曲线叠加起来为一条值为1的水平线，但实际中只能尽量去接近理论的情况。

5．总结与展望

5.1总结

本文首先对数字滤波器的原理进行了分析，并选择了一种实现数字均衡器的方案，在软件的设计上采用软件工程的理论，使用面向数据流的设计方法进行软件设计，并且采用了模块划分的方法。在实现的时候采用MATLAB与VC混合编程的方式，这样不仅利用了MATLAB的强大工具箱，而且使得开发的软件有着运行速度快的特点。界面上采用MATLAB提供的GUIDE工具进行设计，使得软件开发效率较高。对于均衡器所要求的功能，已经全部实现，但限于篇幅，本文只给出了部分功能的实现方法。最后给出了软件的运行结果，并对软件进行了误差分析。

5.2展望

该软件在播放前需要点滤波按钮，这使得操作变得麻烦，所以可以将滤波功能加到播放中，这就意味着需要加快程序的运行速度可者采用多线程机制。带通滤波器频率交界处的误差可以通过优化设计的方法来将其减小，即各带通滤波器不完全以给定的频带为自己的生成参数，而是以附近的某个频率值作为其生成参数，这样就需要用到优化的理论。其实可以将软件的主体用VC进行开发，部分难以实现但运行次数不多的函数则借用MATLAB，实现以VC的主体的混合编程，这样可以实现多线程，并且运行速度会进一步加快，但这种开发方式的开发时间相对较长。