数字滤波器的设计和应用.docx
《数字滤波器的设计和应用.docx》由会员分享，可在线阅读，更多相关《数字滤波器的设计和应用.docx（24页珍藏版）》请在冰豆网上搜索。
数字滤波器的设计和应用
实验四数字滤波器设计及应用
201030021717电信转专业班02号王泽军
一、实验内容
（一）
1，用窗函数法设计如下三种类型的FIR滤波器，滤波器的性能指标：
低通滤波器：
通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB；
高通滤波器：
通带截止频率5000Hz，阻带截止频率4800Hz，通带波纹1dB，最小阻带衰减60dB；
带通滤波器：
通带截止频率fp1=1200Hz,fp2=3000Hz，阻带截止频率fs1=1000Hz，fs2=3200Hz,通带波纹1dB，最小阻带衰减60dB。
（1）设计原理
令希望设计的滤波器的传输函数是
，
是与其对应的单位脉冲响应。
一般情况下，由
求出
，然后由Z变换求出滤波器的系统函数。
但是通常
在边界频率处有不连续点，这使得
是无限长的非因果序列，所以实际是不能实现的。
为了构造一个长度为N的线性相位滤波器，可以将
截取一段来近似，并且根据线性相位的特点，需要保证截取后的序列关于
对称。
设截取的一段为
，则
其中
称为矩形窗函数。
当
的对称中心点取值为
时，就可以保证所设计的滤波器具有线性相位。
常用的窗函数除了矩形窗还有三角形窗、汉宁窗、哈明窗和布莱克曼窗。
在MATLAB中提供了相应的子程序来实现这些窗函数。
w=boxcar（N），在数组w中产生N点的矩形窗函数。
w=bartlett（N），在数组w中产生N点的三角形窗函数。
w=hanning（N），在数组w中产生N点的汉宁窗函数。
w=hamming（N），在数组w中产生N点的哈明窗函数。
w=blackman（N），在数组w中产生N点的布莱克曼窗函数。
这5种窗函数的具体性能见表1。
表1
种窗函数性能比较
窗函数
旁瓣峰值
dB
主瓣宽度
（近似过渡带宽）
精确过渡带宽
阻带最小衰减
dB
矩形窗
-13
4
1.8
-21
三角形窗
-25
8
6.1
-25
汉宁窗
-31
8
6.2
-44
哈明窗
-41
8
6.6
-53
布莱克曼窗
-57
12
11
-74
（2）FIR滤波器的设计
设计低通FIR滤波器，设计条件：
通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB；
最小阻带衰减为60dB，对照表1，只能选择布莱克曼窗函数进行设计，假设采样率
，则归一化的通带截止角频率
。
归一化的阻带截止角频率
。
由此确定FIR低通滤波器的截止频率为
，设计的
，编写如下MATLAB程序：
%%%数字低通FIR滤波器的设计
clear
Fs=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）;%布莱克曼窗的参数M
b=blackman（2*M+1）;
hlp=fir1（2*M,wc/pi,b）;
wvtool（hlp）
程序运行后可以得到低通滤波器的时域和频域特性如下图1：
图1：
低通FIR滤波器的时频特性
设计高通FIR滤波器，设计条件：
通带截止频率5000Hz，阻带截止频率4800Hz，通带波纹1dB，最小阻带衰减60dB；
最小阻带衰减为60dB，所以只能选择布莱克曼窗函数，假设采样率
，编写MATLAB程序如下：
%%%数字高通FIR滤波器的设计
clear
Fs=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）;%布莱克曼窗的参数M
b=blackman（2*M+1）;
hlp=fir1（2*M,wc/pi,'high',b）;
wvtool（hlp）
程序运行后可以得到高通滤波器的时域和频域特性如下图2：
图2：
高通FIR滤波器的时频特性
设计FIR带通滤波器，设计条件：
通带截止频率fp1=1200Hz,fp2=3000Hz，阻带截止频率fs1=1000Hz，fs2=3200Hz,通带波纹1dB，最小阻带衰减60dB。
最小阻带衰减为60dB，所以只能选择布莱克曼窗函数，假设采样率
，编写MATLAB程序如下：
%%%数字带通FIR滤波器的设计
clear
Fs=16000;
fp1=1200;fp2=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;
wc1=（wp1+ws1）/2;wc2=（wp2+ws2）/2;%截止频率的确定
dw=max（abs（wp1-ws1）,abs（wp2-ws2））;
M=ceil（5.56*pi/dw）;%布莱克曼窗的参数M
b=blackman（2*M+1）;
fc=[wc1,wc2]/pi;
hlp=fir1（2*M,fc,'high',b）;
wvtool（hlp）
程序运行后可以得到带通滤波器的时域和频域特性如下图3：
图3：
带通FIR滤波器的时频特性
2，再用双线性变换法设计上面的三种IIR滤波器（可以选择butterworthorchebyshey），并用Freqz函数画出各个滤波器的频率响应。
答：
在实验三中已经详细给出了如何利用双线性变换法设计IIR滤波器，这里不再赘述，只给出相应的程序和结果。
（1）分别用butterworth、chebyshey1型和chebyshey2型设计低通IIR滤波器，设计条件：
通带截止频率1000Hz，阻带截止频率1200Hz，通带波纹1dB，最小阻带衰减60dB；
取采样率
。
编写MATLAB程序如下：
%%%用butteworth设计低通IIR滤波器
clear
a=1;b=60;
Fs=8000;
fp=1000;fs=1200;
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,'s'）;
[p2,p1]=bilinear（num,den,0.5）;%从s域到z域进行双线性变换
[Hlp,w]=freqz（p2,p1）;
figure
（1）
plot（w/2/pi*Fs,abs（Hlp））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用butteworth设计低通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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
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，通带波纹1dB，最小阻带衰减60dB；
取采样率
。
编写MATLAB程序如下：
%%%用butteworth设计高通IIR滤波器
clear
a=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）;%从s域到z域进行双线性变换
[Hhp,w]=freqz（p2,p1）;
figure
（1）
plot（w/2/pi*Fs,20*log10（abs（Hhp）））
gridon
xlabel（'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'）;
[p4,p3]=bilinear（num1,den1,0.5）;%从s域到z域进行双线性变换
[Hhp1,w1]=freqz（p4,p3）;
figure
（2）
plot（w1/2/pi*Fs,abs（Hhp1））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用chebyshey1型设计高通IIR滤波器'）
%%%用chebyshey2型设计高通IIR滤波器
[N2,Wc2]=cheb2ord（Wp,Ws,a,b,'s'）;%确定切比雪夫1型滤波器的阶数和截止频率
[num2,den2]=cheby2（N2,b,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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用chebyshey2型设计高通IIR滤波器'）
程序运行后可以分别得到用butterworth、chebyshey1型和chebyshey2型设计的低通IIR滤波器频域特性如下图5：
图5-1：
用butterworth设计高通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程序如下：
clear
Fs=8000;
a=1;b=60;
fp1=1200;fp2=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;
ifWp>Ws
Ws1=Wp/Ws2;
end;
WWp=1;
WWs=（Wp-Ws1^2）/（Bw*Ws1）;
[N,Wn]=buttord（WWp,WWs,a,b,'s'）;
[BA]=butter（N,Wn,'s'）;
[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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用butteworth设计带通IIR滤波器'）
程序运行后可以得到用butteworth设计带通IIR滤波器的频域特性如图6所示：
图6：
用butteworth设计带通IIR滤波器的频域特性
二、实验内容
（二）
选择3个不同频段的信号对其进行频谱分析，根据信号的频谱特征设计3个不同的数字滤波器。
将三路信号合成一路信号，分析合成信号的时域和频域特点，然后将合成信号分别通过设计好的3个数字滤波器，分离出原来的三路信号，分析得到的三路信号的时域波形和频谱，与原始的三路信号进行比较说明频分复用的特点。
频分复用结构如下图所示：
答：
（1）首先产生3个频率分别为50Hz、150Hz、250Hz的正弦信号，实现的MATLAB程序如下：
%%%产生三个不同频段的信号
clc
clear
t=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）
plot（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的正弦信号'）
程序运行后可以得到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正弦信号的幅频特性'）
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个不同频段信号的幅频特性
（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（'合成信号的幅频特性'）
程序运行后可以得到合成信号的时域和频域特性如图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/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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用butterworth设计低通IIR滤波器'）
%设计IIR数字带通滤波器
%设计条件
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滤波器
ifWp>Ws
Ws1=Wp/Ws2;
end;
Bw=Wp2-Wp1;
WWp=1;
WWs=（Wp-Ws1^2）/（Bw*Ws1）;
[N2,Wc2]=buttord（WWp,WWs,a,b,'s'）;
[BA]=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））
gridon
xlabel（'f/Hz'）,ylabel（'幅度'）
title（'用butterworth设计带通IIR滤波器'）
axis（[0,300,0,1.4]）
%设计IIR数字高通滤波器
%设计条件
Fs=1000;
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]）
gridon
xlabel（'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（fft（y1,N））;
figure（7）
subplot（2,1,1）
plot（y1）,axis（[0,200,-1,1]）
xlabel（'t'）,ylabel（'y1'）
title（'经过低通IIR滤波器后的时域信号'）
subplot（2,1,2）
plot（f,Y1（1:
256））
xlabel（'f/Hz'）,ylabel（'Y1'）
title（'经过低通IIR滤波器后的频域信号'）
%用带通IIR滤波器对信号进行滤波
y2=filter（p4,p3,x）;
Y2=abs（fft（y2,N））;
figure（8）
subplot（2,1,1）
plot（y2）,axis（[0,200,-1,1]）
xlabel（'t'）,ylabel（'y2'）
title（'经过带通IIR滤波器后的时域信号'）
subplot（2,1,2）
plot（f,Y2（1:
256））
xlabel（'f/Hz'）,ylabel（'Y2'）
title（'经过带通IIR滤波器后的频域信号'）
%用高通IIR滤波器对信号进行滤波
y3=filter（p6,p5,x）;
Y3=abs（fft（y3,N））;
figure（9）
subplot（2,1,1）
plot（y3）,axis（[0,200,-1,1]）
xlabel（'t'）,ylabel（'y3'）
titl
升级会员 