数字滤波器的设计共5页文档

模拟滤波器到数字滤波器的转换
一、脉冲响应不变法设计IIR数字滤波器
impinvar
功能：用脉冲响应不变法实现模拟到数字的滤波器变换。

调用格式：
[bd,ad]=impinvar(b,a,Fs);将模拟滤波器系数b,a变换成数字的滤波器系数bd,ad,两者的冲激响应不变。

[bd,ad]=impinvar(b,a)；采用Fs的缺省值1Hz.
例：采用脉冲响应不变法设计一个切比雪夫I型数字带通滤波器，要求：通带w p1=0.3pi, W p2=0.7pi, R p=1dB, 阻带w s1=0.1pi, W s2=0.9pi, A s=15dB, 滤波器采样频率为
F s=2000Hz.
Matlab程序：
%数字滤波器指标
w p1=0.3*pi; w p2=0.7*pi;
w s1=0.1*pi; w s2=0.9*pi;
R p=1; A s=15;
%转换为模拟滤波器指标
Fs=2000; T=1/Fs;
Omgp1=wp1*Fs; Omgp2=wp2*Fs; %模拟滤波器的通带截止频率
Omgp=[Omgp1,Omgp2];
Omgs1=ws1*Fs; Omgs2=ws2*Fs; %模拟滤波器的阻带截止频率
Omgs=[Omgs1,Omgs2];
Bw=Omgp2-Omgp1; w0=sqrt(Omgp1*Omgp2); %模拟通带带宽和中心频率
%模拟原型滤波器计算
[n,omgn]=cheb1ord(omgp,Omgs,Rp,As,’s’);
[z0,p0,k0]=cheb1ap(n,Rp); %设计归一化的模拟原型滤波器（zpk模型）
ba1=k0*real(poly(z0)); %求原型滤波器系统函数分子系数b
aa1=real(poly(p0)); %求原型滤波器系统函数分母系数a
[ba,aa]=lp2bp(ba1,aa1,w0,bw); %变换为模拟带通滤波器
%用脉冲响应不变法计算数字滤波器系数
[bd,ad]=impinvar(ba,aa,Fs);
%求数字系统的频率特性
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H))); %将幅度化为分贝值
%作图
subplot(2,2,3),plot(w/pi,dbH);
axis([0,1,-50,1]); title('实际带通相对幅度');
ylabel('dB');xlabel('数字频率（w/pi）');
set(gca,'Xtick',[0,wp1/pi,ws1/pi,wp2/pi,ws2/pi,1]);
set(gca,'Ytick',[-50,-20,-3,-1]); grid
subplot(2,2,4),plot(w/pi, angle(H)/pi*180);
axis([0,1,-200,200]);title('实际数字带通相位');
set(gca,'Xtick',[ 0,wp1/pi,ws1/pi,wp2/pi,ws2/pi,1]);
set(gca,'Ytick',[-180,-120,0,90,180]);grid
ylabel('\phi');xlabel('数字频率（w/pi）');
二、用双线性变换法设计IIR数字滤波器
bilinear
功能：双线性变换——将s域（模拟域）映射到z域（数字域）的标准方法，将模拟滤波器变换成离散等效滤波器。

调用格式:
1、[numd,dend]=bilinear(num,den,Fs)；将模拟域传递函数（tf模型）变换为数字域传递函数，Fs为取样频率。

2、[numd,dend]=bilinear(num,den,Fs,Fp)；将模拟域传递函数变换为数字域传递函数，Fs 为取样频率。

Fp为通带截止频率。

3、[zd,pd,kd]=bilinear(z,p,k,Fs)；将模拟域传递函数（zpk模型）变换为数字域传递函数，Fs为取样频率。

4、[zd,pd,kd]=bilinear(z,p,k,Fs,Fp)；将模拟域传递函数变换为数字域传递函数，Fs为取样频率。

Fp为通带截止频率。

例：采样双线性变换设计一个巴特沃斯数字低通滤波器，要求：
wp=0.25pi,Rp=1dB,ws=0.4pi,As=15dB.滤波器采样频率Fs=100Hz。

Matlab程序：
%数字滤波器指标
wp=0.25*pi;Rp=1;
ws=0.4*pi;As=15;
%转换为模拟滤波器指标
Fs=100; T=1/Fs;
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
%模拟原型滤波器设计
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,’s’);
[z0,p0,k0]=buttap(n); %归一化原型设计（zpk模型）
ba=k0*real(poly(z0)); %原型系统函数分子系数b
aa=real(poly(p0)); %原型系统函数分母系数a
[ba1,aa1]=lp2lp(ba,aa,Omgc);
%以上4行求滤波器系数ba1,aa1的程序，可以由以下一条程序代替
[ba1,aa1]=butter(n,Omgc,’s’);
%用双线性变换计算数字滤波器系数
[bd,ad]=bilinear(ba1,aa1,Fs);
%求数字系统的频率特性
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
%作图
三、IIR数字滤波器的直接设计
1、butter
功能：巴特沃斯（butterworth）模拟或数字滤波器设计
调用格式：
1）、[b,a]=butter(n,wn),设计截止频率为wn的n阶巴特沃斯数字滤波器（tf模型），
其中，截止频率是幅度下降到3分贝处的频率，wn在0到1之间，1对应0.5Fs(取样频率)。

wn=[w1,w2]时，产生数字带通滤波器。

2）、[b,a]= butter(n,wn,’ftype’),当ftype=high时，设计高通滤波器，当ftype=stop时，设计带阻滤波器，此时wn=[w1,w2]。

3）、[b,a]=butter(n,wn,’s’), 设计截止频率为wn的n阶巴特沃斯模拟低通或带通滤波器
（tf模型），
4)、2）、[b,a]= butter(n,wn,’ftype’,’s’),当ftype=high时，设计模拟高通滤波器，当ftype=stop 时，设计模拟带阻滤波器，此时wn=[w1,w2s]。

2、cheby1
功能：切比雪夫I滤波器设计（通带等波纹）
调用格式：
1）、[b,a]=cheby1(n,Rp,Wn);设计截止频率为wn的n阶切比雪夫I型数字低通和带通滤波器。

2）、[b,a]= cheby1(n, Rp,Wn,’ftype’), 设计截止频率为wn的n阶切比雪夫I型数字高通和带阻滤波器。

3）、1）、[b,a]=cheby1(n,Rp,Wn,’s’);设计截止频率为wn的n阶切比雪夫I型数模拟低通和带通滤波器。

4）、[b,a]= cheby1(n, Rp,Wn,’ftype’,’s’), 设计截止频率为wn的n阶切比雪夫I型模拟高通和带阻滤波器。

3、cheby2
功能：切比雪夫II滤波器设计（阻带等波纹）
调用格式：
1）、[b,a]=cheby2(n,As,Wn);设计截止频率为wn的n阶切比雪夫II型数字低通和带通滤波器。

2）、[b,a]= cheby1(n, As,Wn,’ftype’), 设计截止频率为wn的n阶切比雪夫II型数字高通和带阻滤波器。

3）、1）、[b,a]=cheby1(n,As,Wn,’s’);设计截止频率为wn的n阶切比雪夫II型数模拟低通和带通滤波器。

4）、[b,a]= cheby1(n, As,Wn,’ftype’,’s’), 设计截止频率为wn的n阶切比雪夫II型模拟高通和带阻滤波器。

4、ellip
功能：椭圆滤波器设计
调用格式：
1）、[b,a]=ellip(n,Rp,As,Wn);设计截止频率为wn的n阶椭圆数字低通和带通滤波器。

2）、[b,a]= ellip(n,Rp As,Wn,’ftype’), 设计截止频率为wn的n阶椭圆型数字高通和带阻滤波器。

3）、1）、[b,a]= ellip(n,Rp As,Wn,’s’);设计截止频率为wn的n阶椭圆型数模拟低通和带通滤波器。

4）、[b,a]= ellip(n,Rp As,Wn,’ftype’,’s’), 设计截止频率为wn的n阶椭圆型模拟高通和带阻滤波器。

（重点）例：用直接法设计一个巴特沃斯高通数字滤波器，要求：wp=0.25pi, Rp=1dB, ws=0.4pi, As=20dB,采样频率Fs=200Hz，要求描绘其幅频特性和相频特性曲线。

MATLAB程序：
ws=0.25; wp=0.4;
Rp=1; As=20;
Fs=200;
[n,wc]=buttord(wp,ws,Rp,As);
[b,a]=butter(n,wc,’high’);
[h,w]=freqz(b,a);
plot(w/pi,abs(h));
（重点）例：用直接法设计一个切比雪夫I型数字带通滤波器，要求：
wp1=0.25pi,wp2=0.75pi,Rp=1dB,ws1=0.85pi,ws2=0.15pi,As=20dB,描绘滤波器归一化的绝对
和相对幅频特性，相频特性，零极点分布图，列出系统传达函数式。

MATLAB程序：
ws1=0.15; ws2=0.85; ws=[ws1,ws2];
wp1=0.25; wp2=0.75; wp=[wp1,wp2];
Rp=1; As=20;
[n,wc]=cheb1ord(wp,ws,Rp,As);
[b,a]=cheby1(n,Rp,wc);
[H,w]=freqz(b,a);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1), plot(w/pi,abs(H));
subplot(2,2,2), plot(w/pi,angle(H)) ;
subplot(2,2,3), plot(w/pi,dbH);
subplot(2,2,4), zplane(b,a);
四、线性相位FIR数字滤波器的设计
理想的数字低通滤波器频率响应函数：ideal_lp.m
Function hd=ideal_lp(wc,N)
%hd=点0到Ｎ－１之间的理想脉冲响应
％ｗｃ＝截止频率（弧度）
％Ｎ＝理想滤波器的长度
tao=(N-1)/2;
n=[0(N-1)];
m=n-tao+eps;
hd=sin(wc*m)./(pi*m);
其他选频滤波器则可以由低通频响特性合成，例如一个带通在wc1~wc2之间的带通滤波器，在给定N值的条件下，可以用下列程序实现：
Hd=ideal_lp(wc2,N)-ideal_lp(wc1,N);
例：用矩形窗设计一个FIR数字低通滤波器，要求：N=64，截止频率为wc=0.4pi，描绘理想和实际滤波器的脉冲响应，窗函数及滤波器的幅频响应曲线。

MATLAB程序；
wc=0.4*pi;
N=64;n=0:N-1;
hd=ideal_lp(wc,N); %建立理想低通滤波器
windows=(boxcar(N))’; %使用矩形窗，并将列向量变为行向量
b=hd.*windows; %求FIR系统函数系数
[H,w]=freqz(b,1); %求解频率特性
dbH=20*log10((abs(H)+eps)/max(abs(H)));
%作图
例：选择合适的窗函数设计一个FIR数字低通滤波器，要求：通带截止频率为
wp=0.3pi,Rp=0.05dB;阻带截止频率为ws=0.45pi, As=50dB。

描绘该滤波器的脉冲响应，窗函数及滤波器的幅频响应曲线和相频响应曲线。

wp=0.3*pi;ws=0.45pi;
deltaw=ws-wp; %计数过度带的带宽
N0=ceil(6.6*pi/deltaw); %按表所示选哈明窗，求滤波器的长度N0
N=N0+mod(N0+1,2); %为实现FIR类型I偶对称滤波器，应确保N为奇数windows=(hamming(N))＇；%使用哈明窗
wc=(ws+wp)/2;
hd=ideal_lp(wc,N);
b=hd.*windows;
[db,mag,pha,grd,w]=freqz_m(b,1);
n=0:N-1; dw=2*pi/1000; %dw为频率分辨率，将0~2pi分为1000份Rp=-(min(db(1:wp/dw+1))); %检验通带波动
As=-round(max(db(ws/dw+1:501))); %检验最小阻带衰减。