数字信号处理matlab实验6 无限冲激响应数字滤波器设计

实验目的：

掌握双线性变换法及脉冲相应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。

实验原理：

   在MATLAB中，可以用下列函数辅助设计IIR数字滤波器：1)利用buttord和cheblord可以确定低通原型巴特沃斯和切比雪夫滤波器的阶数和截止频率；2)[num,den]=butter（N,Wn）(巴特沃斯)和[num,den]=cheby1（N,Wn）,[num,den]=cheby2（N,Wn）(切比雪夫1型和2型)可以进行滤波器的设计；3）lp2hp，lp2bp，lp2bs可以完成低通滤波器到高通、带通、带阻滤波器的转换；4)使用bilinear可以对模拟滤波器进行双线性变换，求得数字滤波器的传输函数系数；5)利用impinvar可以完成脉冲响应不变法的模拟滤波器到数字滤波器的转换。

例3-1 设采样周期T=250μs（采样频率fs =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为fc =1kHz。 

           [B,A]=butter(3,2*pi*1000,'s'); 

           [num1,den1]=impinvar(B,A,4000); 

           [h1,w]=freqz(num1,den1); 

           [B,A]=butter(3,2/0.00025,'s'); 

           [num2,den2]=bilinear(B,A,4000); 

           [h2,w]=freqz(num2,den2); 

           f=w/pi*2000; 

           plot(f,abs(h1),'-.',f,abs(h2),'-'); 

           grid; 

           xlabel('频率/Hz ') 

           ylabel('幅值/dB') 

   程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图1给出了这两种设计方法所得到的频响，虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。同时，也看到双线性变换法，在z=-1即Ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在ω=∞处的三阶传输零点通过映射形成的。 

例2  设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。   

            wc=2*1000*tan(2*pi*400/(2*1000)); 

            wt=2*1000*tan(2*pi*317/(2*1000)); 

            [N,wn]=cheb1ord(wc,wt,0.5,19,'s'); 

            [B,A]=cheby1(N,0.5,wn,'high','s'); 

            [num,den]=bilinear(B,A,1000); 

            [h,w]=freqz(num,den); 

            f=w/pi*500; 

            plot(f,20*log10(abs(h))); 

            axis([0,500,-80,10]); 

            grid; 

            xlabel('') 

             ylabel('幅度/dB') 

图2给出了MATLAB计算的结果。     

 例3 设计一巴特沃兹带通滤波器，其３dB边界频率分别为f2=110kHz和f1=90kHz，在阻带f3 = 120kHz处的最小衰减大于１０dB，采样频率fs=400kHz。             

w1=2*400*tan(2*pi*90/(2*400)); 

                  w2=2*400*tan(2*pi*110/(2*400)); 

                  wr=2*400*tan(2*pi*120/(2*400)); 

                  [N,wn]=buttord([w1 w2],[0 wr],3,10,'s'); 

                  [B,A]=butter(N,wn,'s'); 

                  [num,den]=bilinear(B,A,400); 

                  [h,w]=freqz(num,den); 

                  f=w/pi*200; 

                  plot(f,20*log10(abs(h))); 

                  axis([40,160,-30,10]); 

                  grid; 

                  xlabel('频率/kHz') 

                  ylabel('幅度/dB') 

图3给出了MATLAB计算的结果， 

 例4 一数字滤波器采样频率fs = 1kHz，要求滤除100Hz的干扰，其３dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为                    

               w1=95/500; 

               w2=105/500; 

               [B,A]=butter(1,[w1, w2],'stop'); 

               [h,w]=freqz(B,A); 

               f=w/pi*500; 

               plot(f,20*log10(abs(h))); 

               axis([50,150,-30,10]); 

               grid; 

               xlabel('频率/Hz') 

               ylabel('幅度/dB')

            图3.4为MATLAB的计算结果              

实验内容：  利用MATLAB编程设计一个数字带通滤波器，指标要求如下：

通带边缘频率：，，通带峰值起伏：。

阻带边缘频率：，，最小阻带衰减：。

   分别用脉冲响应不变法和双线性变换法进行IIR数字滤波器的设计。

一、双线性变换法（巴特沃兹原型）：

ws1=2*8000*tan(0.3*pi/2);

ws2=2*8000*tan(0.8*pi/2);

wp1=2*8000*tan(0.45*pi/2);

wp2=2*8000*tan(0.65*pi/2);

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=buttord(wp,ws,Rp,Rs,'s');

[num,den]=butter(N,Wn,'s');

[B,A]=bilinear(num,den,8000); 

[h,w]=freqz(B,A);

f=w/pi*4000;

subplot (2,1,1);

plot(f,20*log10(abs(h))); 

axis([0,4000,-60,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

subplot(2,1,2);

plot(f,angle(h));

grid; xlabel('频率/Hz') ;ylabel('相位');

二、脉冲响应不变法（巴特沃兹原型）：

fs=8000;

ws1=0.3*pi*fs;ws2=0.8*pi*fs;

wp1=0.45*pi*fs;wp2=0.65*pi*fs;

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=buttord(wp,ws,Rp,Rs,'s');

[num,den]=butter(N,Wn,'s');

[B,A]=impinvar(num,den,8000); 

[h,w]=freqz(B,A);

f=w/pi*4000;

subplot(2,1,1);

plot(f,20*log10(abs(h))); 

axis([0,4000,-80,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

subplot(2,1,2);

plot(f,angle(h));

grid; xlabel('频率/Hz') ;ylabel('相位');

fs=2;

ws1=2*fs*tan(0.3*pi/2);

ws2=2*fs*tan(0.8*pi/2);

wp1=2*fs*tan(0.45*pi/2);

wp2=2*fs*tan(0.65*pi/2);

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=buttord(wp,ws,Rp,Rs,'s');

[num,den]=butter(N,Wn,'s');

[B,A]=bilinear(num,den,fs); 

[h1,w]=freqz(B,A);

ws11=0.3*pi*fs;ws22=0.8*pi*fs;

wp11=0.45*pi*fs;wp22=0.65*pi*fs;

ws0=[ws11 ws22];wp0=[wp11 wp22];

[N,Wn]=buttord(wp0,ws0,Rp,Rs,'s');

[num,den]=butter(N,Wn,'s');

[B,A]=impinvar(num,den,fs); 

[h2,w]=freqz(B,A);

figure(1)

plot(w/pi,20*log10(abs(h1)),'-.',w/pi,20*log10(abs(h2)),'-');

axis([0,1,-80,0]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

双线性变换法（切比雪夫原型）：

ws1=2*8000*tan(0.3*pi/2);

ws2=2*8000*tan(0.8*pi/2);

wp1=2*8000*tan(0.45*pi/2);

wp2=2*8000*tan(0.65*pi/2);

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');

[num,den]=cheby1(N,1,Wn,'s');

[B,A]=bilinear(num,den,8000); 

[h,w]=freqz(B,A);

f=w/pi*4000;

subplot(2,1,1);

plot(f,20*log10(abs(h))); 

axis([0,4000,-60,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

subplot(2,1,2);

plot(f,angle(h));

grid; xlabel('频率/Hz') ;ylabel('相位');

脉冲响应不变法（切比雪夫原型）：

fs=8000;

ws1=0.3*pi*fs;ws2=0.8*pi*fs;

wp1=0.45*pi*fs;wp2=0.65*pi*fs;

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');

[num,den]=cheby1(N,1,Wn,'s');

[B,A]=impinvar(num,den,8000); 

[h,w]=freqz(B,A);

f=w/pi*4000;

subplot(2,1,1);

plot(f,20*log10(abs(h))); 

axis([0,4000,-90,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

subplot(2,1,2);

plot(f,angle(h));

grid; xlabel('频率/Hz') ;ylabel('相位');

fs=8000;

ws1=2*fs*tan(0.3*pi/2);

ws2=2*fs*tan(0.8*pi/2);

wp1=2*fs*tan(0.45*pi/2);

wp2=2*fs*tan(0.65*pi/2);

ws=[ws1 ws2];wp=[wp1 wp2];

Rp=1;Rs=40;

[N,Wn]=cheb1ord(wp,ws,Rp,Rs,'s');

[num,den]=cheby1(N,1,Wn,'s');

[B,A]=bilinear(num,den,fs); 

[h1,w]=freqz(B,A);

ws11=0.3*pi*fs;ws22=0.8*pi*fs;

wp11=0.45*pi*fs;wp22=0.65*pi*fs;

ws0=[ws11 ws22];wp0=[wp11 wp22];

[N,Wn]=cheb1ord(wp0,ws0,Rp,Rs,'s');

[num,den]=cheby1(N,1,Wn,'s');

[B,A]=impinvar(num,den,fs); 

[h2,w]=freqz(B,A);

figure(1)

plot(w/pi,20*log10(abs(h1)),'-.',w/pi,20*log10(abs(h2)),'-');

axis([0,1,-80,10]);

grid; xlabel('频率/Hz') ;ylabel('幅度/dB');

实验要求：

    给出IIR数字滤波器参数和滤波器的冲激响应，绘出它们的幅度和相位频响曲线，讨论它们各自的实现形式和特点。

    双线性变换法通过将数字频率的取值范围从0到对应到模拟频率，也就对应于模拟域中所有可能的频率值。双线性变换法不会出现频率混叠，但非线性关系却导致数字滤波器的频率响应不能逼真地模仿模拟滤波器的频率响应。脉冲响应不变法通过选择满足设计要求的模拟滤波器冲激响应h（t）的采样值的数字脉冲响应h[n]得到的被采样的冲激响应将给出与原模拟滤波器非常相近的滤波器形状。由于该方法不可避免的要发生频率混叠现象，所以只适合设计低通和带通滤波器。从实验结果可以看出：双线性变换法所设计的巴特沃兹滤波器最符合设计指标，而用脉冲响应不变法设计的滤波器（无论是巴特沃兹还是切比雪夫）有一定的误差，主要是由于混叠所引起。