IIR滤波器的完全设计函数 - 肖伟

2、IIR滤波器的完全设计函数

以上介绍了IIR滤波器设计原理和基本方法步骤，并给出了一些例子说明如何用MATLAB编程分步实现这些步骤。由这些步骤可知我们必须多次调用MATLAB信号处理工具箱中的基本工具函数。

实际上，MATLAB信号处理工具箱还提供了IIR滤波器设计的完全工具函数，用户只要调用这些工具函数即可一次性完成设计，而不需要调用那些基本工具函数分步实现。 IIR滤波器设计的完全工具函数有butter,cheby1,cheby2,ellip。这些工具函数既可用于设计模拟滤波器，也可用于设计数字滤波器。这里介绍这些函数在IIR数字滤波器中的应用。在这两类滤波器设计中，这些工具函数调用格式基本相同，只是在频率处理上有所不同。 在MATLAB滤波器设计工具箱中，数字滤波器采用归一化频率，取值为0~1之间，归一化频率1对应的数字角频率为?，对应的真实频率为采样频率的一半。在应用MATLAB工具函数设计数字滤波器时应注意这一点。 数字IIR滤波器的完全设计函数有： [b,a]=butter(n,wn[,'ftype']) [z,p,k]=butter(n,wn[, 'ftype']) [b,a]=cheby1(n,Rp,wn[,'ftype']) [z,p,k]=cheby1(n, Rp,wn[,'ftype']) [b,a]=cheby2(n,Rs,wn[,'ftype']) [z,p,k]=cheby2(n, Rs,wn[,'ftype']) [b,a]=ellip(n,Rp,Rs,wn[,'ftype']) [z,p,k]=ellip(n, Rp,Rs,wn[,'ftype'])

在上面的调用方式中，n为滤波器的阶数，wn为滤波器的截止频率，取值为0~1。需根据采样频率Fs来定，如滤波器的截止频率为Fc（Hz），则wn的计算公式为：

wn?2*Fc （1） Fs这样就转换为0~1的归一化频率。其中wp,ws等边界频率都要根据此公式进行转换。 'ftype'滤波器的类型为：

‘high’为高通滤波器，截止频率为wn.

‘stop’为带阻滤波器，截止频率为wn=[w1,w2] (w

a,b分别为滤波器传递函数分子和分母多项式系数向量；z,p,k分别为滤波器的零极点和增益。Rp,Rs分别为所设计滤波器的通带波纹和阻带衰减，单位为dB。 设计好的数字滤波器传递函数具有下面形式：

B(z)b(1)?b(2)z?1???b(n?1)z?nH(z)?? (2) ?1?mA(z)a(1)?a(2)z???a(m?1)z上述函数采用双线性变换法和频率的预畸变处理将模拟滤波器离散化为数字滤波器，同时保证模拟滤波器和数字滤波器在wn(或w1,w2)处具有相同的幅频响应。

设计时应注意真实频率和MATLAB归一化数字频率之间的转换，即6-20式的应用。在进行IIR数字滤波器设计之前，请大家注意在模拟滤波器中我们用求取最小阶数和截止频率的函数如buttord,cheb1ord,cheb2ord,ellipord，这些函数完全可以用于IIR数字滤波器设计中，只不过在模拟滤波器设计中需加可选项's'，在数字滤波器中则不加该项。另外，输

出的截止频率也是归一化频率（归一化为0~1）。

【例4】设计一个Butterworth高通数字滤波器，通带边界频率为300Hz,阻带边界频率为200Hz，通带波纹小于1dB，阻带衰减大于20dB,采样频率为1000Hz。假设一个信号

x(t)?sin2?f1t?0.5cos2?f2t，其中f1=10Hz,f2=400Hz。试将原信号与通过该滤波器的输

出信号进行比较。

%Samp6_8

Fs=1000; %采样频率

wp=300*2/Fs;ws=200*2/Fs; %根据采样频率将滤波器边界频率进行转换 Rp=1;Rs=20; %通带波纹和阻带衰减 Nn=128; %显示滤波器频率特性的数据长度

[N,Wn]=buttord(wp,ws,Rp,Rs); %求得数字滤波器的最小阶数和截止频率（归一化频率）

[b,a]=butter(N,Wn,'high'); %设计Butterworth高通数字滤波器 figure(1)

[H,f]=freqz(b,a,Nn,Fs); %用Nn点绘出频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on; n=0:127;

dt=1/Fs;t=n*dt; %时间序列 f1=10;f2=400; %输入信号频率

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %输入信号 figure(2)

subplot(2,1,1),plot(t,x); title('输入信号') %绘制输入信号 y=filter(b,a,x); %对输入信号进行滤波

subplot(2,1,2),plot(t,y),title('输出信号') %绘制输出信号

xlabel('时间/s')

0-50振幅/dB-100-150-200050100150200250300频率/Hz3504004505001000相位/o-100-200-300-400050100150200250300频率/Hz350400450500

图 所设计的高通Butterworth滤波器的频率特性

上图：幅频特性；下图：相频特性

输入信号210-1-200.020.040.060.080.10.120.14输出信号10.50-0.500.020.040.060.08时间/s0.10.120.14

图 所设计滤波器的输入和输出信号

由上图可以看到所设计滤波器在大于300Hz为通带，其衰减均小于1dB；小于200Hz为阻带，其衰减大于20dB，完全符合设计要求，但相频特性是非线性的。由图可以看出，当滤波器输入10Hz和400Hz两种信号后，滤波器滤除了10Hz的信号，使得400Hz的信号通过了滤波器，起到了滤波的效果。 【例5】设计一个带通Chebyshev I型数字滤波器，通带为100Hz~200Hz,过渡带宽均为50Hz,通带波纹小于1dB，阻带衰减大于30dB，采样频率Fs=1000Hz。假设一个信号

x(t)?sin2?f1t?0.3cos2?f2t?0.1sin2?f3t，其中f1=30Hz,f2=100Hz,f3=270Hz。试将原

信号与经过该滤波器滤波后输出信号进行比较。

该滤波器的通带范围为100~200Hz，过渡带宽为50Hz，因此，通带边界频率为100和200Hz，阻带边界频率为50和250Hz。

%Samp6_9

Fs=1000; %采样频率

wp=[100 200]*2/Fs; %通带边界频率（归一化频率）（6-20式） ws=[50 250]*2/Fs; %阻带边界频率（归一化频率）（6-20式）

Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数

[N,Wn]=cheb1ord(wp,ws,Rp,Rs);%求得数字滤波器的最小阶数和归一化截止频率 [b,a]=cheby1(N,Rp,Wn); %按最小阶数、通带波纹和截止频率设计数字滤波器 figure(1)

[H,f]=freqz(b,a,Nn,Fs); %求得滤波器的频率特性

subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))

xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2)

f1=30;f2=100;f3=270; %输入信号的三种频率成分 N=100; %输入信号的数据点数

dt=1/Fs;n=0:N-1;t=n*dt; %时间序列

x=sin(2*pi*f1*t)+0.3*cos(2*pi*f2*t)+0.1*sin(2*pi*f3*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号进行滤波 subplot(2,1,2),plot(t,y) %绘制输出信号 ylim([-0.2 0.3])

title('输出信号'),xlabel('时间/s')

0-100振幅/dB-200-300-400050100150200250300频率/Hz350400450500400200相位/o0-200-400050100150200250300频率/Hz350400450500

图 所设计滤波器的频率特性 上图：幅频特性；下图：相频特性

输入信号210-1-200.010.020.030.040.050.060.070.080.090.1输出信号0.30.20.10-0.1-0.200.010.020.030.040.050.06时间/s0.070.080.090.1

图 所设计滤波器的输入和输出信号

由上图可看出，在100~200Hz为通带，其衰减小于1dB；在50Hz以下和250Hz以上，其衰减大于30dB，完全符合滤波器的设计要求。下图表明滤波器的相位为非线性。将30Hz、100Hz和270Hz的合成振动输入滤波器后，可以看到完全滤除了在阻带范围内的30Hz和270Hz的振动，起到了滤波效果。

【例6】设计一个带通Chebyshev II型滤波器，设计参数与例5相同。假设一个信号

x(t)?sin2?f1t?0.8cos2?f2t，其中f1=10Hz,f2=100Hz。试将原信号与通过该滤波器的输

出信号进行比较。采样频率为1000Hz。

%Samp6_10

Fs=1000; %采样频率

wp=[100 200]*2/Fs;ws=[50 250]*2/Fs; %通带和阻带边界频率（归一化频率）（6-20式） Rp=1;Rs=30;Nn=128; %通带波纹和阻带衰减以及绘制频率特性的数据点数 [N,Wn]=cheb2ord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率 [b,a]=cheby2(N,Rs,Wn); %按最小阶数截止频率和阻带衰减设计数字滤波器 figure(1)

[H,f]=freqz(b,a,Nn,Fs); %按传递函数系数、数据点数和采样频率求频率特性 subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on; subplot(2,1,2),plot(f,180/pi*unwrap(angle(H))) xlabel('频率/Hz');ylabel('相位/^o');grid on;

figure(2)

f1=10;f2=100; %输入信号的频率成分 N=500; %输入信号的数据点数 dt=1/Fs;n=0:N-1;t=n*dt; %时间序列

x=sin(2*pi*f1*t)+0.8*cos(2*pi*f2*t); %输入信号

subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号进行滤波，输出为y. subplot(2,1,2),plot(t,y) %绘制输出信号波形 title('输出信号'),xlabel('时间/s')

0-20振幅/dB-40-60-80050100150200250300频率/Hz350400450500400200相位/o0-200-400050100150200250300频率/Hz350400450500

图 所设计滤波器的频率特性 上图：幅频特性；下图：相频特性

输入信号210-1-200.050.10.150.20.25输出信号0.30.350.40.450.510.50-0.5-100.050.10.150.20.25时间/s0.30.350.40.450.5

图 所设计滤波器的输入和输出信号

可以看出，100~200Hz之间为通带，其衰减不大于1dB；在50Hz以下和250Hz以上，衰减大于30dB，完全符合滤波器的设计要求。但在阻带内有振荡，这是Chebyshev II型滤波器的特点。当滤波器输入10Hz和100Hz两种信号的合成时，滤波器可以滤除处于阻带内的10Hz的低频信号，达到滤波的目的。

【例7】设计一个带阻椭圆滤波器，阻带频率从100Hz~200Hz,通带波纹小于1dB,阻带衰减为50dB,两边过渡带宽为50Hz,采样频率1000Hz。假设一个信号

x(t)?sin2?f1t?0.5cos2?f2t?0.5sin2?f3t，其中f1=50Hz,f2=200Hz,f3=270Hz。试将原

信号与通过该滤波器的输出信号进行比较。

%Samp6_11

Fs=1000; %采样频率

ws=[100 200]*2/Fs;wp=[50 250]*2/Fs; %通带和阻带边界频率（归一化频率） Rp=1;Rs=50;Nn=512; %通带波纹和阻带衰减以及绘制频率特性的数据点数

[N,Wn]=ellipord(wp,ws,Rp,Rs); %求取数字滤波器的最小阶数和归一化截止频率 [b,a]=ellip(N,Rp,Rs,Wn,'stop'); %按最小阶数、截止频率、通带波纹和阻带衰减设计滤波器 figure(1)

[H,f]=freqz(b,a,Nn,Fs); %按传递函数系数、数据点数和采样频率求得滤波器的频率特性

subplot(2,1,1),plot(f,20*log10(abs(H))); xlabel('频率/Hz');ylabel('振幅/dB');grid on;

subplot(2,1,2),plot(f,180/pi*unwrap(angle(H)))

xlabel('频率/Hz');ylabel('相位/^o');grid on; figure(2)

f1=50;f2=200;f3=270; %输入信号的频率成分 N=100; %输入信号的数据点数 dt=1/Fs;n=0:N-1;t=n*dt;%时间序列

x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t)+0.5*sin(2*pi*f3*t); %输入信号 subplot(2,1,1),plot(t,x),title('输入信号') %绘制输入信号 y=filtfilt(b,a,x); %对输入信号滤波 subplot(2,1,2),plot(t,y) %绘制输出信号 title('输出信号'),xlabel('时间/s')

0-20振幅/dB-40-60-80-100050100150200250300频率/Hz350400450500500相位/o0-500050100150200250300频率/Hz350400450500

图 所设计滤波器的频率特性 上图：幅频特性；下图：相频特性

输入信号210-1-200.010.020.030.040.050.060.070.080.090.1输出信号1.510.50-0.5-100.010.020.030.040.050.06时间/s0.070.080.090.1

图 所设计滤波器的输入和输出信号

可以看到，在100~200Hz为阻带，衰减均大于50dB，在其他频率处为通带，衰减小于1dB，完全符合滤波器的设计要求。但在通带和阻带内有振荡，这是椭圆滤波器的特点。当滤波器输入50Hz、200Hz和270Hz的合成信号时，滤波器可以滤除处于阻带内的200Hz信号，只剩下50Hz和270Hz的合成信号，达到滤波的目的。

本文来源：https://www.bwwdw.com/article/a9fa.html