数字信号处理实验报告.docx
《数字信号处理实验报告.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验报告.docx(16页珍藏版)》请在冰豆网上搜索。
数字信号处理实验报告
郑州航空工业管理学院
《数字信号处理》实验报告
专业电子信息工程
学号
姓名
实验一数字滤波器的结构
一、实验目的
(1)加深对数字滤波器分类与结构的了解;
(2)明确数字滤波器的基本结构及其相互间的转换方法;
(3)掌握用MATLAB进行数字滤波器各种结构相互间转换的子函数及程序编写方法。
二、实验原理
一个离散LSI系统可用系统函数来表示;
也可用差分方程来表示:
当
至少有一个不为0时,则在有限z平面上存在极点,表示一个IIR数字滤波器;当
全都为0时,系统不存在极点,表示一个FIR系统。
IIR数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、级联型和并联型。
FIR数字滤波器的基本结构分为横截型、级联型、并联型、、线性相位型和频率抽样型。
三、实验仪器
微型计算机、MATLAB
四、实验内容
(1)已知一个IIR系统的系统函数为
将其从直接型转换为级联型和并联型结构,并画出各种结构的流程图。
(2)已知一个FIR系统的系统函数为
将其从横截型转换为级联型结构,并画出各种结构的流程图。
五、实验程序
num=[0,10,0];
den=[1,-3,2];
[z,p,k]=tf2zp(num,den)
[b,a]=zp2tf(z,p,k)
num=[1.9,2.5,2.5,1.9];
den=[1,-6,5,-0.4];
[sos,g]=tf2sos(num,den)
[b,a]=sos2tf(sos,g)
z=[1,-3]';
p=[2,-4]';
k=5;
[sos,g]=zp2sos(z,p,k)
[z1,p1,k1]=sos2zp(sos,g)
%直接型到并联型的转换
function[C,B,A]=dir2par(num,den)
M=length(num);N=length(den);
[r1,p1,C]=residuez(num,den);%先求系统的单根p1,对应的留数r1及直接项C
p=cplxpair(p1,10000000*eps);%用配对函数cplxpair由p1找共轭复根p
I=cplxcomp(p1,p);%找p1变为p时的排序变化
r=r1(I);%让r1的排序变化为r,保持与极点对应
%变换为二阶子系统
K=floor(N/2);B=zeros(K,2);A=zeros(K,3);%二阶子系统变量的初始化
ifK*2==N;%N为偶数,A(z)的次数为奇,有一个因式是一阶的
fori=1:
2:
N-2
Brow=r(i:
1:
i+1,:
);%取出一对留数
Arow=p(i:
1:
i+1,:
);%取出一对对应的极点
%二个留数极点转为二阶子系统分子分母系数
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:
)=real(Brow);%取Brow的实部,放入系数矩阵B的相应行
A(fix((i+1)/2),:
)=real(Arow);%取Arow的实部,放入系数矩阵A的相应行
end;
[Brow,Arow]=residuez(r(N-1),p(N-1),[]);%处理实单根
B(K,:
)=[real(Brow),0];A(K,:
)=[real(Arow),0];
else%N为奇数,A(z)的次数为偶,所有因式都是二阶的
fori=1:
2:
N-1
Brow=r(i:
1:
i+1,:
);%取出一对留数
Arow=p(i:
1:
i+1,:
);%取出一对对应的极点
%二个留数极点转为二阶子系统分子分母系数
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:
)=real(Brow);%取Brow的实部,放入系数矩阵B的相应行
A(fix((i+1)/2),:
)=real(Arow);%取Arow的实部,放入系数矩阵A的相应行
end
end
num=[8-411-2];
den=[1-1.250.75-0.125];
[C,B,A]=dir2par(num,den)
C=
16
B=
-16.000020.0000
8.00000
A=
1.0000-1.00000.5000
1.0000-0.25000
试验结果分析
实验二用冲激响应不变法设计IIR数字滤波器
一、实验目的
(1)加深对冲激响应不变法设计IIR数字滤波器的基本原理的理解;
(2)掌握用冲激响应不变法设计数字低通、带通滤波器的设计;
(3)了解MATLAB有关冲激响应不变法的常用子函数。
二、实验原理
1.设计思想
是将系统从s平面到z平面的一种映射方法,使数字滤波器的单位冲激响应序列
模仿模拟滤波器的单位冲激响应
,其变换关系式为
。
2.设计步骤
(
(1)确定数字滤波器的性能指标:
通带临界频率
、阻带临界频率
;通带内的最大衰减
;阻带内的最小衰减
;
(2)确定相应的数字角频率
;
(3)将数字滤波器的数字指标转换为模拟滤波器的指标
;
(4)根据
计算模拟低通原型滤波器的阶数N,并求得低通原型的传递函数
;
(5)用冲激响应不变法将模拟滤波器转换为数字滤波器
;
(6)分析滤波器特性,检查其指标是否满足要求。
三、实验仪器
微型计算机、MATLAB
四、实验内容
采用冲激响应不变法设计一个巴特沃斯数字低通滤波器,要求:
通带截止频率
,阻带截止频率
,滤波器的采样频率
,并画出滤波器的幅频响应曲线和相频响应曲线以及零极点图。
五、实验程序
%脉冲响应不变法设计巴特沃斯低通滤波器
%数字滤波器指标
wp=0.25*pi;%滤波器的通带截止频率
ws=0.4*pi;%滤波器的阻带截止频率
Rp=1;As=15;%输入滤波器的通阻带衰减指标
ripple=10^(-Rp/20);%计算通带衰减对应的幅度值
Attn=10^(-As/20);%计算阻带衰减对应的幅度值
%转换为模拟滤波器指标
Fs=2000;T=1/Fs;
Omgp=wp*Fs;Omgs=ws*Fs;
%模拟原型滤波器计算
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s')%计算阶数n和截止频率
[z0,p0,k0]=buttap(n);%设计归一化的模拟原型滤波器
ba1=k0*real(poly(z0));%求原型滤波器系数b
aa1=real(poly(p0));%求原型滤波器系数a
[ba,aa]=lp2lp(ba1,aa1,Omgc);%变换为模拟低通滤波器
%用脉冲响应不变法计算数字滤波器系数
[bd,ad]=impinvar(ba,aa,Fs)
[C,B,A]=dir2par(bd,ad)%转换成并联型
%求数字系统的频率特性
[H,w]=freqz(bd,ad);
dbH=20*log10(abs(H)/max(abs(H)));%化为分贝值
%
subplot(2,2,1),plot(w/pi,abs(H));
ylabel('幅度');xlabel('频率/\pi');axis([0,1,0,1.1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);grid
subplot(2,2,2),plot(w/pi,angle(H)/pi);
ylabel('相位');xlabel('频率/\pi');axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-1,0,1]);grid
subplot(2,2,3),plot(w/pi,dbH);
ylabel('幅度(dB)');xlabel('频率/\pi');axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);
set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);grid
subplot(2,2,4),zplane(bd,ad);
axis([-1.1,1.1,-1.1,1.1]);
ylabel('零极图');
function[C,B,A]=dir2par(num,den)
%直接型到并联型的转换
M=length(num);N=length(den);
[r1,p1,C]=residuez(num,den);%先求系统的单根p1,对应的留数r1及直接项C
p=cplxpair(p1,10000000*eps);%用配对函数cplxpair由p1找共轭复根p
I=cplxcomp(p1,p);%找p1变为p时的排序变化
r=r1(I);%让r1的排序变化为r,保持与极点对应
%变换为二阶子系统
K=floor(N/2);B=zeros(K,2);A=zeros(K,3);%二阶子系统变量的初始化
ifK*2==N;%N为偶数,A(z)的次数为奇,有一个因式是一阶的
fori=1:
2:
N-2
Brow=r(i:
1:
i+1,:
);%取出一对留数
Arow=p(i:
1:
i+1,:
);%取出一对对应的极点
%二个留数极点转为二阶子系统分子分母系数
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:
)=real(Brow);%取Brow的实部,放入系数矩阵B的相应行
A(fix((i+1)/2),:
)=real(Arow);%取Arow的实部,放入系数矩阵A的相应行
end;
[Brow,Arow]=residuez(r(N-1),p(N-1),[]);%处理实单根
B(K,:
)=[real(Brow),0];A(K,:
)=[real(Arow),0];
else%N为奇数,A(z)的次数为偶,所有因式都是二阶的
fori=1:
2:
N-1
Brow=r(i:
1:
i+1,:
);%取出一对留数
Arow=p(i:
1:
i+1,:
);%取出一对对应的极点
%二个留数极点转为二阶子系统分子分母系数
[Brow,Arow]=residuez(Brow,Arow,[]);
B(fix((i+1)/2),:
)=real(Brow);%取Brow的实部,放入系数矩阵B的相应行
A(fix((i+1)/2),:
)=real(Arow);%取Arow的实部,放入系数矩阵A的相应行
end
end
六、试验结果分析
七、思考题
使用冲激响应不变法设计IIR数字滤波器时,有哪些缺点?
该方法在设计数字滤波器时,有何限制?
实验三用双线性变换法设计IIR数字滤波器
一、实验目的
(1)加深对双线性变换法设计FIR数字滤波器的基本原理的理解;
(2)掌握用双线性变换法设计数字低通、高通、带通滤波器的设计;
(3)了解MATLAB有关双线性变换法设计的常用子函数。
二、实验原理
1.设计思想
s平面与z平面之间满足以下映射关系:
s平面的虚轴单值地映射于z平面的单位圆上,s平面的左半平面完全映射到z平面的单位圆内。
双线性变换不存在混叠问题。
双线性变换时一种非线性变换
,这种非线性引起的幅频特性畸变可通过预畸而得到校正。
2.设计步骤
(1)确定数字滤波器的性能指标:
通带临界频率
、阻带临界频率
;通带内的最大衰减
;阻带内的最小衰减
;
(2)确定相应的数字角频率
;
(3)计算经过预畸的相应模拟低通原型的频率
;
(4)根据
计算模拟低通原型滤波器的阶数N,并求得低通原型的传递函数
;
(5)用上面的双线性变换公式代入
,求出所设计的传递函数
;
(6)分析滤波器特性,检查其指标是否满足要求。
三、实验仪器
微型计算机、MATLAB
四、实验内容
(1)采用双线性变换法设计一个巴特沃斯数字低通滤波器,要求:
通带截止频率
,阻带截止频率
,滤波器的采样频率
,并画出滤波器的幅频响应曲线和相频响应曲线以及零极点图。
五、实验程序
%双线性变换法设计数字低通滤波器
%数字滤波器指标
fp=150;fs=250;Fs=800;T=1/Fs;
wp=fp/Fs*2*pi;%数字滤波器的通带截止频率
ws=fs/Fs*2*pi;%数字滤波器的阻带截止频率
Rp=3;As=20;%输入滤波器的通阻带衰减指标
%转换为模拟滤波器指标
Omgp=(2/T)*tan(wp/2);
Omgs=(2/T)*tan(ws/2);
[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s')%计算阶数n和截止频率
%方法1:
模拟原型滤波器计算
[z0,p0,k0]=buttap(n);%归一化原型设计
ba=k0*real(poly(z0));%求原型滤波器系数b
aa=real(poly(p0));%求原型滤波器系数a
[ba1,aa1]=lp2lp(ba,aa,Omgc);%变换为模拟低通滤波器
[bd1,ad1]=bilinear(ba1,aa1,Fs)%双线性变换
[H1,w1]=freqz(bd1,ad1);
dbH1=20*log10(abs(H1)/max(abs(H1)));%化为分贝值
%方法2:
直接求模拟滤波器系数
[ba2,aa2]=butter(n,Omgc,'s');
%用双线性变换法计算数字滤波器系数
[bd2,ad2]=bilinear(ba2,aa2,Fs)%双线性变换
[H2,w2]=freqz(bd2,ad2);
dbH2=20*log10(abs(H2)/max(abs(H2)));%化为分贝值
%方法3:
直接求数字滤波器系数
[n3,wc3]=buttord(wp/pi,ws/pi,Rp,As)%计算数字滤波器n和wc
[bd3,ad3]=butter(n3,wc3)
[H3,w3]=freqz(bd3,ad3);
dbH3=20*log10(abs(H3)/max(abs(H3)));%化为分贝值
%
subplot(3,2,1),plot(w1/2/pi*Fs,dbH1,'k');
title('方法1幅度响应(dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,2),plot(w1/2/pi*Fs,angle(H1)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(3,2,3),plot(w2/2/pi*Fs,dbH2,'k');
title('方法2幅度响应(dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,4),plot(w2/2/pi*Fs,angle(H2)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
subplot(3,2,5),plot(w3/2/pi*Fs,dbH3,'k');
title('直接法幅度响应(dB)');axis([0,Fs/2,-40,5]);
ylabel('dB');xlabel('频率(Hz)');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(3,2,6),plot(w3/2/pi*Fs,angle(H3)/pi*180,'k');
title('相位响应');axis([0,Fs/2,-180,180]);
ylabel('\phi');xlabel('频率(Hz)');
set(gca,'XTickMode','manual','XTick',[0,fp,fs,Fs/2]);
set(gca,'YTickMode','manual','YTick',[-180,0,180]);grid
六、试验结果分析
七、思考题
使用双线性变换法时,数字频率与模拟频率有何关系?
会带来什么影响?
如何解决?
实验四用窗函数法设计FIR数字滤波器
一、实验目的
(4)加深对窗函数法设计FIR数字滤波器的基本原理的理解;
(5)学习用MATLAB语言的窗函数法编写设计FIR数字滤波器的程序;
(6)了解MATLAB有关窗函数法设计的常用子函数。
二、实验原理
1.设计思想
从时域从发,设计
逼近理想
。
设理想滤波器
的单位脉冲响应为
。
一般是无限长的,且是非因果的,不能直接作为FIR滤波器的单位脉冲响应。
要想得到一个因果的有限长的滤波器h(n),最直接的方法是截断
,即截取为有限长因果序列,并用合适的窗函数进行加权作为FIR滤波器的单位脉冲响应。
2.设计步骤
(1)根据过渡带和阻带衰减设计指标选择窗函数的类型,并估算滤波器的阶数
:
(2)由数字滤波器的理想频率响应
求出理想的单位冲激响应
。
(3)计算数字滤波器的单位冲激响应
。
(4)检查设计出的滤波器是否满足技术指标。
若不满足技术指标,则需重新选择或调整窗函数的类型,估算滤波器的阶数
,再重复前面的四个步骤,直到满足指标为止。
三、实验仪器
微型计算机、MATLAB
四、实验内容
选择合适的窗函数设计FIR数字低通滤波器。
要求:
通带截止频率
,阻带截止频率
,描绘实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线和相频响应曲线。
五、实验程序
functionhd=ideal_lp(wc,N)
%hd=点0到N-1之间的理想脉冲响应
%wc=截止频率(弧度)
%N=理想滤波器的长度
tao=(N-1)/2;
n=[0:
(N-1)];
m=n-tao+eps;%加一个小数以避免0作除数
hd=sin(wc*m)./(pi*m);
%freqz_m.m
function[db,mag,pha,grd,w]=freqz_m(b,a);
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:
501))';w=(w(1:
501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
grd=grpdelay(b,a,w);
%FIR低通滤波器
wp=0.3*pi;ws=0.45*pi;%输入设计指标
deltaw=ws-wp;%计算过渡带的宽度
N0=ceil(6.6*pi/deltaw)%按哈明窗计算滤波器长度N0
N=N0+mod(N0+1,2)%为实现FIR类型I偶对称滤波器,应确保N为奇数
windows=(hamming(N))';%使用哈明窗,并将列向量变为行向量
wc=(ws+wp)/2;%截止频率取通阻带频率的平均值
hd=ideal_lp(wc,N);%建立理想低通滤波器
b=hd.*windows%求FIR系统函数系数
[db,mag,pha,grd,w]=freqz_m(b,1);%求解频率特性
n=0:
N-1;dw=2*pi/1000;%dw为频率分辨率,将0~2π分为1000份
Rp=-(min(db(1:
wp/dw+1)))%检验通带波动
As=-round(max(db(ws/dw+1:
501)))%检验最小阻带衰减
subplot(2,2,1),stem(n,b,'filled','k');
axis([0,N,1.1*min(b),1.1*max(b)]);title('实际脉冲响应');
xlabel('n');ylabel('h(n)');
subplot(2,2,2),stem(n,windows,'filled','k');
axis([0,N,0,1.1]);title('窗函数特性');
xlabel('n');ylabel('wd(n)');
subplot(2,2,3),plot(w/pi,db,'k');
axis([0,1,-80,10]);title('幅频响应');
xlabel('频率(单位:
\pi)');ylabel('H(e^{j\omega})');
set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[-50,-20,-3,0]);grid
subplot(2,2,4),plot(w/pi,pha,'k');
axis([0,1,-4,4]);
title('相频响应');
xlabel('频率(单位:
\pi)');ylabel('\phi(\omega)');
set(gca,'XTickMode','manual','XTick',[0,wp/pi,ws/pi,1]);
set(gca,'YTickMode','manual','YTick',[-pi,0,pi]);grid
%FIR低通滤波器
wp=0.3*pi;ws=0.45*pi;%输入设计指标
deltaw=ws-wp;%计算过渡带的宽度
N0=ceil(6.6*pi/deltaw);%按哈明窗计算滤波器长度N0
N=N0+mo