数字信号处理实验报告
姓名:寇新颖 学号:20100304026 专业:电子信息科学与技术
实验四 无限冲激响应(IIR)数字滤波器的设计
一、实验目的
1.掌握双线性变换法及冲激响应不变法设计IIR数字滤波器的具体设计方法及其原理,熟悉用双线性变换法及冲激响应不变法设计低通IIR数字滤波器的计算机编程。
2.观察双线性变换及冲激响应不变法设计的滤波器的频域特性,了解双线性变换法及冲激响应不变法的特点。
3.熟悉Butterworth滤波器的频率特性。
二、实验原理
1.利用模拟滤波器设计IIR数字滤波器方法
(1)根据所给出的数字滤波器性能指标计算出相应的模拟滤波器的设计指标。
(2)根据得出的滤波器性能指标设计出相应的模拟滤波器的系统函数H(S)。
(3)根据得出的模拟滤波器的系统函数H(S),经某种变换得到对该模拟滤波器相应的数字仿真系统——数字滤波器。
将模拟滤波器转换成数字滤波器的实质是,用一种从s平面到z平面的映射函数将Ha(s)转换成H(z)。对这种映射函数的要求是:(1) 因果稳定的模拟滤波器转换成数字滤波器,仍是因果稳定的。 (2)数字滤波器的频率响应模仿模拟滤波器的频响,s平面的虚轴映射z平面的单位圆,相应的频率之间成线性关系。冲激响应不变法和双线性变换法都满足如上要求。
2.冲激响应不变法
用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应ha(t),让h(n)正好等于ha(t)的采样值,即h(n)=ha(nT),其中T为采样间隔。
3.双线性变换法
s平面与z平面之间满足以下映射关系:
s平面的虚轴单值地映射于z平面的单位圆上,s平面的左半平面完全映射到z平面的单位圆内。双线性变换不存在混叠问题。
双线性变换时一种非线性变换,这种非线性引起的幅频特性畸变可通过预畸而得到校正。
以低通数字滤波器为例,将设计步骤归纳如下:
(1)确定数字滤波器的性能指标:通带临界频率fp、阻带临界频率fs;通带内的最大衰减Ap;阻带内的最小衰减As;
(2)确定相应的数字角频率,ωp=2πfp;ωs=2πfs;
(3)计算经过预畸的相应模拟低通原型的频率,;
(4)根据Ωp和Ωs计算模拟低通原型滤波器的阶数N,并求得低通原型的传递函数Ha(s);
(5)用上面的双线性变换公式代入Ha(s),求出所设计的传递函数H(z);
(6)分析滤波器特性,检查其指标是否满足要求。
三、主要实验仪器及材料
微型计算机、Matlab7.0教学版、TC编程环境。
四、实验内容
1.知识准备
在实验前复习数字信号处理理论课中有关滤波器设计的知识,认真阅读本实验的原理部分。
2.编制用冲激响应不变法和用双线性变换法设计IIR数字滤波器的程序。
例1 采用冲激响应不变法设计一个巴特沃斯数字低通滤波器,要求
,滤波器采样频率。
例2 采用双线性变换法设计一个巴特沃斯数字低通滤波器,要求:,滤波器采样频率。
五、实验程序及结果
例1 冲激响应不变法设计
程序如下:
wp=0.25*pi;
ws=0.4*pi;
Rp=1;As=15;
ripple=10^(-Rp/20);%计算通带衰减对应的幅度值
Attn=10^(-As/20);%计算阻带衰减对应的幅度值
Fs=2000;T=1/Fs;
Omgp=wp*Fs;Omgs=ws*Fs;
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s') %计算阶数N和截止频率
[z0,p0,k0]=buttap(n);%设计归一化的巴特沃斯模拟原型滤波器
ba1=k0*real(poly(z0));%求原型滤波器系数b
aa1=real(poly(p0));%求原型滤波器系数a
[ba,aa]=lp2lp(ba1,aa1,Omgc);变换为模拟低通滤波器
[bd,ad]=impinvar(ba,aa,Fs)
%[C,B,A]=dir2par(bd,ad)
%[r,p,k]=residuez(bd,ad)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('H模');title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);
ylabel('dB');title('频率\pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
bd=[0.0000 0.0031 0.0419 0.0569 0.0125 0.0003 0.0000];
ad=[1.0000 -2.5418 3.1813 -2.3124 1.0072 -0.2457 0.0260];
subplot(2,2,4),zplane(bd,ad); title('零极图');
例2 双线性变换法设计
程序如下:
wp=0.25*pi;
ws=0.4*pi;
Rp=1;As=15;
ripple=10^(-Rp/20);
Attn=10^(-As/20);
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') %计算阶数n和截止频率
[ba1,aa1]=butter(n,Omgc,'s');%模拟低通滤波器系数b,a
%[z0,p0,k0]=buttap(n);
%ba=k0*real(poly(z0));
%aa=real(poly(p0));
%[ba1,aa1]=lp2lp(ba,aa,Omgc);
[bd,ad]=bilinear(ba1,aa1,Fs)%用双线性变换法求数字滤波器系数b,a
[sos,g]=tf2sos(bd,ad)
[H,w]=freqz(bd,ad);
dbH=20*log10((abs(H)+eps)/max(abs(H)));
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('H模');title('幅度响应');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('\phi');title('相位响应');axis([0,1,-1,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);
ylabel('dB');title('频率\pi)');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);title('零极图');