数字滤波器的设计和应用.docx
《数字滤波器的设计和应用.docx》由会员分享,可在线阅读,更多相关《数字滤波器的设计和应用.docx(24页珍藏版)》请在冰豆网上搜索。
![数字滤波器的设计和应用.docx](https://file1.bdocx.com/fileroot1/2023-2/21/64dd11e5-a484-4f04-962c-7c9606eaf2cf/64dd11e5-a484-4f04-962c-7c9606eaf2cf1.gif)
数字滤波器的设计和应用
实验四数字滤波器设计及应用
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