基于MATLAB的数字均衡器设计

1、基于MATLAB的数字滤波器设计摘要:通过Matlab 强大的信号处理功能,分析数字均衡器的设计要求,对各种数字音频信号进行滤波处理,设计出一种比较合理的数字滤波器,并在此基础上设计一种均衡器,最后对该数字滤波器和均衡器进行综合测试并改进,使其达到要求的指标。关键词: MATLAB;数字滤波器;均衡器一、概述随着数字化技术的快速、深入发展,人们对数字化电子产品所产生的图像、图形以及声音等质量的要求越来越高。在实时数字处理过程中,与D/A 和A/D 转换相关的模拟信号重构过程是决定数字系统输出质量的关键。在声音的拾取过程及通过音响设备的传送过程中,由于设备或器件的原因,其幅度对频率的响应往往不一

2、致,这样就达不到原来的效果,往往需要对目标信号进行滤波处理,以满足用户对信号的要求。MATLAB 是一个数据分析和处理功能十分强大的工程实用软件,它的滤波器设计工具箱为实现声音信号的数字滤波提供了十分方便的函数和命令。本文将介绍基于MATLAB设计出的一种实用的数字滤波器,并对其功能进行扩展,设计出一种均衡器。滤波器的种类很多,按所通过信号的频段分为低通、高通、带通和带阻滤波器四种。低通滤波器:它允许信号中的低频或直流分量通过,抑制高频分量或干扰和噪声。高通滤波器:它允许信号中的高频分量通过,抑制低频或直流分量。带通滤波器:它允许一定频段的信号通过,抑制低于或高于该频段的信号、干扰和噪声。带阻

3、滤波器:它抑制一定频段内的信号,允许该频段以外的信号通过。上述每种滤波器又可以分为模拟滤波器和数字滤波器。如果滤波器的输入输出都是数字信号,则这样的滤波器称之为数字滤波器,它通常通过一定的运算关系改变输入信号所含频率成分的相对比例或者滤除某些频率成分来实现滤波。根据数字滤波器冲激响应的时域特性,可将数字滤波器分为两种,即无限长冲激响应(IIR滤波器和有限长冲激响应(FIR滤波器。根据数字信号处理的一般理论, IIR滤波器的特征是具有无限持续时间的冲激响应,而FIR滤波器使冲激响应只能持续一定的时间。从设计方法上讲,数字滤波器又可分为:巴特沃斯响应(最平坦响应巴特沃斯响应能够最大化滤波器的通带平

4、坦度。该响应非常平坦,非常接近DC信号,然后慢慢衰减至截止频率点为-3dB,最终逼近-20ndB/decade的衰减率,其中n为滤波器的阶数。巴特沃斯滤波器特别适用于低频应用,其对于维护增益的平坦性来说非常重要。贝塞尔响应除了会改变依赖于频率的输入信号的幅度外,滤波器还会为其引入了一个延迟。延迟使得基于频率的相移产生非正弦信号失真。就像巴特沃斯响应利用通带最大化了幅度的平坦度一样,贝塞尔响应最小化了通带的相位非线性。切贝雪夫响应在一些应用当中,最为重要的因素是滤波器截断不必要信号的速度。如果你可以接受通带具有一些纹波,就可以得到比巴特沃斯滤波器更快速的衰减。附录A 包含了设计多达8阶的具巴特沃

5、斯、贝塞尔和切贝雪夫响应滤波器所需参数的表格。其中两个表格用于切贝雪夫响应一个用于0.1dB 最大通带纹波;另一个用于1dB 最大通带纹波。二、设计原理根据数字滤波器冲激响应的时域特性的比较,FIR 数字滤波器优点是对有限字长效应不敏感,严格线性相位;但是实现同样指标,所需计算量远大于IIR 滤波器 。为了说明两种滤波器的设计方法,本设计中在数字滤波器部分使用FIR 方法,而在均衡器部分使用IIR 的方法分别进行设计。2.1 用窗函数法设计FIR 滤波器根据过渡带宽及阻带衰减要求,选择窗函数的类型并估计窗口长度N (或阶数M=N-1,窗函数类型可根据最小阻带衰减As 独立选择,因为窗口长度N

6、对最小阻带衰减As 没有影响,在确定窗函数类型以后,可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N ,设待求滤波器的过渡带宽为w ,它与窗口长度N 近似成反比,窗函数类型确定后,其计算公式也确定了,不过这些公式是近似的,得出的窗口长度还要在计算中逐步修正,原则是在保证阻带衰减满足要求的情况下,尽量选择较小的N ,在N 和窗函数类型确定后,即可调用MATLAB 中的窗函数求出窗函数wd (n 。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n,如果给出待求滤波器频率应为Hd ,则理想的单位脉冲响应可以用下面的傅里叶反变换式求出:在一般情况下,hd(n是不能用封闭公式表示的,需要

7、采用数值方法表示;从w=0到w=2采样N 点,采用离散傅里叶反变换(IDFT 即可求出。用窗函数wd(n将hd(n截断,并进行加权处理,得到如果要求线性相位特性, 则h(n还必须满足:根据上式中的正、 负号和长度N 的奇偶性又将线性相位FIR 滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如, 要设计线性相位低通特性可选择h(n=h(N-1-n一类,而不能选h(n=-h(N-1-n一类。验算技术指标是否满足要求,为了计算数字滤波器在频域中的特性,可调用freqz 子程序,如果不满足要求,可根据具体情况,调整窗函数类型或长度,直到满足要求为止。1(2j j n d d h n H

8、 e e d -=(d h n h n n =(1h n h N n =±-2.2 用双线性变换法设计IIR滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从S平面到Z平面是多值的映射关系所造成的。为了克服这一缺点,可以采用非线性频率压缩方法,将整个频率轴上的频率范围压缩到-/T/T之间,再用z=esT转换到Z平面上。也就是说,第一步先将整个S平面压缩映射到S1平面的-/T/T一条横带里;第二步再通过标准变换关系z=es1T将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系,消除了多值变换性,也就消除了频谱混叠现象,映射关系如图1所示。 图1

9、双线性变换的映射关系IIR数字滤波器的最通用的方法是借助于模拟滤波器的设计方法。模拟滤波器设计已经有了相当成熟的技术和方法,有完整的设计公式,还有比较完整的图表可以查询,因此设计数字滤波器可以充分利用这些丰富的资源来进行。为了将S平面的整个虚轴j压缩到S1平面j1轴上的-/T到/T段上,可以通过以下的正切变换实现 (1式中,T仍是采样间隔。当1由-/T经过0变化到/T时,由-经过0变化到+,也即映射了整个j轴再将S1平面通过以下标准变换关系映射到Z平面z=es1T从而得到S平面和Z平面的单值映射关系为: (2 (3式(2与式(3是S平面与Z平面之间的单值映射关系,这种变换都是两个线性函数之比,

10、因此称为双线性变换。对于IIR数字滤波器的设计具体步骤如下:(1按照一定的规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。(2根据转换后的技术指标设计模拟低通滤波器G(s(G(s是低通滤波器的传递函数。(3再按照一定的规则将G(s转换成H(z(H(z是数字滤波器的传递函数。若设计的数字滤波器是低通的,上述的过程可以结束,若设计的是高通、带通或者是带阻滤波器,那么还需要下面的步骤:将高通、带通或带阻数字滤波器的技术指标转换为低通模拟滤波器的技术指标,然后设计出低通G(s,再将G(s转换为H(z。Matlab信号工具箱提供了几个直接设计IIR数字滤波器的函数,直接调用这些函数就可以

11、很方便地对滤波器进行设计。这里选取巴特沃斯法、切比雪夫I、切比雪夫、椭圆法四种方法进行比较。三、设计步骤3.1模块设计模块设计如下图所示。 图2 模块设计3.2 GUI面板设计根据设计的模块,设计面板控件如下: 图3 GUI面板设计其中,在输入方式的电脑录音中,用户可以设定录音的采样频率和录音时间。均衡器是分为八段,分别是100 200 500 1000 2000 4000 8000 16000 20000。均衡器其实是多个带通滤波器输出的信号经过加权后混迭,其通带就是上面的各个数字之间的带宽。各带通滤波器的输入信号均为原始信号,而不是“串联”地滤波。设原始输入信号为x(n,第i 路的输出信号

12、为yi(n,第i 路的权值为ki,均衡器的输出信号为y(n,则有 (1对两种输入方式的信号采样分别编写一下程序,保存采样信号,以便后面滤波器或均衡器使用。导入wav文件:handles.y handles.fs bits=wavread('语音信号.wav'nn=length(handles.y;handles.yy=handles.y(fix(nn/2:fix(nn*3/4;axes(handles.axes1;plot(handles.axes1,handles.yy;title('滤波前波形'wavplay(handles.yy,handles.fs;电脑

13、录音:handles.fs=handles.soundf;handles.yy=wavrecord(handles.time*handles.soundf;wavplay(2*handles.yy,handles.soundf;plot(handles.axes1,handles.yy;title('滤波前波形'3.3用窗函数方法设计滤波器在MATLAB中,可以用b=fir1(N,Wn,ftype,taper 等函数辅助设计FIR 数字滤波器。N代表滤波器阶数;Wn代表滤波器的截止频率(归一化频率,当设计带通和带阻滤波器时,Wn为双元素相量;ftype代表滤波器类型,如high

14、高通,stop带阻等;taper为窗函数类型,默认为海明窗,窗系数需要实现用窗函数blackman, hamming,hanning chebwin, kaiser产生。本设计中采用的是kaiser窗进行滤波。由于滤波器有四种功能,每种滤波所用的参数都不一样,故用下面的一段程序,对用户的选择进行判断,调用相应的程序段,实现相应的滤波功能。其中低通滤波器的截止频率是3200Hz,带通滤波器的带宽是2000-8000Hz,带阻的带宽也是2000-8000Hz,高通的截止频率是5000Hz。Y1=fft(handles.yy;switch handles.typecase 1 %低通fp=3000;

15、fc=3200;As=100;Ap=1;Fs=17000;wc=2*fc/Fs; wp=2*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7;N=ceil(As-8/2.285/wdel;ws=(wp+wc/2/pi;wn= kaiser(N+1,beta;b=fir1(N,ws,wn;case 2 %带通fp1=2200 ;fp2=8000 ;fc1=2000 ;fc2=8200 ;As=100 ;Ap=1 ;Fs=17000 ;wp1=2*fp1/Fs; wc1=2*fc1/Fs; wp2=2*fp2/Fs; wc2=2*fc2/Fs;wdel=wp1-wc1;bet

16、a=0.112*(As-8.7;N=ceil(As-8/2.285/wdel;ws =(wp1+wc1/2/pi,(wp2+wc2/2/pi;wn= kaiser(N+1,beta;b=fir1(N,ws,wn;case 3%带阻fp1=2200 ;fp2=8000 ;fc1=2000 ;fc2=8200 ;As=100 ;Ap=1 ;Fs=17000 ;wp1=2*fp1/Fs; wc1=2*fc1/Fs; wp2=2*fp2/Fs; wc2=2*fc2/Fs;wdel=wp1-wc1;beta=0.112*(As-8.7;N=ceil(As-8/2.285/wdel;ws =(wp1+wc

17、1/2/pi,(wp2+wc2/2/pi;wn= kaiser(N+1,beta;b=fir1(N,ws,'stop',wn;case 4%高通fp=4800;fc=5000;As=100;Ap=1;Fs=17000;wc=2*fc/Fs; wp=2*fp/Fs;wdel=wc-wp;beta=0.112*(As-8.7;N=ceil(As-8/2.285/wdel;wn= kaiser(N+1,beta;ws=(wp+wc/2/pi;b=fir1(N,ws,'high',wn;下面以导入wav文件的输入方式为例,说明程序执行的结果。1、截止频率是3200Hz的

18、低通滤波器利用freqz函数查看低通滤波器的响应特性,如下图。 图4 低通滤波器相应特性进行低通滤波后波形及频谱显示: 图5 低通滤波后波形和频谱显示2、带宽是2000-8000Hz的带通滤波器相应特性: 图6 带通滤波器响应特性进行带通滤波后波形及频谱显示: 图7 带通滤波后波形和频谱显示带阻和高通滤波器结果就不再赘述了,其结果可参看附带程序文件。3.3用双线性变换法设计均衡器本设计中的均衡器有八段,分别用八个滚动条来表示各个频率段经过带通滤波后信号幅值的加权值大小(范围是0-10。用户可选取四种分别基于巴特沃斯法、切比雪夫I、切比雪夫、椭圆法的均衡器。用户的选择功能可用一个下拉菜单来实现。

19、对于滤波器的阶数本应该是有一定的函数进行计算的,本设计中由于各个滤波器的频率大小有很大的跨越,经过计算和实验,得出了满足要求的阶数的范围,为1-4Y1=fft(handles.yy;handles.fband =100 200 500 1000 2000 4000 8000 16000 20000; handles.Fs=80000;num=8;fband=handles.fband;y=zeros(size(handles.yy;N=handles.N;for i=1:num;if get(handles.popupmenu1,'Value'=1;b a=butter(N,2*fband(i:i+1/handles.Fs;elseif get(handles.popupmenu1,'Value'=2;b a=c