[H,w]=freqz(bz,az);%计算0~pi上的响应
dbHx=-20*log10(abs(Hx)/max(abs(H)))%归一化并求dB dbH=20*log10(abs(H)/max(abs(H)));%归一化的频率响应 plot(w/pi,dbH);%画图
xlabel('数字频率域频率\\Omega/\\pi');ylabel('dB'); axis([0,1,-150,5]); grid
运行,得
N = 3;
OmegaC = 1.0e+04 * 0.3430 2.7528 w0 = 0.3142 1.8850 0.6283 1.2566
bz =0.0946 -0.3508 0.7174 -0.8802 0.7174 -0.3508 0.0946 az = 1.0000 -1.4600 0.3178 -0.3506 0.6885 0.2289 -0.3824 dbHx = 0.2548 3.0000 39.8925 39.8925
0-50dB-100-15000.10.20.30.40.50.60.7数字频率域频率?/?0.80.91
图3双线性变换法设计IIR带阻滤波器
由程序返回得到的数值可以得知,这是一个2阶带阻滤波器。转换为模拟原型滤波器的频率见OmegaC返回的值。数字滤波器的边界频率见w0返回的值。 由bz,az返回值可得
0.0946-0.3508z-1+0.7174z-2-0.8802z-3?0.7174z-4?0.3508z-5?0.0946z-6H(z)?
1?1.4600z-1?0.3718z-2?0.3506z-3?0.6885z-4?0.2289z-5?0.3824z-6由dbHx返回值可知在四个边界频率处幅度响应大小: 对应到模拟频率,即有
f/Hz 500 0.2548 3000 3.0000 1000 39.8925 2000 39.8925 H(f) 满足题目要求的设计指标,通带波纹不大于3dB, 截止频率处衰减大于30dB。
7.15
设计一个数字切比雪夫Ⅱ型带阻滤波器,给定指标为: (4) 衰减As?30dB,当10kHz?f?20kHz (5) 波纹Rp?3dB,当f?5kHz,f?30kHz (6) 抽样频率fs?100kHz
试用双线性变换法进行设计,最后写出H(z)的表达式,并画出系统的幅频响应特性 (dB)。设计时请先想一想,这一题和上一题有什么相似处。由此应该得出什么结论。解:
分析:这一题的各频率指标要求均是7.14的10倍,抽样频率也是10倍。此时对应的数字频率是一样的。如果用切比雪夫Ⅰ型滤波器,双线性变换法设计出来的与7.14中的肯定是一样的(见图4).这说明只要给定的各数字频率参数一样,用同种方法设计出来的数字滤波器是一样的。程序如下:
%设计数字切比雪夫Ⅱ型带阻滤波器 %双线性变换法 clc;clear all
OmegaP1=2*pi*5000;OmegaP2=2*pi*30000;%带通截止频率 OmegaS1=2*pi*10000;OmegaS2=2*pi*20000;%阻带截止频率 Rp=3;%通带波纹dB As=30;%阻带衰减dB
Fs=100*10^3;%抽样频率10khz
OmegaP=[OmegaP1,OmegaP2];OmegaS=[OmegaS1,OmegaS2]; wp=OmegaP/Fs;ws=OmegaS/Fs;%等效数字频率
OmegaP_t=2*Fs*tan(wp/2);OmegaS_t=2*Fs*tan(ws/2);
[N,OmegaC]=cheb2ord(OmegaP_t,OmegaS_t,Rp,As,'s')ˉ阶数和截至频率 [b,a]=cheby2(N,As,OmegaS,'stop','s');ˉ系统函数的分子分母 [bz,az]=bilinear(b,a,Fs)%双线性变换法F to DF w0=[wp,ws];%四个频点
Hx=freqz(bz,az,w0);%计算两个频点上对应的幅度响应 [H,w]=freqz(bz,az);%计算0~pi上的响应
dbHx=-20*log10(abs(Hx)/max(abs(H)))%归一化并求dB dbH=20*log10(abs(H)/max(abs(H)));%归一化的频率响应 plot(w/pi,dbH);%画图
xlabel('数字频率域频率\\Omega/\\pi');ylabel('dB'); axis([0,1,-150,5]); grid
0-50dB-100-15000.10.20.30.40.50.60.7数字频率域频率?/?0.80.91
图4双线性变换法设计IIR带阻滤波器(fs=100kHz)
7.17
要求设计一个数字带通滤波器,其抽样频率fs=25kHz,通达截止频率为
fp1=5kHz,fp2=7kHz,通带衰减Rp=0.5dB,阻带截止频率为fst1=3.5kHz, fst2=8.5kHz,阻带衰减为As=45dB。
(1) 利用MATLAB工具箱中ellipord及ellip设计椭圆函数滤波器; (2) 利用MATLAB工具箱中cheblord及cheb1设计切比雪夫Ⅰ型滤波器; (3) 利用MATLAB工具箱中cheb2ord及cheb2设计切比雪夫Ⅱ型滤波器; (4) 利用buttord及butter设计巴特沃思型滤波器。
要求每种设计都给出系统函数并画出幅频特性(dB)、相频特性以及单位冲激响应。 解:
(1) 利用MATLAB工具箱中ellipord及ellip设计椭圆函数滤波器:
核心语句为:
%设计椭圆带通滤波器 %双线性变换法 clc;clear all
OmegaP1=2*pi*5000;OmegaP2=2*pi*7000;%带通截止频率 OmegaS1=2*pi*3500;OmegaS2=2*pi*8500;%阻带截止频率 Rp=0.5;%通带波纹dB As=45;%阻带衰减dB
Fs=35*10^3;%抽样频率35khz
OmegaP=[OmegaP1,OmegaP2];OmegaS=[OmegaS1,OmegaS2]; wp=OmegaP/Fs/pi;ws=OmegaS/Fs/pi;%等效数字频率 w0=[wp,ws];%四个频点
[N,wc]=ellipord(wp,ws,Rp,As)?阶数和截至频率
[bz,az]=ellip(N,Rp,As,wc)?系统函数的分子分母 Hx=freqz(bz,az,w0*pi);%计算四个频点上对应的幅度响应 [H,w]=freqz(bz,az);%计算0~1上的幅频响应
dbHx=-20*log10(abs(Hx)/max(abs(H)))%归一化并求dB
dbH=20*log10(abs(H)/max(abs(H)));%归一化的频率响应
pha=unwrap(angle(H));%计算0~1上的相频响应
figure(1)
plot(w/pi,dbH);%画幅频图
xlabel('数字频率域频率\\Omega/\\pi');ylabel('dB'); axis([0,1,-80,5]); grid figure(2)
plot(w/pi,pha);%画相频图
xlabel('数字频率域频率\\Omega/\\pi');ylabel('Phase\\degree'); % axis([0,1,-150,5]); grid
disp('系统传递函数H(z)'); printsys(bz,az,'z'); figure(3)
h=dimpulse(bz,az); stem(h)
xlabel('n');ylabel('Impulse response');