实验五：双线性变换法设计IIR 数字滤波器

一、实验目的：

1、熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。

2、掌握数字滤波器的计算机仿真方法。

3、熟悉Batterworth滤波器设计方法及特点

二、实验原理

（一）、IIR数字滤波器的设计步骤：

① 按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标；

② 根据模拟滤波器技术指标设计为响应的模拟低通滤波器；

③ 跟据脉冲响应不变法和双线性不变法把模拟滤波器转换为数字滤波器；

④ 如果要设计的滤波器是高通、带通或带阻滤波器，则首先把它们的技术指标转化为模拟低通滤波器的技术指标，设计为数字低通滤波器，最后通过频率转换的方法来得到所要的滤波器。

在MATLAB中，经典法设计数字滤波器主要采用以下步骤：

IIR数字滤波器设计步骤

（二）、用模拟滤波器设计数字滤波器的方法

1、冲激响应不变法：

冲激响应不变法是从时域出发，要求数字滤波器的冲激响应h(n) 对应于模拟滤波器h(t) 的等间隔抽样。

优点：时域逼近良好；保持线性关系。

缺点：频域响应混叠。只适用于限带低通滤波器和带通滤波器

2、双线性变换法

优点：克服了频域混叠

缺点：高频时会引起畸变

1）冲激响应不变法impinvar

格式：[BZ,AZ]= impinvar（B,A,Fs）

功能：把具有[B,A]模拟滤波器传递函数模型转换为采样频率为Fs的数字滤波器的传递函数模型[BZ,AZ]，Fs默认值为1。

例：一个4阶的Butterworth模拟低通滤波器的系统函数如下：

试用冲激响应不变法求出Butterworth模拟低通数字滤波器的系统函数。

num=1;

den=[1,sqrt(5),2,sqrt(2),1];

[num1,den1]=impinvar(num,den)

2）双线性变换法bilinear

格式一：[Zd,Pd,Kd]= bilinear(Z,P,K,Fs)

功能：把模拟滤波器的零极点模型转换成数字滤波器的零极点模型，Fs是采样频率

格式二：[numd,dend]= bilinear(num,den,Fs)

功能：把模拟滤波器的传递函数模型转换为数字滤波器的传递函数模型。

例：一个三阶的模拟Butterworth模拟低通滤波器的系统函数如下：

，试用双线性变换法求出数字Butterworth数字低通滤波器的系统函数。

num=1;

den=[1,sqrt(3),sqrt(2),1];

[num1,den1]=bilinear(num,den,1)

3）IIR数字滤波器的频率变换实现

步骤：

1按一定的规则将数字滤波器的技术指标转换为模拟低通滤波器的技术指标

2根据转换后的技术指标使用滤波器阶数函数，确定滤波器的最小阶数N和截止频率Wc

3利用最小阶数N产生模拟低通原型

4利用截止频率Wc把模拟低通滤波器原型转换为模拟低通、高通、带通、带阻滤波器

5利用冲激响应不变法或双线性变换法把模拟滤波器转换为数字滤波器

表一 IIR滤波器阶次估计

1．数字滤波器的设计参数

滤波器的4个重要的通带、阻带参数为：

：通带截止频率（Hz）      ：阻带起始频率（Hz）

：通带内波动（dB），即通带内所允许的最大衰减；

：阻带内最小衰减

设采样速率（即奈奎斯特速率）为，将上述参数中的频率参数转化为归一化角频率参数：

：通带截止角频率（rad/s） ，；     

：阻带起始角频率（rad/s） ， 

通过以上参数就可以进行离散滤波器的设计。

2、巴特沃斯滤波器设计

1）巴特沃斯滤波器阶数的选择：

在已知设计参数，，，之后，可利用“buttord”命令可求出所需要的滤波器的阶数和3dB截止频率，其格式为：

[n，Wn]=buttord[Wp，Ws，Rp，Rs]，其中Wp，Ws，Rp，Rs分别为通带截止频率、阻带起始频率、通带内波动、阻带内最小衰减。返回值n为滤波器的最低阶数，Wn为3dB截止频率。

2）巴特沃斯滤波器系数计算：

由巴特沃斯滤波器的阶数n以及3dB截止频率Wn可以计算出对应传递函数H(z)的分

子分母系数，MATLAB提供的命令如下：

(a)巴特沃斯低通滤波器系数计算：

[b，a]=butter(n,Wn)，其中b为H(z)的分子多项式系数，a为H(z)的分母多项式系数

(b)巴特沃斯高通滤波器系数计算：[b，a]=butter(n，Wn，’High’) 

(c)巴特沃斯带通滤波器系数计算：[b，a]=butter(n，[W1，W2])，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的是2*n阶滤波器系数。

(d)巴特沃斯带阻滤波器系数计算：[b，a]=butter(ceil(n/2)，[W1，W2]，’stop’)，其中[W1，W2]为截止频率，是2元向量，需要注意的是该函数返回的也是2*n阶滤波器系数。

三、巴特沃斯滤波器设计实例：

例题1：采样速率为8000Hz，要求设计一个低通滤波器， =2100Hz， =2500Hz， =3dB， =25dB。

1、采样速率为10000Hz，要求设计一个巴特沃斯带阻滤波器， =[1000Hz，1500Hz]， =[1200Hz，1300Hz]， =3dB， =30dB。

用直接设计法

程序如下：

fn=8000;%采样频率 

fp=2100;  %通带截止频率

fs=2500; %阻带起始频率

Rp=3; %通带最大衰减

Rs=25;%阻带最小衰减

Wp=fp/(fn/2);%计算归一化角频率

Ws=fs/(fn/2);

[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率

[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数

[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)

subplot(2,1,1)

plot(F,20*log10(abs(H))) %画出幅频特性图

xlabel('Frequency(Hz)'); ylabel('Magnitude(dB)') 

title('低通滤波器')

axis([0 4000 -30 3]);grid on

pha=angle(H)*180/pi;

subplot(2,1,2)

plot(F,pha);grid on %画出相频特性图

xlabel('Frequency(Hz)'); 

ylabel('phase');

用双线性变换法

wp=2100*2*pi; %利用 

ws=2500*2*pi;

Rp=3;

Rs=25;

Fs=8000;

Ts=1/Fs;

%选择滤波器的最小阶数

[N,Wn]=buttord(wp,ws,Rp,Rs,'s');%创建butterworth模拟滤波器

[Z,P,K]=buttap(N);

%把滤波器零极点模型转化为传递函数模型

[Bap,Aap]=zp2tf(Z,P,K); 

%把模拟滤波器原型转换成截至频率为Wn的低通滤波器

[b,a]=lp2lp(Bap,Aap,Wn);

%用双线性变换法实现模拟滤波器到数字滤波器的转换

[bz,az]=bilinear(b,a,Fs);

%绘制频率响应曲线

[H,W]=freqz(bz,az);

plot(W*Fs/(2*pi),abs(H));

grid

xlabel('频率/Hz')

ylabel('幅度')

例题2：模拟原型直接变换法设计数字滤波器：

已知四阶归一化低通巴特沃斯模拟滤波器系统函数为，编写MATLAB程序实现从设计3dB截止频率为，设采样周期为T=1，的四阶低通巴特沃斯数字滤波器。

程序如下：步骤一：将设计内容题所给归一化巴特沃斯低通滤波器以3dB截止频率为进行去归一化。

步骤二：用双线性变化法将低通模拟滤波器变换为低通数字滤波器

3、已知四阶归一化低通巴特沃斯模拟滤波器系统函数为，编写MATLAB程序实现从设计3dB截止频率为的四阶高通巴特沃斯数字滤波器。

设计程序如下：

clear;       

T=1;    fs=1/T;   N=4;

wc=pi/2;  

omegach=2*tan(wc/2)/T;%模拟滤波器的截止频率

M=1;  

 N=[1,2.6131,3.4142,2.6131,1];

[h,w]=freqs(M,N,512);     %模拟滤波器的幅频响应

subplot(2,1,1);

plot(w,20*log10(abs(h))); 

axis([0,10,-90,0]),grid on;    

xlabel('Hz');ylabel('幅度'); title('归一化模拟低通滤波器');

[Ms,Ns]=lp2lp(M,N,omegach);           %对低通滤波器进行频率变换

[hs,ws]=freqs(Ms,Ns,512);               %模拟滤波器的幅频响应

subplot(2,1,2);plot(ws,20*log10(abs(hs))); 

grid;

axis([0,10,-90,0]);

xlabel('Hz');ylabel('幅度'); title('去归一化模拟低通滤波器');

[Mz,Nz]=bilinear(Ms,Ns,1/T);            %对模拟滤波器双线性变换

[h1,w1]=freqz(Mz,Nz);                  %数字滤波器的幅频响应

figure

plot(w1/pi,20*log10(abs(h1)));   grid;

xlabel('ω/π');ylabel('幅度(dB)'); title('数字低通滤波器');

axis([0,1,-160,0])

四、实验内容

1、采样速率为10000Hz，要求设计一个巴特沃斯带阻滤波器， =[1000Hz，1500Hz]， =[1200Hz，1300Hz]， =3dB， =30dB。

程序：

fn=10000;%采样频率 

fp=[1000,1500];  %通带截止频率

fs=[1200,1300]; %阻带起始频率

Rp=3; %通带最大衰减

Rs=30;%阻带最小衰减

Wp=fp/(fn/2);%计算归一化角频率

Ws=fs/(fn/2);

[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率

[b,a]=butter(n,Wn,'stop');%计算H(z)分子、分母多项式系数

[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)

subplot(2,1,1)

plot(F,20*log10(abs(H))) %画出幅频特性图

xlabel('Frequency(Hz)'); ylabel('Magnitude(dB)') 

title('带阻滤波器')

axis([0 4000 -30 3]);grid on

pha=angle(H)*180/pi;

subplot(2,1,2)

plot(F,pha);grid on %画出相频特性图

xlabel('Frequency(Hz)'); 

ylabel('phase');

提示：[b,a]=butter(N,Wc,'stop')

2、 采样速率为10000Hz，要求设计一个带通滤波器， =[1000Hz，1500Hz]， =[600Hz，1900Hz]， =3dB， =20dB。

程序：

fn=10000;%采样频率 

fp=[1000,1500];  %通带截止频率

fs=[600,1900]; %阻带起始频率

Rp=3; %通带最大衰减

Rs=20;%阻带最小衰减

Wp=fp/(fn/2);%计算归一化角频率

Ws=fs/(fn/2);

[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率

[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数

[H,F]=freqz(b,a,1000,8000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)

subplot(2,1,1)

plot(F,20*log10(abs(H))) %画出幅频特性图

xlabel('Frequency(Hz)'); ylabel('Magnitude(dB)') 

title('带通滤波器')

axis([0 4000 -30 3]);grid on

pha=angle(H)*180/pi;

subplot(2,1,2)

plot(F,pha);grid on %画出相频特性图

xlabel('Frequency(Hz)'); 

ylabel('phase')；

3、已知四阶归一化低通巴特沃斯模拟滤波器系统函数为，编写MATLAB程序实现从设计3dB截止频率为的四阶高通巴特沃斯数字滤波器。

clear;       

T=1;    fs=1/T;   N=4;

wc=pi/4;  

omegach=2*tan(wc/2)/T;%模拟滤波器的截止频率

M=1;  

 N=[1,sqrt(5),2,sqrt(2),1];

[h,w]=freqs(M,N,512);     %模拟滤波器的幅频响应

subplot(2,1,1);

plot(w,20*log10(abs(h))); 

axis([0,10,-90,0]),grid on;    

xlabel('Hz');ylabel('幅度'); title('归一化模拟低通滤波器');

[Ms,Ns]=lp2lp(M,N,omegach);           %对低通滤波器进行频率变换

[hs,ws]=freqs(Ms,Ns,512);               %模拟滤波器的幅频响应

subplot(2,1,2);plot(ws,20*log10(abs(hs))); 

grid;

axis([0,10,-90,0]);

xlabel('Hz');ylabel('幅度'); title('去归一化模拟低通滤波器');

[Mz,Nz]=bilinear(Ms,Ns,1/T);            %对模拟滤波器双线性变换

[h1,w1]=freqz(Mz,Nz);                  %数字滤波器的幅频响应

figure

plot(w1/pi,20*log10(abs(h1)));   grid;

xlabel('ω/π');ylabel('幅度(dB)'); title('数字低通滤波器');

axis([0,1,-160,0])

4、设计低通滤波器，把输入信号中的三个信号分离出来。

要求：画出滤波前后信号的波形及频谱及低通滤波器的幅频响应。

提示：[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数

y=filter(b,a,x);%对输入的信号进行滤波

程序：

fn=2000;%采样频率

N=2000;%数据点数

n=0:N-1;

t=0:1/fs:200/fs;%采样时间序列

f0=100;%信号频率

x=cos(2*pi*f0*t)+cos(2*pi*200*t)+cos(2*pi*400*t);

subplot(3,1,1);

plot(t,x);

xlabel('t');

ylabel('sin(2*pi*100*t)');

title('时域信号');

Y=fft(x,N);%对信号进行FFT变换

magY=abs(Y);%求得FFT变换后的幅度

f=n*fs/N;%频率序列

subplot(3,1,2);

plot(f(1:N/2),magY(1:N/2));%画出幅频响应

xlabel('f');

ylabel('幅度');

title('N=2000');

grid;

%滤波器

fp=90;  %通带截止频率

fs=110; %阻带起始频率

Rp=3; %通带最大衰减

Rs=25;%阻带最小衰减

Wp=fp/(fn/2);%计算归一化角频率

Ws=fs/(fn/2);

[n,Wn]=buttord(Wp,Ws,Rp,Rs);%计算阶数和截止频率

[b,a]=butter(n,Wn);%计算H(z)分子、分母多项式系数

[H,F]=freqz(b,a,1000,2000);%计算H(z)的幅频响应,freqz(b,a,计算点数,采样速率)

subplot(3,1,3)

plot(F,20*log10(abs(H))) %画出幅频特性图

axis([0,1000,-200,3]);

y=filter(b,a,x);

figure;

plot(t,y);