一、 实验目的
1. 掌握利用双线性变换法设计IIR数字滤波器的原理和具体方法。 2. 深入理解利用双线性变换法设计IIR数字滤波器的优缺点及使用范围。 二、 实验内容
1. 利用巴特沃思模拟滤波器,通过双线性变换法设计巴特沃思数字滤波器,数字滤波器的技
术指标为
采样周期为T=1。 程序代码
T=2; %设置采样周期为2 fs=1/T; Wp=0.25*pi/T;
Ws=0.35*pi/T; Ap=20*log10(1/0.9);
As=20*log10(1/0.18); [B,A]=butter(N,Wc,'s'); W=linspace(0,pi,400*pi); hf=freqs(B,A,W); subplot(2,1,1);
plot(W/pi,abs(hf)/abs(hf(1))); grid on;
title('巴特沃斯模拟滤波器'); xlabel('Frequency/Hz'); ylabel('Magnitude');
[D,C]= bilinear (B,A,fs); Hz=freqz(D,C,W); subplot(2,1,2);
plot(W/pi,abs(Hz)/abs(Hz(1))); %绘出巴特沃斯数字低通滤波器的幅频特性曲线 grid on;
title('巴特沃斯数字滤波器'); xlabel('Frequency/Hz');
1
%采样频率为周期倒数
%设置归一化通带和阻带截止频率 %设置通带最大和最小衰减 %调用butter函数设计巴特沃斯滤波器 %指定一段频率值
%计算模拟滤波器的幅频响应
%绘出巴特沃斯模拟滤波器的幅频特性曲线
[N,Wc]=buttord(Wp,Ws,Ap,As,'s'); %调用butter函数确定巴特沃斯滤波器阶数
%调用双线性变换法
%返回频率响应
ylabel('Magnitude');
进行试验
输出图像：
巴特沃斯模拟滤波器1Magnitude0.5000.10.20.30.40.50.60.7Frequency/Hz巴特沃斯数字滤波器0.80.911.5Magnitude10.5000.10.20.30.40.50.6Frequency/Hz0.70.80.91 实验分析
通过查workspace可以得出脉冲响应为：
滤波器幅度响应：滤波器相位响应：
2. 利用巴特沃思模拟滤波器，通过双线性变换法设计数字带阻滤波器，数字滤波器的技术
指标为
2
采样周期为T=1。 程序代码
T=1; %设置采样周期为1 fs=1/T; %采样频率为周期倒数 wp=[0.25*pi,0.75*pi]; ws=[0.35*pi,0.65*pi]; Wp=(2/T)*tan(wp/2);
Ws=(2/T)*tan(ws/2); %设置归一化通带和阻带截止频率 Ap=20*log10(1/0.9);
As=20*log10(1/0.18); %设置通带最大和最小衰减
[N,Wc]=buttord(Wp,Ws,Ap,As,'s');%调用butter函数确定巴特沃斯滤波器阶数 [B,A]=butter(N,Wc, 'stop','s'); %调用butter函数设计巴特沃斯滤波器 W=linspace(0,2*pi,400*pi); %指定一段频率值
hf=freqs(B,A,W); %计算模拟滤波器的幅频响应 subplot(2,1,1); plot(W/pi,abs(hf)); grid on;
title('巴特沃斯模拟滤波器'); xlabel('Frequency/Hz'); ylabel('Magnitude');
[D,C]=bilinear(B,A,fs); %调用双线性变换法 Hz=freqz(D,C,W); %返回频率响应 subplot(2,1,2);
plot(W/pi,abs(Hz)); %绘出巴特沃斯数字带阻滤波器的幅频特性曲线 grid on;
title('巴特沃斯数字滤波器'); xlabel('Frequency/Hz'); ylabel('Magnitude');
进行试验
输出图像：
%绘出巴特沃斯模拟滤波器的幅频特性曲线
3
巴特沃斯模拟滤波器1Magnitude0.5000.20.40.60.811.21.4Frequency/Hz巴特沃斯数字滤波器1.61.821Magnitude0.5000.20.40.60.811.2Frequency/Hz1.41.61.82 实验分析
通过查workspace可以得出脉冲响应为：
滤波器幅度响应：滤波器相位响应：
3. 通过切比雪夫模拟滤波器，利用双线性变换法设计数字高通滤波器，数字滤波器的技术
指标为
采样周期为T=1，并观察T分别为0.01，0.1，0.25，0.5，0.8时数字滤波器的频率响应。
4
程序代码
T=1; fs=1/T; Wp=(2/T)*tan(0.75*pi/2);
Ws=(2/T)*tan(0.65*pi/2); %设置归一化通带和阻带截止频率 Ap=20*log10(1/0.9);
As=20*log10(1/0.18); %设置通带最大和最小衰减
[N,Wc]=cheb1ord(Wp,Ws,Ap,As,'s');%调用cheb1ord函数确定切比雪夫I型滤波器阶数 [B,A]=cheby1(N,Ap,Wc, 'high','s'); %调用cheby1函数设计切比雪夫I型滤波器 W=linspace(0,7*pi,400*pi); %指定一段频率值
hf=freqs(B,A,W); %计算模拟滤波器的幅频响应 subplot(2,1,1);
plot(W/pi,abs(hf)); %绘出切比雪夫I型模拟滤波器的幅频特性曲线 grid on;
title('切比雪夫I型模拟滤波器'); xlabel('Frequency/Hz'); ylabel('Magnitude');
w=linspace(0,pi,400*pi); %指定一段频率值
[D,C]=bilinear(B,A,fs); %调用双线性变换法 Hz=freqz(D,C,w); %返回频率响应 subplot(2,1,2); plot(w/pi,abs(Hz)); grid on;
title('切比雪夫I型数字滤波器'); xlabel('Frequency/Hz'); ylabel('Magnitude');
进行试验 输出图像如下：
%绘出切比雪夫I型数字带阻滤波器的幅频特性曲线
%设置采样周期为1
%采样频率为周期倒数
5
切比雪夫I型模拟滤波器1Magnitude0.50012345Frequency/Hz切比雪夫I型数字滤波器671Magnitude0.5000.10.20.30.40.50.6Frequency/Hz0.70.80.91 实验分析
通过查workspace可以得出脉冲响应为：
滤波器幅度响应：滤波器相位响应：
当T分别为0.01，0.1，0.25，0.5，0.8时的数字滤波器的频率响应如下图像所示：
6
切比雪夫I型数字滤波器10.90.80.7Magnitude0.60.50.40.30.20.1000.10.20.30.40.50.6Frequency/Hz0.70.80.91
其中， T=0.01
由图像可以看出，当采样周期从0.01s增长到0.8s时，数字滤波器的幅度响应没有任何
T=0. 1 T=0.25 T=0.5 T=0.8
“ ” “-------------” “—·—·—” “————” “ ”
变化。故对于用双线性变换法设计的切比雪夫高通数字滤波器，由于事先做了预畸处理，消除了频率映射时引起的非线性关系，所得数字滤波器的幅度响应较为稳定。
对比用双线性变换法和脉冲响应不变法得到的数字高通滤波器幅度响应形状可知：双线性变换法不会出现由于高频部分超过折叠频率而混叠到低频部分去的现象。这意味着，模拟滤波器的全部频率响应被压缩于等效的数字频率范围（0，π）之内。这是双线性变换法的主要优点，虽然这是以频率的严重非线性关系为代价的，但因为可以用预畸的方法进行补偿，所以仍可以得到比较理想的数字滤波器幅度响应
7
各种设计方法的MATLAB实现
下面以一带通数字滤波器的设计为例,说明基于MATLAB的三种实现方法。 3.1 基于模拟低通原型的MATLAB实现
%通过模拟低通原型滤波器进行数字带通设计程序： fp=480 %模拟低通通带上限频率 fs=520 %模拟低通阻带下限频率
wp=2*pi*fp %模拟低通通带上限角频率 ws=2*pi*fs %模拟低通阻带下限角频率 rp=3 %通带波动 rs=20 %阻带衰减
%巴特沃斯模拟低通原型滤波器设计 [n,wn]=buttord(wp,WS,rp ,rs,'s')
[z,p,k]=buttap(n) %模拟低通原型零、极点系数
[b1,a1]=zp2tf(z,p,k) %零、极点系数转换为传递函数 %巴特沃斯模拟低通原型滤波器频率响应 [hl,w1]=freqs(b1,a1) mag1=abs(h1)
%模拟低通原型滤波器幅频特性曲线 subplot(221);semilogx(w1,mag1) fw=40 %模拟带通滤波器带宽频率
bw=2*pi*fw %模拟带通滤波器带宽角频率 % 由模拟低通原型变换为模拟带通滤波器
[b2,a2]=lp2bp(b1,a1,wn,bw) %模拟带通滤波函数系数 %巴特沃斯模拟带通滤波器频率响应 [h2,w2]= freqs(b2,a2) mag2 =abs(1l2)
%模拟带通滤波器幅频特性曲线(db) subplot(222);plot(20*log10(mag2)) % 冲击响应不变法进行离散化设计 fo=2000 %采样频率
[bz,az]=impinvar(b2,a2,2000) %数字带通滤波函数系数 %巴特沃斯型数字带通滤波器频率响应 [hz,w]= freqz(bz,az) magz=abs(hz)
phz=unwrap(angle(hz))
subplot(223);plot(magz) %数字带通滤波器幅频特性曲线 subplot(224);flot(plot) %数字带通滤波器相频特性曲线
8
3.2 基于合适类型模拟滤波器的MATLAB实现 %通过合适类型模拟滤波器进行数字带通设计程序 fp= [480,520];fs=[450,550] %模拟通带、阻带频率 wp=[480,520]*pi*2 %模拟通带角频率 ws=[450,550]*pi*2 %模拟阻带角频率 rp=3;rs=20 %通带波动、阻带衰减 % 巴特沃斯型模拟带通滤波器设计 [n,wn]=buttord (wp,ws,rp,rs,'s')
[b,a]=butter(n,wn,'s') %模拟带通滤波函数系数 %巴特沃斯型模拟带通滤波器频率响应 [ha,w]= freqs(b,a)
ma=abs(ha);pha=unwrap(angle(ha))
subplot(421);plot(w/(2*pi),ma) %模拟幅频曲线 subplot(423);plot(w/(2 pi),pha) %模拟相频曲线 % 冲击响应不变法进行离散化设计 fo=5000 %采样频率
[bn,an]=impinvar(b,a,5000) %数字带通滤波函数系数 %巴特沃斯型数字带通滤波器频率响应 [hz,w]=freqz(bn,an)
mz=abs(hz);phz=unwrap(angle(hz))
subplot(422);plot(w,mz) %数字滤波器幅频曲线 subplot(424);plot(w,phz) %数字滤波器相频曲线 hi=impz(bn,an) %数字滤波器冲击响应 subplot(425),plot(hi) %冲击响应曲线 n=0:300;t=n/fo
xl=2*square(2*pi*500*t) %500Hz方波信号 subplot(426);plot(x1) %500Hz方波波形 yi=conv(hi,x1) %时域卷积输出
subplot(427);plot(yi) %卷积输出波形 y1=filter(bn,an,x1) %数字滤波函数输出
subplot(428);plot(y1) %数字滤波器输出波形
9
3.3 基于直接原型变换法的MATLAB实现 %数字带通滤波器直接设计程序
fp= [480,520];fs=[450,550] %模拟通带、阻带频率 rp=3;rs=20 %通带波动、阻带衰减 fo=10000 %采样频率 %频率指标变换
wp=2*pi*fp/f0 %数字通带频率 ws=2*pi*fs/fo %数字阻带频率
%切比雪夫1型数字带通滤波器直接设计 [n,wn]=cbeblord(wp/pi,ws/pi,rp,rs)
[b,a]=chebyl(n,rp,wn) %数字带通滤波器系数 %切比雪夫1型数字带通滤波器频率响应 [h,w]= freqz(b,a,128,10000)
mag=abs(h;pha=unwrqp(angle(h)) subplot(321);plot(w,mag) %幅频曲线 subplot(322);plot(w,pha) %相频曲线 hi=impz(b,a) %冲击响应
subplot(324);plot(hi) %响应曲线 n=0:500;t=n/fc
x1=2*square(2*pi*500*t) %500Hz方波信号 subplot(323);plot(t,x1) %500Hz方波波形 yi=conv(hi,x1) %时域卷积输出
subplot(326);plot(yi) %卷积输出波形
10
y1=filter(b,a,x1) %数字滤波函数输出
subplot(325); stem(y1) %数字滤波器输出波形
4 结束语
基于MATLAB的信号处理工具箱为数字滤波器设计带来了全新的实现手段,设计快捷方便,仿真波形直观。上述三种设计方案均可实现设计指标,但以直接原型变换法最为简便。实际应用中,数字滤波器也可以对连续时间信号进行处理,但需要先对连续信号进行A/D变换,经数字滤波后,再经D/A转换得到所需要的连续信号(此处略)。 参考文献
1 倪养华,等．数字信号处理与实现．上海：上海交通大学出版社,1998 2 施阳．等．MATLAB语言工具箱．西安：西北工业大学出版
11
因篇幅问题不能全部显示，请点此查看更多更全内容