数字信号处理实验Word格式文档下载.docx
《数字信号处理实验Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验Word格式文档下载.docx（13页珍藏版）》请在冰豆网上搜索。
subplot（4,3,4）;
plot（n,y1）;
直接滤波的结果'
subplot（4,3,5）;
plot（n,abs（y1）,'
g'
直接滤波幅频特性曲线'
subplot（4,3,6）;
plot（n,angle（y1）,'
m'
直接滤波相频特性曲线'
subplot（4,3,7）;
plot（n,y2,'
cascade滤波后的波形'
subplot（4,3,8）;
plot（n,abs（y2）,'
r'
cascade滤波幅频特性曲线'
subplot（4,3,9）;
plot（n,angle（y2）,'
b'
cascade滤波相频特性曲线'
subplot（4,3,10）;
plot（n,y3,'
parallel滤波的波形'
subplot（4,3,11）;
plot（n,abs（y3）,'
parallel滤波幅频特性曲线'
subplot（4,3,12）;
plot（n,angle（y3）,'
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];
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:
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;
K
w（i+1,:
）=filter（B（i,:
）,A（i,:
）,w（i,:
））;
hcas=b0*w（K+1,:
function[C,B,A]=dir2par（b,a）
[r1,p1,C]=residuez（b,a）;
p=cplxpair（p1,10000000*eps）;
I=cplxcomp（p1,p）;
r=r1（I）;
B=zeros（K,2）;
fori=1:
N-2
Brow=r（i:
Arow=p（i:
[Brow,Arow]=residuez（Brow,Arow,[]）;
）=real（Brow）;
）=real（Arow）;
end
[Brow,Arow]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（Brow）0];
A（K,:
）=[real（Arow）0];
N-1
functionhpar=parfilter（C,B,A,delta）
）=filter（C,1,delta）;
）,delta）;
hpar=sum（w）;
程序运行的结果如下：
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）;
plot（n,x2,'
x2的图像'
直接滤波波形'
cascade滤波波形'
plot（n,y3）;
parallel滤波图像'
plot（n,abs（y3））;
parallel滤波幅频特性'
plot（n,angle（y3））;
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）;
幅度响应'
gridon;
axis（[0,0.8,0,1]）;
ylabel（'
|幅度|'
xlabel（'
以\pi为单位的频率'
subplot（2,2,3）;
plot（w/pi,db）;
幅度响应（dB）'
axis（[0,0.8,-80,0]）;
对数幅度/dB'
subplot（2,2,2）;
plot（w/pi,pha）;
相位响应'
axis（[0,0.8,-4,4]）;
相位'
subplot（2,2,4）;
plot（w/pi,grd）;
群延迟'
axis（[0,0.8,0,10]）;
样本'
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'
[Brow,Arow]=residues（Brow,Arow,[]）;
function[db,mag,pha,grd,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'
whole'
H=（H（1:
501））'
;
w=（w（1:
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
grd=grpdelay（b,a,w）;
双线性变换：
%双线性变换设计数字巴特沃斯低通滤波器
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）
axis（[0,0.8,-60,0]）;
%将模拟滤波器直接形式转换为级联形式的函数
function[C,B,A]=sdir2cas（b,a）
%计算增益系数C
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:
Na
Arow=p（n:
n+1,:
Arow=poly（Arow）;
A（fix（（n+1）/2）,:
elseifNa==1
A=[0real（poly（p））];
A=zeros（K+1,3）;
A（K+1,:
）=[0real（poly（p（Na）））];
%计算分子的二阶因子部分
z=cplxpair（roots（b））;
K=floor（Nb/2）;
ifNb==0
B=[00poly（z）];
elseifK*2==Nb
B=zeros（K,3）;
Nb
Brow=z（n:
Brow=poly（Brow）;
B（fix（（n+1）/2）,:
elseifNb==1
B=[0real（poly（z））];
B=zeros（K+1,3）;
B（K+1,:
）=[0real（poly（z（Nb）））];
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（p）;
增益系数：
disp（k）;
二阶节：
disp（real（sos））;
并联型：
[r,p,k]=residuez（a,b）;
留数：
disp（r'
disp（p'
常数：
disp（k'
升级会员 