数字滤波器设计和应用

1、实验四 数字滤波器设计及应用201030021717 电信转专业班02号 王泽军一、实验内容（一）1，用窗函数法设计如下三种类型的FIR滤波器，滤波器的性能指标： 低通滤波器：通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB； 高通滤波器：通带截止频率5000Hz，阻带截止频率4800Hz，通带波纹1dB，最小阻带衰减60dB； 带通滤波器：通带截止频率fp1=1200Hz, fp2=3000Hz，阻带截止频率fs1=1000Hz， fs2=3200Hz, 通带波纹1dB，最小阻带衰减60dB。 (1)设计原理令希望设计的滤波器的传输函数是，是与其对应的

2、单位脉冲响应。一般情况下，由求出，然后由Z变换求出滤波器的系统函数。但是通常在边界频率处有不连续点，这使得是无限长的非因果序列，所以实际是不能实现的。为了构造一个长度为N的线性相位滤波器，可以将截取一段来近似，并且根据线性相位的特点，需要保证截取后的序列关于对称。设截取的一段为，则其中称为矩形窗函数。当的对称中心点取值为时，就可以保证所设计的滤波器具有线性相位。常用的窗函数除了矩形窗还有三角形窗、汉宁窗、哈明窗和布莱克曼窗。在MATLAB中提供了相应的子程序来实现这些窗函数。 w=boxcar(N)，在数组w中产生N点的矩形窗函数。 w=bartlett(N)，在数组w中产生N点的三角形窗函数

3、。 w=hanning(N)，在数组w中产生N点的汉宁窗函数。 w=hamming(N)，在数组w中产生N点的哈明窗函数。 w=blackman(N)，在数组w中产生N点的布莱克曼窗函数。这5种窗函数的具体性能见表1。表1 种窗函数性能比较窗函数旁瓣峰值dB主瓣宽度(近似过渡带宽)精确过渡带宽阻带最小衰减dB矩形窗-1341.8-21三角形窗-2586.1-25汉宁窗-3186.2-44哈明窗-4186.6-53布莱克曼窗-571211-74(2) FIR滤波器的设计 设计低通FIR滤波器，设计条件：通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB；最小

4、阻带衰减为60dB，对照表1，只能选择布莱克曼窗函数进行设计，假设采样率，则归一化的通带截止角频率。归一化的阻带截止角频率。由此确定FIR低通滤波器的截止频率为，设计的，编写如下MATLAB程序：% 数字低通FIR滤波器的设计clearFs=16000;fp=1000;fs=1200;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;wc=(wp+ws)/2; % 截止频率的确定dw=ws-wp;M=ceil(5.56*pi/dw); % 布莱克曼窗的参数Mb=blackman(2*M+1);hlp=fir1(2*M,wc/pi,b);wvtool(hlp)程序运行后可以得到低通滤波器的

5、时域和频域特性如下图1：图1：低通FIR滤波器的时频特性 设计高通FIR滤波器，设计条件：通带截止频率5000Hz，阻带截止频率4800Hz，通带波纹1dB，最小阻带衰减60dB；最小阻带衰减为60dB，所以只能选择布莱克曼窗函数，假设采样率，编写MATLAB程序如下：% 数字高通FIR滤波器的设计clearFs=16000;fp=5000;fs=4800;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;wc=(wp+ws)/2; % 截止频率的确定dw=wp-ws;M=ceil(5.56*pi/dw); % 布莱克曼窗的参数Mb=blackman(2*M+1);hlp=fir1(2*

6、M,wc/pi,'high',b);wvtool(hlp)程序运行后可以得到高通滤波器的时域和频域特性如下图2：图2：高通FIR滤波器的时频特性 设计FIR带通滤波器，设计条件：通带截止频率fp1=1200Hz, fp2=3000Hz，阻带截止频率fs1=1000Hz，fs2=3200Hz,通带波纹1dB，最小阻带衰减60dB。最小阻带衰减为60dB，所以只能选择布莱克曼窗函数，假设采样率，编写MATLAB程序如下：% 数字带通FIR滤波器的设计clearFs=16000;fp1=1200;fp2=3000;fs1=1000;fs2=3200;wp1=2*pi*fp1/Fs;w

7、p2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;wc1=(wp1+ws1)/2;wc2=(wp2+ws2)/2; % 截止频率的确定dw=max(abs(wp1-ws1),abs(wp2-ws2);M=ceil(5.56*pi/dw); % 布莱克曼窗的参数Mb=blackman(2*M+1);fc=wc1,wc2/pi;hlp=fir1(2*M,fc,'high',b);wvtool(hlp)程序运行后可以得到带通滤波器的时域和频域特性如下图3：图3：带通FIR滤波器的时频特性2，再用双线性变换法设计上面的三种IIR滤波器（可以选

8、择butterworth or chebyshey），并用Freqz函数画出各个滤波器的频率响应。答：在实验三中已经详细给出了如何利用双线性变换法设计IIR滤波器，这里不再赘述，只给出相应的程序和结果。(1)分别用butterworth、chebyshey1型和chebyshey2型设计低通IIR滤波器，设计条件：通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB；取采样率。编写MATLAB程序如下：% 用butteworth设计低通IIR滤波器cleara=1;b=60;Fs=8000;fp=1000;fs=1200;wp=2*pi*fp/Fs;ws=2

9、*pi*fs/Fs;Wp=tan(wp/2); % 模拟频带边界频率的确定Ws=tan(ws/2);N,Wn=buttord(Wp,Ws,a,b,'s'); % 确定巴特沃兹低通滤波器的截止频率num,den=butter(N,Wn,'s'); p2,p1= bilinear(num,den,0.5); % 从s域到z域进行双线性变换Hlp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,abs(Hlp)grid onxlabel('f/Hz'),ylabel('幅度')title('用but

10、teworth设计低通IIR滤波器')% 用chebyshey1型设计低通IIR滤波器N1,Wc1=cheb1ord(Wp,Ws,a,b,'s'); % 确定切比雪夫1型滤波器的阶数和截止频率num1,den1=cheby1(N1,a,Wc1,'s'); p4,p3=bilinear(num1,den1,0.5); % 从s域到z域进行双线性变换Hlp1,w1=freqz(p4,p3);figure(2)plot(w1/2/pi*Fs,abs(Hlp1)grid onxlabel('f/Hz'),ylabel('幅度')t

11、itle('用chebyshey1型设计低通IIR滤波器')% 用chebyshey2型设计低通IIR滤波器N2,Wc2=cheb2ord(Wp,Ws,a,b,'s'); % 确定切比雪夫2型滤波器的阶数和截止频率num2,den2=cheby2(N2,b,Wc2,'s');p6,p5=bilinear(num2,den2,0.5); % 从s域到z域进行双线性变换Hlp2,w2=freqz(p6,p5);figure(3)plot(w2/2/pi*Fs,abs(Hlp2)grid onxlabel('f/Hz'),ylabel(

12、'幅度')title('用chebyshey2型设计低通IIR滤波器')程序运行后可以分别得到用butterworth、chebyshey1型和chebyshey2型设计的低通IIR滤波器频域特性如下图4：图4-1：用butterworth设计低通IIR滤波器的频域特性图4-2：用chebyshey1型设计低通IIR滤波器的频域特性图4-3：用chebyshey2型设计低通IIR滤波器的频域特性(2) 分别用butterworth、chebyshey1型和chebyshey2型设计高通IIR滤波器，设计条件：通带截止频率5000Hz，阻带截止频率4800Hz，通

13、带波纹1dB，最小阻带衰减60dB；取采样率。编写MATLAB程序如下：% 用butteworth设计高通IIR滤波器cleara=1;b=60;Fs=16000;fp=5000;fs=4800;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=tan(wp/2); % 模拟频带边界频率的确定Ws=tan(ws/2);N,Wn=buttord(Wp,Ws,a,b,'s'); % 确定巴特沃兹高通滤波器的截止频率num,den=butter(N,Wn,'high','s'); p2,p1= bilinear(num,den,0.5);

14、% 从s域到z域进行双线性变换Hhp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,20*log10(abs(Hhp)grid onxlabel('f/Hz'),ylabel('幅度')title('用butteworth设计高通IIR滤波器')% 用chebyshey1型设计高通IIR滤波器N1,Wc1=cheb1ord(Wp,Ws,a,b,'s'); % 确定切比雪夫1型滤波器的阶数和截止频率num1,den1=cheby1(N1,a,Wc1,'high','s'

15、); p4,p3= bilinear(num1,den1,0.5); % 从s域到z域进行双线性变换Hhp1,w1=freqz(p4,p3);figure(2)plot(w1/2/pi*Fs,abs(Hhp1)grid onxlabel('f/Hz'),ylabel('幅度')title('用chebyshey1型设计高通IIR滤波器')% 用chebyshey2型设计高通IIR滤波器N2,Wc2=cheb2ord(Wp,Ws,a,b,'s'); % 确定切比雪夫1型滤波器的阶数和截止频率num2,den2=cheby2(N2,b

16、,Wc2,'high','s'); p6,p5=bilinear(num2,den2,0.5); % 从s域到z域进行双线性变换Hhp2,w2=freqz(p6,p5);figure(3)plot(w2/2/pi*Fs,abs(Hhp2)grid onxlabel('f/Hz'),ylabel('幅度')title('用chebyshey2型设计高通IIR滤波器')程序运行后可以分别得到用butterworth、chebyshey1型和chebyshey2型设计的低通IIR滤波器频域特性如下图5：图5-1：用but

17、terworth设计高通IIR滤波器的频域特性图5-2：用chebyshey1型设计高通IIR滤波器的频域特性图5-3：用chebyshey2型设计高通IIR滤波器的频域特性(3) 用butterworth设计带通IIR滤波器，设计条件：通带截止频率fp1=1200Hz, fp2=3000Hz，阻带截止频率fs1=1000Hz，fs2=3200Hz,通带波纹1dB，最小阻带衰减60dB。设计思路是：先用butterworth设计低通的IIR滤波器，然后使用lp2bp函数完成从低通到带通的转换。取采样率，编写MATLAB程序如下：clearFs=8000;a=1;b=60;fp1=1200;fp

18、2=3000;fs1=1000;fs2=3200;wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Wp=Wp1*Wp2;Ws=Ws1*Ws2;Bw=Wp2-Wp1;if Wp>Ws Ws1=Wp/Ws2;end;WWp=1;WWs=(Wp-Ws12)/(Bw*Ws1);N,Wn=buttord(WWp,WWs,a,b,'s');B A=butter(N,Wn,'s&#

19、39;);num,den=lp2bp(B,A,sqrt(Wp),Bw);p2,p1=bilinear(num,den,0.5);Hbp,w=freqz(p2,p1);figure(1)plot(w/2/pi*Fs,abs(Hbp)grid onxlabel('f/Hz'),ylabel('幅度')title('用butteworth设计带通IIR滤波器')程序运行后可以得到用butteworth设计带通IIR滤波器的频域特性如图6所示：图6：用butteworth设计带通IIR滤波器的频域特性二、实验内容（二）选择3个不同频段的信号对其进行频谱分

20、析，根据信号的频谱特征设计3个不同的数字滤波器。将三路信号合成一路信号，分析合成信号的时域和频域特点，然后将合成信号分别通过设计好的3个数字滤波器，分离出原来的三路信号，分析得到的三路信号的时域波形和频谱，与原始的三路信号进行比较说明频分复用的特点。频分复用结构如下图所示：答：(1) 首先产生3个频率分别为50Hz、150Hz、250Hz的正弦信号，实现的MATLAB程序如下：% 产生三个不同频段的信号clccleart=0:.001:0.2;x1=cos(100*pi*t);x2=cos(300*pi*t);x3=cos(500*pi*t);figure(1)subplot(3,1,1)pl

21、ot(x1),axis(0,200,-1,1)xlabel('t'),ylabel('x1')title('频率为50Hz的正弦信号')subplot(3,1,2)plot(x2),axis(0,200,-1,1)xlabel('t'),ylabel('x2')title('频率为150Hz的正弦信号')subplot(3,1,3)plot(x3),axis(0,200,-1,1)xlabel('t'),ylabel('x3')title('频率为250Hz的

22、正弦信号')程序运行后可以得到3个不同频率的正弦信号如图7：图7：3个不同频段信号的时域波形 (2) 对3个不同频段的正弦信号进行频域分析，实现的MATLAB程序如下：% 对三个不同频率信号进行频谱分析N=512;X1=abs(fft(x1,N);X2=abs(fft(x2,N);X3=abs(fft(x3,N);f=1000*(0:255)/512;figure(2)subplot(3,1,1)plot(f,X1(1:256),axis(0,500,0,100)xlabel('f/Hz'),ylabel('X1')title('50Hz正弦信号

23、的幅频特性')subplot(3,1,2)plot(f,X2(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X2')title('150Hz正弦信号的幅频特性')subplot(3,1,3)plot(f,X3(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X3')title('250Hz正弦信号的幅频特性')程序运行后可以得到3个不同频率的正弦信号的幅频特性如图8：图8：3个不同频段信号的幅频特性

24、 (3) 将三路不同频段正弦信号合成一路信号并进行时频分析，实现的MATLAB程序如下：% 将三路信号合成一路信号并进行时频分析x=x1+x2+x3;figure(3)subplot(2,1,1)plot(t,x)xlabel('t'),ylabel('x')title('合成信号的时域波形')X=abs(fft(x,N);subplot(2,1,2)plot(f,X(1:256),axis(0,500,0,120)xlabel('f/Hz'),ylabel('X')title('合成信号的幅频特性'

25、;)程序运行后可以得到合成信号的时域和频域特性如图9：图9：合成信号的时域和频域特性 (4) 设计三个数字IIR滤波器，它们分别是低通、带通和高通。其中50Hz在低通IIR滤波器的通带内，150Hz、250Hz在其阻带内；150Hz在带通滤波器的通带内，50Hz、250Hz在其阻带内；250Hz在高通滤波器的通带内，50Hz、150Hz在其阻带内。实现的MATLAB程序如下：% 设计三个数字滤波器(IIR或FIR)% 设计IIR数字低通滤波器% 设计条件Fs=600;fp=70;fs=100;wp=2*pi*fp/Fs;ws=2*pi*fs/Fs;Wp=tan(wp/2);Ws=tan(ws/

26、2);a=1;b=40;% 用butterworth设计低通IIR滤波器N1,Wc1=buttord(Wp,Ws,a,b,'s');num1,den1=butter(N1,Wc1,'s');p2,p1=bilinear(num1,den1,0.5);H1,w1=freqz(p2,p1);figure(4)plot(w1/2/pi*Fs,abs(H1)grid onxlabel('f/Hz'),ylabel('幅度')title('用butterworth设计低通IIR滤波器')% 设计IIR数字带通滤波器% 设计条

27、件Fs=800;fp1=100;fp2=180;fs1=80;fs2=200;wp1=2*pi*fp1/Fs;wp2=2*pi*fp2/Fs;ws1=2*pi*fs1/Fs;ws2=2*pi*fs2/Fs;a=1;b=30;Wp1=tan(wp1/2);Wp2=tan(wp2/2);Ws1=tan(ws1/2);Ws2=tan(ws2/2);Wp=Wp1*Wp2;Ws=Ws1*Ws2;% 用butterworth设计带通IIR滤波器if Wp>Ws Ws1=Wp/Ws2;end;Bw=Wp2-Wp1;WWp=1;WWs=(Wp-Ws12)/(Bw*Ws1);N2,Wc2=buttord(

28、WWp,WWs,a,b,'s');B A=butter(N2,Wc2,'s');num,den=lp2bp(B,A,sqrt(Wp),Bw);p4,p3=bilinear(num,den,0.5);H2,w2=freqz(p4,p3);figure(5)plot(w2/2/pi*Fs,abs(H2)grid onxlabel('f/Hz'),ylabel('幅度')title('用butterworth设计带通IIR滤波器')axis(0,300,0,1.4)% 设计IIR数字高通滤波器% 设计条件Fs=1000;

29、fs=200;fp=220;ws=2*pi*fs/Fs;wp=2*pi*fp/Fs;a=1;b=60;Wp=tan(wp/2);Ws=tan(ws/2);% 用chebyshey2型设计高通IIR滤波器N3,Wc3=cheb2ord(Wp,Ws,a,b,'s');num,den=cheby2(N3,b,Wc3,'high','s');p6,p5=bilinear(num,den,0.5);H3,w3=freqz(p6,p5);figure(6)plot(w3/2/pi*Fs,abs(H3),axis(0,500,0,1.4)grid onxlab

30、el('f/Hz'),ylabel('幅度')title('用chebyshey2型设计高通IIR滤波器')程序运行后可以得到3个IIR滤波器的频域特性如图10所示：图10-1：用butterworth设计低通IIR滤波器的幅频特性 图10-2：用butterworth设计带通IIR滤波器的幅频特性 图10-3：用chebyshey2型设计高通IIR滤波器的幅频特性 (5) 将合成信号依次通过低通、带通、高通IIR滤波器，滤波后得到了3路信号。实现的MATLAB程序如下：% 将合成信号通过三个IIR滤波器% 用低通IIR滤波器对信号进行滤波y1=filter(p2,p1,x);Y1=abs(f