基于Matlab的IIR数字滤波器的设计别人实验.doc

基于MatlabIIR数字滤波器设计
一 绪论
数字滤波是数字信号处理重要内容，是由乘法器、加法器和单位延时器组成一种运算过程，其功能是对输入离散信号进行运算处理，以达到改变信号频谱目。数字滤波器根据频域特性可分为低、高通、带通和带阻四个基本类型；根据时域特性可分为无限长单位冲激响应IIR(Infinite Impulse Response)滤波器和有限长单位冲激响应FIR(Finite Impulse Response)滤波器。
数字滤波在通信、图像编码、语音编码、雷达等许多领域中有着十分广泛应用。目前，数字信号滤波器设计图像处理、数据压缩等方面应用取得了令人瞩目进展和成就。鉴于此，数字滤波器设计就显得尤为重要。
MATLAB是美国MathWorks公司推出一套用于工程计算可视化高性能语言及软件环境。MATLAB为数字滤波研究和应用提供了一个直观、高效、便捷利器。它以矩阵运算为基础，把计算、可视化、程序设计融合到了一个交互式工作环境中。MATLAB推出工具箱使各个领域研究人员可以直观方便地进行科学研究、工程应用，其中信号处理(signal processing)、图像处理(image processing),小波(wavelet)等工具箱为数字滤波研究蓬勃发展提供了有力工具。
二 数字滤波器
2.1什么是数字滤波器
滤波器是指用来对输入信号进行滤波硬件和软件。所谓数字滤波器，是指输入、输出均为数字信号，通过一定运算关系改变输入信号所含频率成分相对比例或者滤除某些频率成分器件。数字滤波器和模拟滤波器相比，因为信号形式和实现滤波方法不同，数字滤波器具有比模拟滤波器精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配等优点。
一般用两种方法来实现数字滤波器:一是采用通用计算机，把滤波器所要完成运算编成程序通过计算机来执行，也就是采用计算机软件来实现;二是设计专用数字处理硬件。
MATLAB信号处理工具箱是专门应用于信号处理领域专用工具箱，它两个基本组成就是滤波器设计及实现部分以及谱分析部分。工具箱提供了丰富而简便设计，使原来繁琐程序设计简化成函数调用。只要以正确指标参数调用相应滤波器设计程序或工具箱函数，便可以得到正确设计结果，使用非常方便。
2.2数字滤波器分类
数字滤波器从功能上分类:可分为低通滤波器、高通滤波器、带通滤波器、带阻滤波器。
从滤波器网络结构或者从单位脉冲响应分类:可分为IIR滤波器(即无限长单位冲激响应滤波器)和FIR滤波器(即有限长单位冲激响应滤波器)。它们函数分别为:

第一个公式中H (z)称为N阶IIR滤波器函数，第二个公式中H (z)称为(N-1)阶FIR滤波器函数。
2.3数字滤波器设计要求
滤波器指标常常在频域给出。数字滤波器频响特性函数一般为复函数，所以通常表示为：

其中，||称为幅频特性函数，Φ(w)称为相频特性函数。幅频特性表示信号通过该滤波器后各频率成分衰减情况，而相频特性反映各频率通过滤波器后在时间上延时情况。一般IIR数字滤波器，通常只用幅频响应函数||来描述设计指标，相频特性一般不作要求。
IIR滤波器指标参数如下图所示。图中，ωp和ωs分别为通带边界频率和阻带边界频率；δ1和δ2分别为通带波纹和阻带波纹；允许衰减一般用dB数表示，通带内所允许最大衰减（dB）和阻带内允许最小衰减（dB）分别为αp和αs表示：

一般要求：

低通滤波器技术要求

2.4数字滤波器设计方法概述
IIR数字滤波器设计步骤流程图如下：

步骤流程图
IIR滤波器设计方法有两类，经常用到一类设计方法是借助于模拟滤波器设计方法进行。其设计思路是:先设计模拟滤波器得到传输函数Ha(s)，然后将Ha(s)按某种方法转换成数字滤波器系统函数H (Z)。这一类方法是基于模拟滤波器设计方法相对比较成熟，它不仅有完整设计公式，也有完整图表供查阅。更可以直接调用MATLAB中对应函数进行设计。另一类是直接在频域或者时域中进行设计，设计时必须用计算机作辅助设计，直接调用MATLAB中一些程序或者函数可以很方便地设计出所需要滤波器。

三 IIR滤波器设计
3.1典型IIR数字滤波器设计
模拟滤波器理论和设计方法已发展得相当成熟，且有一些典型模拟滤波器供我们选择，如巴特沃斯(Butterworth )滤波器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer )滤波器、贝塞尔(Bessel )滤波器等，这些典型滤波器各有特点。
用MATLAB进行典型数字滤波器设计，一般步骤如下:
（1） 将设计指标归一化处理。如果采用双线性变换法，还需进行预畸变。
（2） 根据归一化频率，确定最小阶数N和频率参数Wn。可供选用阶数择函数有:buttord，cheblord，cheb2ord，ellipord等。
（3） 运用最小阶数N设计模拟低通滤波器原型。模拟低通滤波器创建函数有:buttap,cheblap, cheb2ap,ellipap和besselap，这些函数输出是零极点式形式，还要用zp2tf函数转换成分子分母多项式形式。如果想根据最小阶数直接设计模拟低通滤波器原型，可用butter,chebyl,cheby2,ellip,bessel等函数，只是注意要将函数中Wn设为1。
（4） 根据第2步频率参数Wn，模拟低通滤波原型转换模拟低通、高通、带通、带阻滤波器，可用函数分别是：lp21p,lp2hp,lp2bp,lp2bs。
（5） 运用脉冲响应不变法或双线性变法把模滤波器转数字滤波器，调用函数是impinvar和bilinear。脉冲响应不变法适用于采样频率大于4倍截止频率锐截止低通带通滤波器，而双线性变换法适合于相位特性要求不高各型滤波器。
（6） 根据输出分子分母系数，用tf函数生成H（z）表达式，再用freqz函数验证设计结果。
设计巴特沃思数字低通滤波器和椭圆数字低通滤波器，要求通带边界频率fp=2.1kHZ，通带最大衰减Rp=0.5dB；阻带边界频率fs=8kHZ，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。

用脉冲响应不变法设计巴特沃思数字低通滤波器M程序如下：
fp=2100;
fs=8000;
Fs=20000;
Rp=0.5;
Rs=30;
T=1/Fs; %设计指标
W1p=fp/Fs*2;W1s=fs/Fs*2;%求归一化频率
[N,Wn]=buttord(W1p,W1s,Rp,Rs,'s');
%确定butterworth最小介数N和频率参数Wn
[z,p,k]=buttap(N); %设计模拟低通原型零极点增益参数
[bp,ap]=zp2tf(z,p,k); %将零极点增益转换成分子分母参数
[bs,as]=lp2lp(bp,ap,Wn*pi*Fs);%将低通原型转换为模拟低通
[bz,az]=impinvar(bs,as,Fs); %用脉冲响应不变法进行模数变换
sys=tf(bz,az,T); %给出传输函数H(Z)
[H,W]=freqz(bz,az,512,Fs); %生成频率响应参数
subplot(2,1,1);
plot(W,20*log10(abs(H))); %绘制幅频响应
grid on; %加坐标网格
xlabel('频率/Hz');
ylabel('振幅/dB');
subplot(2,1,2);
plot(W,abs(H)); grid on;
xlabel('频率/Hz');
ylabel('振幅/H');

运行后波形如下：

运行结果：
N =4
bz = 0.0000 0.0999 0.1914 0.0252
az= 1.0000 -1.4336 1.0984 -0.4115 0.0627

用双线性变换法设计椭圆数字低通滤波器M程序如下：
fs=20000;
wp=2*pi*2100/fs;
ws=2*pi*8000/fs;
Rp=0.5;
Rs=30;
Ts=1/fs;
Wp=2/Ts*tan(wp/2);Ws=2/Ts*tan(ws/2); %按频率转换公式进行转换
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器最小阶数
[z,p,k]=ellipap(N,Rp,Rs);%设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器频率转化
[bz,az]=bilinear(b,a,fs); %运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,512,fs);
subplot(2,1,1);
plot(f,20*log10(abs(H)));
title('N=2 频率响应');
grid on;
xlabel('频率/Hz');
ylabel('振幅/dB');
subplot(2,1,2);
plot(f,abs(H)); grid on;
xlabel('频率/Hz');
ylabel('振幅/H');

运行结果：
N=2
bz= 0.1213 0.1662 0.1213
az= 1.0000 -0.9889 0.4218

3.2完全滤波器设计
除了典型设计以外，MATLAB信号处理工具箱提供了几个直接设计IIR数字滤波器函数，直接调用就可以设计滤波器，这为设计通用滤波器提供了方便。
设计Butterworth滤波器用函数butter()，可以设计低通、高通、带通和带阻数字和模拟滤波器，其特性是通带内幅度响应最大限度平滑，但损失了截止频率处下降斜度。
设计Chebyshev I型滤波器用函数chebyl()。可以设计低通、高通、带通和带阻数字和模拟Chebyshev I型滤波器，其通带内为等波纹，阻带内为单调。Chebyshev I型滤波器下降斜度比II型大，但其代价目是通带内波纹较大。
设计Chebyshev II型滤波器用函数cheby2()。可以设计低通、高通、带通和带阻数字和模拟Chebyshev II型滤波器，其通带内为单调，阻带内等波纹。Chebyshev II型滤波器下降斜度比I型小，但其阻带内波纹较大。
设计椭圆滤波器用函数ellip()，及chebyl, cheby2类似，可以设计低通、高通、带通和带阻数字和模拟滤波器。及Butterworth和chebyshev滤波器相比，ellip函数可以得到下降斜度更大滤波器，得通带和阻带均为等波纹。一般情况下，椭圆滤波器能以最低阶实现指定性能指标。
在使用各类滤波器函数时应当注意以下重点:
A、阶数和固有频率选择:[N,Wn]=buttord(Wp,Ws,Rp,Rs)可得到符合要求性质滤波器最小阶数N以及数字Butterworth滤波器固有频率Wn(即3dB )。设计要求是在通带内衰减不超过Rp,在阻带内衰减不小于Rs，通带和阻带有截止频率分别是Wp, Ws,它们是归一化频率，范围是[0, 1]，对应π弧度。
B、有关滤波器设计当中频率归一化问题:信号处理工具箱中经常使用频率是Nyquist频率，它被定义为采样频率一半，在滤波器阶数选择和设计中截止频率均使用Nyquist频率进行归一化处理。例如对于一个采样频率为1000 Hz系统，400Hz归一化即为400/500=0.8。归一化频率范围在[0, 1]之间。如果要将归一化频率转换为角频率，则将归一化频率乘以π;如果要将归一化频率转换为Hz,则将归一化频率乘以采样频率一半。
C、设计一个N阶低通Butterworth滤波器使用函数[B,A]=butter(N, Wn)，返回滤波器系数矩阵[B,A]。其中固有频率Wn必须是归一化频率。它最大值是采样频率一半。Fs缺省时默认为2Hz。如果Wn=[Wl,W2]是一个两元素向量，则函数将设计出一个2N阶带通滤波器，通带为[W1,W2]。
设计Chebyshev I型和Chebyshev II型数字低通滤波器，要求通带边界频率fp=2.1kHZ，通带最大衰减Rp=0.5dB；阻带边界频率fs=8kHZ，阻带最小衰减Rs=30dB,采样频率为Fs=20kHZ。
Chebyshev I型M程序如下：
Fs=20000; %抽样频率20KHz
Flp=2100;
Fls=8000;
Wp=2*Flp/Fs; %归一化通带截止频率
Ws=2*Fls/Fs; %归一化阻带截止频率
Rp=0.5; %通带最大衰减（单位：dB）
Rs=30; %阻带最小衰减（单位：dB）
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率
[b,a]=cheby1(N,Rp,Wn); %返回H(z)分子分母系数
[hw,w]=freqz(b,a);
subplot(2,1,1);
plot(w/pi,20*log10(abs(hw)));grid on;
xlabel('ω/π');ylabel('幅度（dB）')
title('切比雪夫I型幅频响应');
subplot(2,1,2);plot(w/pi,abs(hw));
grid on;
xlabel('ω/π');ylabel('幅度（H）');
运行后波形如下：

运行结果：
N=2
b= 0.1007 0.2014 0.1007
a= 1.0000 -0.9872 0.4140

Chebyshev II型M程序如下：
Fs=20000; %抽样频率20KHz
Flp=2100;
Fls=8000;
Wp=2*Flp/Fs; %归一化通带截止频率
Ws=2*Fls/Fs; %归一化阻带截止频率
Rp=0.5; %通带最大衰减（单位：dB）
Rs=30; %阻带最小衰减（单位：dB）
[N,Wn]=cheb2ord(Wp,Ws,Rp,Rs); %返回最小阶数和截止频率
[b,a]=cheby2(N,Rs,Wn); %返回H(z)分子分母系数
[hw,w]=freqz(b,a);
subplot(2,1,1);
plot(w/pi,20*log10(abs(hw)));grid on;
xlabel('ω/π');ylabel('幅度（dB）')
title('切比雪夫II型幅频响应');
subplot(2,1,2);plot(w/pi,abs(hw));
grid on;
xlabel('ω/π');ylabel('幅度（H）');
运行后波形：

运行结果：
N=2
b= 0.2357 0.4241 0.2357
a= 1.0000 -0.2996 0.1950

3.3 结果分析
从频率响应图中可以看出:巴特沃斯滤波器具有单调下降幅频特性，通带内平滑;切比雪I型滤波器幅频特性在通带内有波动，阻带内单调;chebyshev II型滤波器幅频特性在阻带内有波动，通带内单调;椭圆滤波器选择性相对前三种是最好，下降斜度比较大，通带和阻带内均为等波纹，同样性能指标，椭圆滤波器可以最低阶数来实现。这样根据不同要求可以选用不同类型滤波器。
四 对语音信号滤波
利用函数wavread对语音信号采集：
i=1;
[x,fs,bits]=wavread('1.wav');
%x：语音数据；fs：采样频率；bits：采样点数
sound(x,fs,bits); %话音回放
N=length(x);
n=0:N-1;
figure(i);
subplot(2,1,1);
plot(n,x); %画出原始语音信号波形
xlabel('n');
ylabel('x(n)');
title('原始语音信号');
subplot(2,1,2);
[H,f]=freqz(x,1,512,fs);
plot(f,20*log10(abs(H))); %画出原始语音信号频谱
xlabel('n');
ylabel('x(n)');
title('原始语音信号频谱');

运行结果：
fs=44100
bits=16
加入上面设计椭圆滤波器对其进行滤波：
fs=20000;
wp=2*pi*2100/fs;
ws=2*pi*8000/fs;
Rp=0.5;
Rs=30;
Ts=1/fs;
Wp=2/Ts*tan(wp/2);Ws=2/Ts*tan(ws/2); %按频率转换公式进行转换
[N,Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %计算模拟滤波器最小阶数
[z,p,k]=ellipap(N,Rp,Rs);%设计模拟原型滤波器
[Bap,Aap]=zp2tf(z,p,k); %零点极点增益形式转换为传递函数形式
[b,a]=lp2lp(Bap,Aap,Wn); %低通转换为低通滤波器频率转化
[bz,az]=bilinear(b,a,fs); %运用双线性变换法得到数字滤波器传递函数
[H,f]=freqz(bz,az,512,fs);
plot(f,20*log10(abs(H)));
title('N=2 频率响应');
grid on;
xlabel('频率/Hz');
ylabel('振幅/dB');

%-------滤波前后波形比较-----%
i =i+1;
figure(i);
y=filter(bz,az,x);
subplot(2,2,1);
plot(n,x); %画出原语音信号波形
xlabel('n');
ylabel('x(n)');
title('原始语音信号');
subplot(2,2,3);
plot(n,y); %IIRLowPass后波形
title('IIRLowPass后波形');
%---------滤波前后语音信号频谱-----%
subplot(2,2,2);
[H,f]=freqz(x,1,512,fs);
plot(f,20*log10(abs(H))); %滤波前语音信号频谱
xlabel('频率/Hz');
ylabel('振幅/dB');
title('IIR低通滤波前语音信号频谱');
grid on;
subplot(2,2,4);
[H,f]=freqz(y,1,512,fs);
plot(f,20*log10(abs(H))); %滤波后语音信号频谱
xlabel('频率/Hz');
ylabel('振幅/dB');
title('IIR低通滤波后语音信号频谱');
grid on;
sound(y,fs,bits);

在matlab上运行后，发现语音信号在经过上面设计椭圆数字低通滤波器后，声音发生了明显改变。观察波形，因为导入语音信号波形频率分布绝大多数分布在低频区，故没有被滤波器过滤掉多少，所以语音信号波形变化不大。
但频谱发生了明显得改变。
同理，只要把之前设计其他三个滤波器程序替换掉如上椭圆滤波器程序，就可得到语音信号在不同滤波器滤波后波形变化情况。

参考文献
[1]赵树杰等，数字信号处理，西安电子科技大学出版社，1997.10
[2] 陈怀琛等，MATLAB及在电子信息课程中应用，电子工业出版社出版，2002.4
[3]李丽 王振领，MATLAB工程计算及应用，人民邮电出版社，2001.9
[4] 陈怀琛等译，《数字信号处理及其MATLAB实现》，电子工业出版社；
18 / 18