数字信号处理实验报告.docx
《数字信号处理实验报告.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验报告.docx(13页珍藏版)》请在冰豆网上搜索。
![数字信号处理实验报告.docx](https://file1.bdocx.com/fileroot1/2022-11/21/6c7038f8-02de-4afc-9a34-dd0d8e921c6a/6c7038f8-02de-4afc-9a34-dd0d8e921c6a1.gif)
数字信号处理实验报告
课程名称:
数字信号处理实验名称:
IIR数字滤波器设计
班级:
姓名:
1、实验目的
1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法;
2、熟悉用脉冲响应不变法设计IIR数字滤波器的原理与方法。
2、实验内容
1、已知低通滤波器的指标为:
通带边缘频率:
0.4π,Ap=0.5dB;
阻带边缘频率:
0.6π,As=50dB;
a、采用脉冲响应不变法设计巴特沃斯,T=1.画出幅度响应和数字滤波器的脉冲响应h(n);
b、采用脉冲响应不变法设计巴特沃斯,T=1.画出幅度响应和数字滤波器的脉冲响应h(n)。
2、用双线性变换法设计低通滤波器,满足技术指标
wp=0.2π,Ap=0.25dB;
ws=0.4π,As=50dB,
并对方波信号进行滤波,画出滤波前后的波形图并进行简要分析。
3、设计一个数字高通滤波器H(z),它用在结构
xa(t)A/DH(z)D/Aya(t)
中,满足下列要求:
a、采样速率为10kHZ;
b、阻带边缘频率为1.5kHZ,衰减为40dB;
c、通带边缘频率为2kHZ,衰减为3dB;
d、单调的通带和阻带。
4、设计一个带阻滤波器,要求通带上下截止频率为0.8π,0.2π,通带内衰减不大于1dB,阻带起始频率为0.7π,0.4π,阻带内衰减不小于30dB。
设计巴特沃斯带阻滤波器并画出该数字高通滤波器的幅度响应和脉冲响应。
3、实验程序及解释和实验分析及图形
1a、clear;closeall;%清屏
wp=0.4*pi;%通带边缘频率设为0.4*pi
ws=0.6*pi;%阻带边缘频率设为0.6*pi
Ap=0.5;%通带最大衰减为0.5dB
As=50;%阻带最小衰减为50dB
%原型指标的频率逆映射
T=1;%设周期T为1
Fs=1/T;
OmegaP=wp/T;%通带截止频率原型
OmegaS=ws/T;%阻带截止频率原型
[cs,ds]=afd_butt(OmegaP,OmegaS,Ap,As);%模拟巴特沃斯原型滤波器计算
[b,a]=impinvar(cs,ds,Fs);%脉冲响应不变变换
[h,w]=freqz(b,a);%计算数字滤波器的Z变换频率响应
subplot(2,2,1);plot(w/pi,abs(h));%画出幅频响应曲线图
title('幅度响应');%命名为“幅频响应”
grid;%给图形加上网格线
subplot(2,2,2);plot(w/pi,angle(h));%画出相频响应曲线图
title('相位响应');%给图形命名
grid;%加网格线
subplot(2,2,3);
plot(w/pi,20*log(abs(h)));%画出取其对数的幅频响应曲线
title('幅度响应dB');grid;%给图形命名,加网格线
n=[0:
1:
59];%设置n值
imp=[1;zeros(59,1)];%输入矩形波
y=filter(b,a,imp);%对imp进行滤波
subplot(2,2,4);plot(n,y);%画出y的图形
title('脉冲响应');%给图形命名
grid;%加网格线
实验分析:
;用脉冲响应不变法设计低通巴特沃斯滤波器时,因为w=Ωt,w与Ω呈线性关系,所以其相位响应图是呈线性的,如上图的“相位响应”图所示;因为其设计的是低通滤波器,所以会把其高频的波滤掉,而留下低频的,从而在变化的时候表现的平缓些,而不是特别的陡峭,如上图的“脉冲响应”图所示,从而实现低通滤波的功能。
1b、clear;closeall;
wp=0.4*pi;%通带边缘频率设为0.4*pi
ws=0.6*pi;%阻带边缘频率设为0.6*pi
Ap=0.5;%通带最大衰减为0.5dB
As=50;%阻带最小衰减为50dB
%原型指标的频率逆映射
T=1;%把周期T设为1
Fs=1/T;
OmegaP=(2/T)*tan(wp/2);%计算通带截止频率原型
OmegaS=(2/T)*tan(ws/2);%计算阻带截止频率原型
[cs,ds]=afd_butt(OmegaP,OmegaS,Ap,As);%巴特沃斯原型滤波器计算
[b,a]=bilinear(cs,ds,Fs);%双线性不变法变换
[h,w]=freqz(b,a);%计算数字滤波器的Z变换频率响应
subplot(2,2,1);plot(w/pi,abs(h));%画幅频响应曲线图
title('幅度响应');
grid;
subplot(2,2,2);plot(w/pi,angle(h));%画出相频响应曲线图
title('相位响应');grid;
subplot(2,2,3);plot(w/pi,20*log(abs(h)));%画出取其对数的幅频响应曲线
title('幅度响应dB');
axis([0,1,-80,5]);grid;%设置坐标轴,加网格线
n=[0:
1:
59];%设n值
imp=[1;zeros(59,1)];%输入矩形波
y=filter(b,a,imp);%对imp进行滤波
subplot(2,2,4);plot(n,y);%画出y的图形
title('脉冲响应');
grid;
实验分析:
用双线性变换法实现低通巴特沃斯滤波器时,Ω=(2/T)tan(w/2),可以知道w与Ω是呈非线性关系的,所以其图形如上图的“相位响应”图所示,呈现的是不是直线,而是曲线。
由于其设计的是低通滤波器,所以把其高频的部分滤掉了,留下低频的部分,从而在变化的时候表现的平缓些,而不是特别的陡峭,如上图的“脉冲响应”图所示,从而实现低通滤波的功能。
2、wp=0.2*pi;%通带边缘频率设为0.2*pi
ws=0.4*pi;%阻带边缘频率设为0.4*pi
Ap=0.25;%通带最大衰减为0.25dB
As=50;%阻带最小衰减为50dB
%原型指标的频率逆映射
T=1;%把周期T设为1
Fs=1/T;
OmegaP=(2/T)*tan(wp/2);%计算通带截止频率原型
OmegaS=(2/T)*tan(ws/2);%计算阻带截止频率原型
N=ceil((log10((10^(Ap/10)-1)/(10^(As/10)-1)))/
(2*log10(OmegaP/OmegaS)));%确定N值
fprintf('\n***ButterworthFilterOrder=%2.0f\n',N);
OmegaC=OmegaP/((10^(Ap/10)-1)^(1/(2*N)));
wn=2*atan((OmegaC*T)/2);%计算截止频率
wn=wn/pi;
[b,a]=butter(N,wn);%计算截止频率为wn的N阶低通巴特沃斯滤波器
[h,w]=freqz(b,a);%计算数字滤波器的Z变换频率响应
subplot(2,2,1);
plot(w/pi,abs(h));%画幅频响应曲线图
title('幅度响应');grid;
subplot(2,2,3);
plot(w/pi,20*log10(abs(h)));%画出取其对数的幅频响应曲线
title('幅度响应dB');
axis([0,1,-80,5]);grid;%设置坐标轴,加网格线
n=[0:
1:
59];
imp=[1;zeros(59,1)];
y=filter(b,a,imp);%对imp进行滤波
subplot(2,2,4);plot(n,y);%画出y的图形
title('脉冲响应');grid;
n=[0:
100];
t=0.001*n;%fs=1/t=1000HZ
x=2*sin(2*pi*400*t)+5*sin(2*pi*50*t);%输入含有两种频率的波形
x1=ones(1,101);%输入单位抽样序列x1
y1=filter(b,a,x1);%对x1进行滤波
figure
(2);
subplot(2,1,1);
plot(n,x1,'p-');%绘制一条由实线串起来的五角星形线x1
subplot(2,1,2);
plot(n,y1,'p-');%绘制y1曲线
实验分析:
用双线性法设计这个低通滤波器,由实验程序中可知,其fs=1000HZ,当输入的是方波函数时,把其高频的部分滤掉,留下低频部分,所以其滤波后的图形在变化时表现的稍微平缓些,而不是特别的陡峭,如上图的“脉冲响应”图所示,从而实现低通滤波的功能。
3、f=10;fp=2;fs=1.5;
wp=2*pi*fp*(1/f);%设置通带边缘频率
ws=2*pi*fs*(1/f);%设置阻带边缘频率
Ap=3;
As=40;
[N,wn]=buttord(wp/pi,ws/pi,Ap,As);%确定N值和截止频率wn
[b,a]=butter(N,wn,'high');%设计截止频率为wn的高通滤波器
[h,w]=freqz(b,a);%计算数字滤波器的Z变换频率响应
subplot(2,2,1);plot(w/pi,abs(h));%画幅频响应曲线图
title('幅度响应');grid;
subplot(2,2,3);plot(w/pi,20*log10(abs(h)));%画幅频的对数的幅频曲线
title('幅度响应dB');grid;
n=[0:
1:
59];imp=[1;zeros(59,1)];
y=filter(b,a,imp);%对imp进行滤波
subplot(2,2,4);plot(n,y);%画出y的图形
title('脉冲响应');grid;
n=[0:
1000];t=0.0001*n;%fs=1/t=10000HZ
x=2*sin(2*pi*4000*t)+5*sin(2*pi*50*t);%输入混有两种频率的波形
y1=filter(b,a,x);%对x进行滤波
figure
(2);
subplot(2,1,1);
plot(n,x,'p-');%绘制一条由实线串起来的五角星形x
subplot(2,1,2);
plot(n,y1,'p-');%绘制一条由实线串起来的五角星形线y1
实验分析:
这里设计的是一个高通数字滤波器,由实验程序中可知,其fs=10000HZ,而输入信号x=x1+x2,(其中x1=2*sin(2*pi*4000*t),x2=5*sin(2*pi*50*t))是由两种不同频率的sin函数混合而成的,fs/2=5000HZ,f1=4000HZ,相对来说是高频,可以通过,f2=50HZ,相对来说是低频,会被滤波器滤掉,所以当x经过高通滤波器滤波后,只会留下x1的频谱,x2会被滤掉;从上图中可以看出,没经滤波前,其波形的幅度变化范围为7到-7,x1、x2两种频率的波是混合在一起的,经滤波后其波形的幅度范围变,2到-2,只剩下了x1波形,即x2被高通滤波器滤掉了,实现了高通滤波器的功能。
4、clear;closeall;
wp=[0.2*pi,0.8*pi];
ws=[0.4*pi,0.7*pi];
Ap=1;
As=30;
[N,wn]=buttord(wp/pi,ws/pi,Ap,As);%确定N和截止频率wn的值
[b,a]=butter(N,wn,'stop');%设计截止频率为wn的带阻滤器
[h,w]=freqz(b,a);%计算数字滤波器的Z变换频率响应
subplot(2,2,1);
plot(w/pi,abs(h));%画幅频响应曲线图
title('幅度响应');%加图名
grid;
subplot(2,2,3);
plot(w/pi,20*log(abs(h)));%画幅频的对数的幅频曲线
title('幅度响应dB');%加图名
axis([0,1,-80,5]);%设置坐标
grid;%加网格线
n=[0:
1:
59];
imp=[1;zeros(59,1)];
y=filter(b,a,imp);%对imp进行滤波
subplot(2,2,4);
plot(n,y);%画出滤波后的曲线y
title('脉冲响应');%加图名
grid;
n=[0:
100];
t=0.001*n;%fs=1/t=1000HZ
x=2*sin(2*pi*400*t)+5*sin(2*pi*50*t);
x1=ones(1,101);
y=filter(b,a,x);%对x进行滤波,得到y
figure
(2);
subplot(2,1,1);
plot(n,x,'p-');%绘制一条由实线串起来的五角星形线x
subplot(2,1,2);
plot(n,y,'p-');%绘制一条由实线串起来的五角星形线y
实验分析:
这是设计的一个带阻滤波器,从上面的实验程序中可知,其fs=1000HZ,输入信号x=x1+x2,其中x1=2*sin(2*pi*4000*t),x2=5*sin(2*pi*50*t),是由两种不同的频率混合而成的,fs/2=500,而wp的变化范围为(0.2*pi)到(0.8*pi),由于带阻滤波器的幅度响应是关于π对称的,所以分析w从0到π可知,在π处的f=fs/2,根据比列可知fp1=100,fp2=400,其通带范围为0到100和400到500,所以x1、x2两种波形都能通过该滤波器;从上图中可以看出,没经滤波前的波形与滤波后的波形是一样的,即两种波形都没有被阻隔。
function[b,a]=afd_butt(Omegp,Omegs,Ap,As);
%b=Ha(s)分子多项式的系数
%a=Ha(s)分母多项式的系数
%Omegp=通带截止频率单位为rad/s;Ωp>0
%Omegs=阻带起始频率单位为rad/s;Ωs>Ωp>0
%Ap=Passbandripplein+dB;(Ap>0)
%As=Stopbandattenuationin+dB;(As>0)
ifOmegp<=0
error('Passbandedgemustbelargerthan0')
end
ifOmegs<=Omegp
error('StopbandedgemustbelargerthanPassbandedge')
end
if(Ap<=0)|(As<0)
error('PBrippleand/orSBattenuationmustbelargerthan0')
end
N=ceil((log10((10^(Ap/10)-1)/(10^(As/10)-1)))/(2*log10(Omegp/Omegs)));
fprintf('\n***ButterworthFilterOrder=%2.0f\n',N)
OmegaC=Omegp/((10^(Ap/10)-1)^(1/(2*N)));
[b,a]=u_buttap(N,OmegaC);
end
function[b,a]=u_buttap(N,Omegac);
%未归一化的Butterworth模拟低通滤波器原型
%[b,a]=u_buttap(N,Omegac);
%b=Ha(s)分子多项式的系数
%a=Ha(s)分母多项式的系数
%N=滤波器的阶数r
%Omegac=截止频率单位rad/s
[z,p,k]=buttap(N);
p=p*Omegac;
k=k*Omegac^N;
B=real(poly(z));
b=k*B;
b0=k;
a=real(poly(p));
end
5、实验小结
通过本次实验,我熟悉了用双线性变换法和用脉冲响应不变法设计IIR数字滤波器的原理与方法,也进一步加深了我对各种滤波器的原理及其工作的理解;此外,对于不懂的地方,不但要善于向别人请教,而且也要勤加思考。