巴特沃斯带阻IIR数字滤波器设计

2.巴特沃斯带阻IIR数字滤波器设计

1.设计思路—基于冲激响应不变法的IIR数字滤波器设计

冲激响应不变法的设计原理是利用数字滤波器的单位抽样响应序列H(z)来逼近模拟滤波器的冲激响应g(t)。

按照冲激响应不变法的原理，通过模拟滤波器的系统传递函数G(s)，可以直接求得数字滤波器的系统函数H(z)，其转换步骤如:

(1)利用ω=ΩT(可由关系式Z=e sT推导出)，将ωp, ωs转换成Ωp, Ωs ，而αp，αs不变；

(2)

求解低通模拟滤波器的传递函数G(s)；

(3)将模拟滤波器的传递函数G(s)转换为数字滤波器的传递函数H(z)。

尽管通过冲激响应不变法求取数字滤波器的系统传递函数比较方便，并具有良好的时域逼特性，但若G(s)不是带限的，或是抽样频率不高，那么在H(e jω)中将发生混叠失真，数字滤波器的频率响应不能重现模拟滤波器的频率响应。只有当模拟滤波器的频率响应在超过折叠频率后的衰减很大时，混叠失真才很小，此时采样脉冲响应不变法设计的数字滤波器才能满足设计的要求，这是冲激响应不变法的一个严重的缺点。

2.设计要求及方案

设计一带阻巴斯沃特IIR滤波器，要求如下：

带纹波为Rp＝1dB，

通带上、下限角频率为0.11π、0.81π，

阻带上、下限角频率为0.31π、0.61π，

阻带最小衰减αs=40dB，采样频率f s=15000Hz

3．用MTALAB算法设计巴特沃斯带阻IIR数字滤波器

fs=15000;T=1/fs;

rp=1;rs=40;

wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi; %数字带阻滤波器技术指标

wc1=(2/T)*tan(wp1/2); %频率预畸变

wc2=(2/T)*tan(wp2/2);wr1=(2/T)*tan(ws1/2);wr2=(2/T)*tan(ws2/2);

w0=sqrt(wc1*wc2);B=wc2-wc1;

wp=1; %归一化通带截止频率

ws=wp*(wr1*B)/(w0^2-wr1^2); %归一化阻带截止频率

[N,wc]=buttord(wp,ws,rp,rs,'s'); %求滤波器阶数和3dB截止频率

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

[Md,Nd]=zp2tf(Z,P,K); %将零极点形式转换为传输形式[M,N]=lp2bs(Md,Nd,w0,B); %对低通滤波器进行频率变换，转换为带阻滤波器[h,w]=freqs(M,N); %模拟带阻滤波器的幅频响应

plot(w/(2*pi,abs(h)));grid; >> xlabel('频率/Hz');ylabel('幅度');title('模拟带阻滤波器');

[b,a]=bilinear(M,N,15000); %对模拟滤波器双线性变换

figure(1);

freqz(b,a);[H,W]=freqz(b,a); %绘出频率响应

axis([0,1,-100,20]);

figure(2);

plot(W*fs/(2*pi),abs(H));grid on;

xlabel('频率/Hz');ylabel('幅值');

title('数字滤波器幅频响应|H(ejOmega)| ');

仿真出的幅频响应曲线图如下图2.1所示：

图2.1：幅频响应曲线

相频特性及幅度特性曲线如下图2.2所示：

图2.2：相频特性与幅度特性曲线

fs=15000;T= 1/fs; rp=1;rs=40;

wp1=0.11*pi;wp2=0.81*pi;ws1=0.31*pi;ws2=0.61*pi;%数字带阻滤波器技术指标wc1=(2/T)*tan(wp1/2);%频率预畸变

wc2=(2/T)*tan(wp2/2);

wr1=(2/T)*tan(ws1/2);

wr2=(2/T)*tan(ws2/2);

w0=sqrt(wc1*wc2);

B=wc2-wc1;

wp=1;%归一化通带截止频率

ws=wp*(wr1*B) / (w0^2-wr1^2) %归一化阻带截止频率

[N,wc]=buttord(wp,ws,rp,rs,'s')%求滤波器阶数和3dB截止频率

[Z,P,K]=buttap(N)%设计模拟低通滤波器

[Md,Nd]=zp2tf(Z,P,K)%将零极点形式转换为传输函数形式

[M,N]=lp2bs(Md,Nd,w0,B)%对低通滤波器进行频率变换，转换为带阻滤波器[h,w]=freqs(M,N);%模拟带阻滤波器的幅频响应

plot(w/(2*pi),abs(h));grid;

xlabel('频率Hz');ylabel('幅度');title('模拟带阻滤波器');

[b,a]=bilinear(M,N,15000)%对模拟滤波器双线性变换

figure(1);

reqz(b,a);[H,W]=freqz(b,a); %绘出频率响应；

axis([0,1,-100,20]);

figure(2);

plot(W*fs/(2*pi),abs(H));grid on;

xlabel('频率/Hz');

ylabel('幅值');

n=0:199;t=n/fs;

x=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);

figure(3);

subplot(311);plot(t,x);axis([0,0.01,-5,5]);

title('输入信号');grid on;

y=filter(b,a,x);

subplot(312);stem(y,'.');title('输出序列');grid on;

ya=y*sinc(fs*(ones(length(n),1)*t-(n/fs)'*ones(1,length(t))));

subplot(313);plot(t,ya);axis([0,0.01,-3,3]);title('输出波形');grid on;

t=(0:100)/fs;

figure(4)；

fs=1.5*10000;

n=(0:100)/fs;

f=sin(2*pi*400*t)+3*sin(2*pi*3000*t)+2*sin(2*pi*5000*t);

y=fftfilt(b,x);

[H1,f1]=freqz(f,[1]);

[H2,f2]=freqz(y,[1]);

f1=f1/pi*fs/2;

f2=f2/pi*fs/2;