数字信号处理实验.docx
《数字信号处理实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验.docx（13页珍藏版）》请在冰豆网上搜索。
数字信号处理实验
1.设三阶的DF的系统函数
，当输入信号
或
，用三种结构实现IIRDF的结构并给出输出信号
，判断该DF的类型。
当输入
时，
b=[1,0.5814,0.2114,0];a=[1,-0.3984,0.2475,-0.04322];
delta=impseq（0,-10,10）;hdir=filter（b,a,delta）;[b0,B,A]=dir2cas（b,a）;
hcas=casfilter（b0,B,A,delta）;
[C,B,A]=dir2par（b,a）;
hpar=parfilter（C,B,A,delta）;
a=0.5;w0=pi./2;N=40;n=0:
N-1;x1=cos（n.*w0）;
y1=filter（b,a,x1）;y2=casfilter（b0,B,A,x1）;y3=parfilter（C,B,A,x1）;
subplot（4,3,2）;plot（n,x1）;title（'x1的图像'）;
subplot（4,3,4）;plot（n,y1）;title（'直接滤波的结果'）;
subplot（4,3,5）;plot（n,abs（y1）,'g'）;title（'直接滤波幅频特性曲线'）;
subplot（4,3,6）;plot（n,angle（y1）,'m'）;title（'直接滤波相频特性曲线'）;
subplot（4,3,7）;plot（n,y2,'m'）;title（'cascade滤波后的波形'）;
subplot（4,3,8）;plot（n,abs（y2）,'r'）;title（'cascade滤波幅频特性曲线'）;
subplot（4,3,9）;plot（n,angle（y2）,'b'）;title（'cascade滤波相频特性曲线'）;
subplot（4,3,10）;plot（n,y3,'g'）;title（'parallel滤波的波形'）;
subplot（4,3,11）;plot（n,abs（y3）,'m'）;title（'parallel滤波幅频特性曲线'）;
subplot（4,3,12）;plot（n,angle（y3）,'g'）;title（'parallel滤波相频特性曲线'）
调用的功能函数：
%单位取样序列函数
function[x,n]=impseq（n0,n1,n2）
n=[n1:
n2];x=[（n-n0）==0];
%IIFDF直接型到级联型的转换函数
function[b0,B,A]=dir2cas（b,a）
%计算增益系数
b0=b
（1）;b=b/b0;
a0=a
（1）;a=a/a0;
b0=b0/a0;
%
M=length（b）;N=length（a）;
ifN>M
b=[b,zeros（1,N-M）];
elseifM>N
a=[a,zeros（1,M-N）];N=M;
else
NM=0;
end
%
K=floor（N/2）;B=zeros（K,3）;A=zeros（K,3）;
ifK*2==N
b=[b,0];a=[a,0];
end
%
broots=cplxpair（roots（b））;
aroots=cplxpair（roots（a））;
fori=1:
2:
2*K
Brow=broots（i:
1:
i+1,:
）;
Brow=real（poly（Brow））;
B（fix（（i+1）/2）,:
）=Brow;
Arow=aroots（i:
1:
i+1,:
）;
Arow=real（poly（Arow））;
A（fix（（i+1）/2）,:
）=Arow;
End
functionhcas=casfilter（b0,B,A,delta）
[K,L]=size（B）;N=length（delta）;
w=zeros（K+1,N）;
w（1,:
）=delta;
fori=1:
1:
K
w（i+1,:
）=filter（B（i,:
）,A（i,:
）,w（i,:
））;
end
hcas=b0*w（K+1,:
）;
function[C,B,A]=dir2par（b,a）
M=length（b）;N=length（a）;
[r1,p1,C]=residuez（b,a）;
p=cplxpair（p1,10000000*eps）;
I=cplxcomp（p1,p）;r=r1（I）;
K=floor（N/2）;B=zeros（K,2）;A=zeros（K,3）;
ifK*2==N
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）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
[Brow,Arow]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（Brow）0];A（K,:
）=[real（Arow）0];
else
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）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
end
functionhpar=parfilter（C,B,A,delta）
[K,L]=size（B）;N=length（delta）;
w=zeros（K+1,N）;
w（1,:
）=filter（C,1,delta）;
fori=1:
1:
K
w（i+1,:
）=filter（B（i,:
）,A（i,:
）,delta）;
end
hpar=sum（w）;
程序运行的结果如下：
当输入
时，
b=[1,0.5814,0.2114,0];a=[1,-0.3984,0.2475,-0.04322];
delta=impseq（0,-10,10）;hdir=filter（b,a,delta）;[b0,B,A]=dir2cas（b,a）;
hcas=casfilter（b0,B,A,delta）;
[C,B,A]=dir2par（b,a）;hpar=parfilter（C,B,A,delta）;
n=[-20,20];[x2,n]=impseq（0,-20,20）;
y1=filter（b,a,x2）;y2=casfilter（b0,B,A,x2）;y3=parfilter（C,B,A,x2）;
subplot（4,3,2）;plot（n,x2,'b'）;title（'x2的图像'）;
subplot（4,3,4）;plot（n,y1）;title（'直接滤波波形'）;
subplot（4,3,5）;plot（n,abs（y1）,'b'）;title（'直接滤波幅频特性曲线'）;
subplot（4,3,6）;plot（n,angle（y1）,'r'）;title（'直接滤波相频特性曲线'）;
subplot（4,3,7）;plot（n,y2,'g'）;title（'cascade滤波波形'）;
subplot（4,3,8）;plot（n,abs（y2）,'m'）;title（'cascade滤波幅频特性曲线'）;
subplot（4,3,9）;plot（n,angle（y2）,'g'）;title（'cascade滤波相频特性曲线'）;
subplot（4,3,10）;plot（n,y3）;title（'parallel滤波图像'）;
subplot（4,3,11）;plot（n,abs（y3））;title（'parallel滤波幅频特性'）;
subplot（4,3,12）;plot（n,angle（y3））;title（'parallel滤波相频特性'）;
程序运行的结果如下：
2.用Matlab设计一个模拟巴特沃斯低通滤波器，指标
，采用冲激响应不变法和双线性变换设计。
冲激响应不变法：
%冲激响应不变法设计数字巴特沃斯低通滤波器
%数字滤波器指标
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
%转换为模拟域指标
T=1;
omegap=wp/T;omegas=ws/T;
%模拟巴特沃斯滤波器的计算
[cs,ds]=afd_buttap（omegap,omegas,Rp,As）;
%冲激不变法
[b,a]=imp_invr（cs,ds,T）;
[C,B,A]=dir2par（b,a）
%计算数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度响应'）;gridon;
axis（[0,0.8,0,1]）;ylabel（'|幅度|'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;gridon;
axis（[0,0.8,-80,0]）;
ylabel（'对数幅度/dB'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,2）;plot（w/pi,pha）;title（'相位响应'）;gridon;
axis（[0,0.8,-4,4]）;ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,4）;plot（w/pi,grd）;title（'群延迟'）;gridon;
axis（[0,0.8,0,10]）;xlabel（'以\pi为单位的频率'）;ylabel（'样本'）;
调用的功能函数：
function[b,a]=afd_buttap（wp,ws,Rp,As）
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（wp/ws）））;
fprintf（'\n***ButterworthFilterOrder=%2.0f\n',N）
omega=wp/（（10^（Rp/10）-1）^（1/（2*N）））;
[b,a]=buttap_o（N,omega）;
function[b,a]=imp_invr（c,d,T）
[R,p,k]=residue（c,d）;
p=exp（p*T）;
[b,a]=residuez（R,p,k）;
b=real（b'）;a=real（a'）;
function[C,B,A]=dir2par（b,a）
M=length（b）;N=length（a）;
[r1,p1,C]=residuez（b,a）;
p=cplxpair（p1,10000000*eps）;
I=cplxcomp（p1,p）;r=r1（I）;
K=floor（N/2）;B=zeros（K,2）;A=zeros（K,3）;
ifK*2==N
fori=1:
2:
N-2
Brow=r（i:
1:
i+1,:
）;
Arow=p（i:
1:
i+1,:
）;
[Brow,Arow]=residues（Brow,Arow,[]）;
B（fix（（i+1）/2）,:
）=real（Brow）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
[Brow,Arow]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（Brow）0];A（K,:
）=[real（Arow）0];
else
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）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
end
function[db,mag,pha,grd,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
1:
501））';w=（w（1:
1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;grd=grpdelay（b,a,w）;
程序运行的结果如下：
双线性变换：
%双线性变换设计数字巴特沃斯低通滤波器
%数字滤波器指标
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
%转换为模拟域指标
T=1;Fs=1/T;
omegap=（2/T）*tan（wp/2）;omegas=（2/T）*tan（ws/2）;
%模拟巴特沃斯滤波器的计算
[cs,ds]=afd_buttap（omegap,omegas,Rp,As）
%双线性变换
[b,a]=bilinear（cs,ds,Fs）;
[C,B,A]=sdir2cas（b,a）
%计算数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度响应'）;gridon;
ylabel（'|幅度|'）;xlabel（'以\pi为单位的频率'）;
axis（[0,0.8,0,1]）;
subplot（2,2,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;gridon;
ylabel（'对数幅度/dB'）;xlabel（'以\pi为单位的频率'）;
axis（[0,0.8,-60,0]）;
subplot（2,2,2）;plot（w/pi,pha）;title（'相位响应'）;gridon;
axis（[0,0.8,-4,4]）;ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,4）;plot（w/pi,grd）;title（'群延迟'）;gridon;
axis（[0,0.8,0,10]）;xlabel（'以\pi为单位的频率'）;ylabel（'样本'）;
调用的功能函数：
function[b,a]=afd_buttap（wp,ws,Rp,As）
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（wp/ws）））;
fprintf（'\n***ButterworthFilterOrder=%2.0f\n',N）
omega=wp/（（10^（Rp/10）-1）^（1/（2*N）））;
[b,a]=buttap_o（N,omega）;
%将模拟滤波器直接形式转换为级联形式的函数
function[C,B,A]=sdir2cas（b,a）
%计算增益系数C
b0=b
（1）;b=b/b0;a0=a
（1）;a=a/a0;C=b0/a0
Na=length（a）-1;Nb=length（b）-1;
%计算分母的二阶因子部分
p=cplxpair（roots（a））;K=floor（Na/2）;
ifK*2==Na
A=zeros（K,3）;
forn=1:
2:
Na
Arow=p（n:
1:
n+1,:
）;
Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
elseifNa==1
A=[0real（poly（p））];
else
A=zeros（K+1,3）;
forn=1:
2:
2*K
Arow=p（n:
1:
n+1,:
）;
Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
A（K+1,:
）=[0real（poly（p（Na）））];
end
%计算分子的二阶因子部分
z=cplxpair（roots（b））;
K=floor（Nb/2）;
ifNb==0
B=[00poly（z）];
elseifK*2==Nb
B=zeros（K,3）;
forn=1:
2:
Nb
Brow=z（n:
1:
n+1,:
）;
Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
elseifNb==1
B=[0real（poly（z））];
else
B=zeros（K+1,3）;
forn=1:
2:
2*K
Brow=z（n:
1:
n+1,:
）;
Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
B（K+1,:
）=[0real（poly（z（Nb）））];
end
function[db,mag,pha,grd,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
1:
501））';w=（w（1:
1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;grd=grpdelay（b,a,w）;
程序运行的结果如下：
3.用级联型和并联型结构实现以下系统函数，每个二阶节都采用直接II型结构。
级联型：
a=[5,-2.4412*5,2.4412*5,-5];%输入分子系数向量
b=[1,-1.7728,1.4464,-0.405];%输入分母系数向量
[z,p,k]=tf2zp（a,b）;%求出各子系统的零极点
[sos,A]=zp2sos（z,p,k）;%求出各二阶节乘系数
disp（'零点：
'）;disp（z）;
disp（'极点：
'）;disp（p）;
disp（'增益系数：
'）;disp（k）;
disp（'二阶节：
'）;disp（real（sos））;
程序运行的结果如下：
并联型：
a=[5,-2.4412*5,2.4412*5,-5];%输入分子系数向量
b=[1,-1.7728,1.4464,-0.405];%输入分母系数向量
[r,p,k]=residuez（a,b）;%求出各子系统的零极点
disp（'留数：
'）;disp（r'）;
disp（'极点：
'）;disp（p'）;
disp（'常数：
'）;disp（k'）;
程序运行的结果如下：
升级会员 