matlab上机实验 数字信号处理实验课必备.docx
《matlab上机实验 数字信号处理实验课必备.docx》由会员分享，可在线阅读，更多相关《matlab上机实验 数字信号处理实验课必备.docx（31页珍藏版）》请在冰点文库上搜索。
matlab上机实验数字信号处理实验课必备
MATLAB上机实验
上机实验一
主要内容：
系统频率响应；系统函数零、极点图；系统的单位脉冲响应。
10-1
，求系统频率响应；系统函数零、极点图；系统的单位脉冲响应。
b=[0.20.31];%分子系数
a=[10.41];%分母系数
figure
（1）,zplane（b,a）;%零极点图
%%求频率响应
b=[0.20.31];%分子系数
a=[10.41];%分母系数
figure
（2）,freqz（b,a）;%数字系统频响
%%求系统的单位脉冲响应
b=[0.20.31];%分子系数
a=[10.41];%分母系数
figure（3）,impz（b,a）;%系统的单位脉冲响应
10-2已知
分析其频率特性
，并作
、
图。
w=[0:
1:
500]*2*pi/500;%[0，2pi]区域分为501点
X1=1;
X2=1-0.9.*exp（-1*j*w）;
X=X1./X2;
magX=abs（X）;
angX=angle（X）.*180./pi;
subplot（211）;plot（w/pi,magX）;
xlabel（'以pi为单位的频率'）;
title（'幅度部分'）;
ylabel（'幅度'）;
subplot（212）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;
ylabel（'相位'）;
10-3已知
，画出其幅频特性。
w=[0:
1:
500]*2*pi/500;%[0,pi]区域分为501点
X1=1-exp（-6*j*w）;%分子多项式
magX=abs（X1）;
angX=angle（X1）.*180/pi;
subplot（211）;plot（w/pi,magX）;
title（'幅度部分'）;
ylabel（'振幅'）;
xlabel（'单位：
pi'）;
subplot（212）;
plot（w/pi,angX）;
xlabel（'单位：
pi'）;
title（'相位部分'）;
ylabel（'相位'）;
10-4分析矩形序列的频率特性，并作出幅频特性图。
N=6
w=[0:
1:
500]*pi/500;%[0，pi]区域分为501点
X1=1-exp（-1*j*w*N）;
X2=1-exp（-1*j*w）;
X22=X2+（X2==0）*eps;%逻辑数组参加运算，使“0”被“机器零“代替
X=X1./X22;
magX=abs（X）;
angX=angle（X）.*180./pi;
figure
（1）,subplot（211）;plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'幅度'）;
subplot（212）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
%%N=9
N=9;
w=[0:
1:
500]*pi/500;%[0，pi]区域分为501点
X1=1-exp（-1*j*w*N）;
X2=1-exp（-1*j*w）;
X22=X2+（X2==0）*eps;%逻辑数组参加运算，使“0”被“机器零“代替
X=X1./X22;
magX=abs（X）;
angX=angle（X）.*180./pi;
figure
（2）,subplot（211）;plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'幅度'）;
subplot（212）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
10-5求横向结构网络
的频响图。
w=[0:
1:
500]*.2*pi/500;%[0,2pi]区域分为501点
X1=1-0.9^8.*exp（-8*j*w）;
X2=1-0.9.*exp（-1*j*w）;
X=X1./X2;
magX=abs（X）;
angX=angle（X）.*180./pi;
subplot（211）;plot（w/pi,magX）;
xlabel（'以pi为单位的频率'）;
title（'幅度部分'）;
ylabel（'幅度'）;
subplot（212）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;
ylabel（'相位'）;
10-6已知系统函数
；求零、极点图；单位脉冲响应；系统频响。
b=[1000001];
a=[1000000];
w=[0:
1:
500]*pi/500;
X1=1+exp（-6*j*w）;
X2=1%分母多项式
X=X1./X2;
magX=abs（X）;
angX=angle（X）.*180./pi;
figure
（1）,subplot（221）;impz（b,a,15）;%系统的单]位脉冲响应
title（'冲激响应'）;
subplot（223）;
zplane（b,a）;%零极点图
title（'零极点图'）;
subplot（222）;
plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'振幅'）;
subplot（224）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
10-7已知系统函数
；求零、极点图；单位脉冲响应；系统频响。
b=[1-0.50];
a=[1-11];
w=[0:
1:
500]*pi/500;%[0，pi]区域分为501点
x1=1-0.5*exp（-1*j*w）;%分子多项式
x2=1-exp（-1*j*w）+exp（-2*j*w）;%分母多项式
x22=x2+（x2==0）*eps;%逻辑数组参加运算，使“0”被“机器零“代替
x=x1./x22;
magX=abs（x）;
angX=angle（x）.*180./pi;
figure
（1）,subplot（221）;impz（b,a,27）;%系统的单]位脉冲响应
title（'冲激响应'）;
subplot（223）;
zplane（b,a）;%零极点图
title（'零极点图'）;
subplot（222）;
plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'振幅'）;
subplot（224）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
10-8已知系统函数
；求零、极点图；单位脉冲响应；系统频响。
b=[1-0.500001-0.5];
a=[1-1100000];
w=[0:
1:
500]*pi/500;%[0，pi]区域分为501点
x1=1-0.5*exp（-1*j*w）+exp（-6*j*w）-0.5*exp（-7*j*w）;
x2=1-exp（-1*j*w）+exp（-2*j*w）;
x22=x2+（x2==0）*eps;%逻辑数组参加运算，使“0”被“机器零“代替
x=x1./x22;
magX=abs（x）;
angX=angle（x）.*180./pi;
figure
（1）,subplot（221）;impz（b,a,27）;%系统的单]位脉冲响应
title（'冲激响应'）;
subplot（223）;
zplane（b,a）;%零极点图
title（'零极点图'）;
subplot（222）;
plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'振幅'）;
subplot（224）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
10-9已知系统函数
；求零、极点图；单位脉冲响应；系统频响。
b=[10.5-0.5-1-0.51];
a=[100000];
w=[0:
1:
500]*pi/500;%[0，pi]区域分为501点
x1=1+0.5*exp（-1*j*w）-exp（-3*j*w）-0.5*exp（-2*j*w）-0.5*exp（-4*j*w）+exp（-5*j*w）;
x2=1;
x=x1./x2;
magX=abs（x）;
angX=angle（x）.*180./pi;
figure
（1）,subplot（221）;impz（b,a,15）;%系统的单]位脉冲响应
title（'冲激响应'）;
subplot（223）;
zplane（b,a）;%零极点图
title（'零极点图'）;
subplot（222）;
plot（w/pi,magX）;
title（'幅度部分'）;ylabel（'振幅'）;
subplot（224）;plot（w/pi,angX）;
xlabel（'以pi为单位的频率'）;
title（'相位部分'）;ylabel（'相位'）;
上机实验二
主要内容：
离散傅里叶级数；离散傅里叶变换；数字滤波器结构。
10-10离散傅里叶级数
，绘制下列情况下的
幅度
（1）L=5，N=20；
（2）L=5，N=40；（3）L=5，N=60；（4）L=7，N=60
L=5;N=20;
n=1:
N;
xn=[ones（1,L）,zeros（1,N-L）];
Xk=dfs（xn,N）;
magXk=abs（[Xk（N/2+1:
N）Xk（1:
N/2+1）]）;
k=[-N/2:
N/2];
figure
（1）
subplot（211）
stem（n,xn）;
xlabel（'n'）;ylabel（'xtilde（n）'）;
title（'DFSofSQ.wave:
L=50,N=20'）;
subplot（212）
stem（k,magXk）;xlabel（'n'）;
axis（[-N/2,N/2,-0.5,5.5]）;
xlabel（'k'）;ylabel（'xtilde（k）'）;
%part2
L=5;N=40;
n=1:
N;
xn=[ones（1,L）,zeros（1,N-L）];
Xk=dfs（xn,N）;
magXk=abs（[Xk（N/2+1:
N）Xk（1:
N/2+1）]）;
k=[-N/2:
N/2];
figure
（2）
subplot（211）
stem（n,xn）;
xlabel（'n'）;ylabel（'xtilde（n）'）;
title（'DFSofSQ.wave:
L=5,N=40'）;
subplot（212）
stem（k,magXk）;xlabel（'n'）;
axis（[-N/2,N/2,-0.5,5.5]）;
xlabel（'k'）;ylabel（'xtilde（k）'）;
%part3
L=5;N=60;
n=1:
N;
xn=[ones（1,L）,zeros（1,N-L）];
Xk=dfs（xn,N）;
magXk=abs（[Xk（N/2+1:
N）Xk（1:
N/2+1）]）;
k=[-N/2:
N/2];
figure（3）
subplot（211）
stem（n,xn）;
xlabel（'n'）;ylabel（'xtilde（n）'）;
title（'DFSofSQ.wave:
L=5,N=60'）;
subplot（212）
stem（k,magXk）;xlabel（'n'）;
axis（[-N/2,N/2,-0.5,5.5]）;
xlabel（'k'）;ylabel（'xtilde（k）'）;
%part4
L=7;N=60
n=1:
N;
xn=[ones（1,L）,zeros（1,N-L）];
Xk=dfs（xn,N）;
magXk=abs（[Xk（N/2+1:
N）Xk（1:
N/2+1）]）;
k=[-N/2:
N/2];
figure（4）
subplot（211）
stem（n,xn）;
xlabel（'n'）;ylabel（'xtilde（n）'）;
title（'DFSofSQ.wave:
L=7,N=60'）;
subplot（212）
stem（k,magXk）;xlabel（'n'）;
axis（[-N/2,N/2,-0.5,7.5]）;
xlabel（'k'）;ylabel（'xtilde（k）'）;
扩展函数dfs
function[Xk]=dfs（xn,N）
n=[0:
1:
N-1];%n的行矢量
k=[0:
1:
N-1];
WN=exp（-j*2*pi/N）;%WN因子
nk=n'*k;%得到nk的N*N的矩阵值
WNnk=WN.^nk;%DFS矩阵
Xk=xn*WNnk;%DFS系数的行矢量
10-11离散傅里叶
（1）取x（n）的前10点数据，求N=10点的X（k）。
（2）将a中的x（n）补零至100点，求N=100点的X（k）。
（3）取x（n）的前100点数据，求N=100点的X（k）。
figure
（1）;
n1=[0:
1:
9];
y1=cos（0.48*pi*n1）+cos（0.52*pi*n1）;
subplot（211）;stem（n1,y1）;title（'signalx（n）,0<=n<=9'）;xlabel（'n'）;axis（[0,10,-2.5,2.5]）;
Y1=fft（y1）;
magY1=abs（Y1（1:
1:
6））;
k1=0:
1:
5;
w1=2*pi/10*k1;
subplot（212）;stem（w1/pi,magY1）;title（'sampleofDTFTMagnitude'）;xlabel（'frequency'）;axis（[0,1,0,10]）;
%%part2,取x（n）前10点的高密度频谱（100点）
figure
（2）;
n3=[0:
1:
99];
y3=[y1（1:
1:
10）zeros（1,90）];
subplot（211）;stem（n3,y3）;title（'signalx（n）,0<=n<=9+90zeros'）;xlabel（'n'）;axis（[0,100,-2.5,2.5]）;
Y3=fft（y3）;
magY3=abs（Y3（1:
1:
51））;
k3=0:
1:
50;
w3=2*pi/100*k3;
subplot（212）;stem（w3/pi,magY3）;title（'DTFTMagnitude'）;xlabel（'frequency'）;axis（[0,1,0,10]）;
%%part3
figure（3）;
n=[0:
1:
99];
x=cos（0.48*pi*n）+cos（0.52*pi*n）;
subplot（211）;stem（n,x）;title（'signalx（n）,0<=n<=99'）;xlabel（'n'）;axis（[0,100,-2.5,2.5]）;
X=fft（x）;
magX=abs（X（1:
1:
51））;
k=0:
1:
50;
w1=2*pi/100*k;
subplot（212）;plot（w1/pi,magX）;title（'DTFTMagnitude'）;
10-12已知系统的直接形式：
转换成级联和并联形式，并求3种形式实的单位冲激响应。
clf;
b=[1,-3,11,-27,18];
a=[16,12,2,-4,-1];
[b0,B,A]=dir2cas（b,a）;
[C,B1,A1]=dir2par（b,a）;
N=24;n=1:
N+1;
delta=impseq（0,0,N）;
h1=filter（b,a,delta）;
h2=casfiltr（b0,B,A,delta）;
h3=parfiltr（C,B1,A1,delta）;
figure
（1）,subplot（311）,plot（n,h1）,title（'Directfromstructure'）;
subplot（312）,plot（n,h2）,title（'Cascadestructure'）;
subplot（313）,plot（n,h3）,title（'parallelfromstructure'）;
扩展函数impseq
function[x,n]=impseq（n0,n1,n2）
if（（n0n2）|（n2error（'argumentsmustsatisfyn1<=n0<=n2'）
end
n=[n1:
n2];
x=[（n-n0）==0]
扩展函数dir2cas
文件头:
function[b0,B,A]=dir2cas（b,a）;
%变直接形式为级联形式
%[b0,B,A]=dir2cas（b,a）
%b0=增益系数
%B=包含各因子系数bk的K行3列矩阵
%A=包含各因子系数ak的K行3列矩阵
%a=直接型分子多项式系数
%b=直接型分母多项式系数
%计算增益系数
b0=b
（1）;b=b/b0;
a0=a
（1）;a=a/a0;
b0=b0/a0;
%将分子、分母多项式系数的长度补齐进行计算
M=length（b）;N=length（a）;
ifN>M
b=[bzeros（1,N-M）];
elseifM>N
a=[azeros（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=[b0];
a=[a0];
end
%根据多项式系数利用函数roots求出所有的根
%利用函数cplxpair进行按实部从小到大的成对排序
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
扩展函数dir2par
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
扩展函数casfiltr
functiony=casfiltr（b0,B,A,x）
[K,L]=size（B）;
N=length（x）;
w=zeros（K+1,N）;
w（1,:
）=x;
fori=1:
1:
K
w（i+1,:
）=filter（B（i,:
）,A（i,:
）,w（i,:
））;
end
y=b0*w（K+1,:
）;
扩展函数parfiltr
functiony=parfiltr（C,B,A,x）
[K,L]=size（B）;
N=length（x）;
w=zeros（K+1,N）;
w（1,:
）=filter（C,1,x）;
fori=1:
1:
K
w（i+1,:
）=filter（B（i,:
）,A（i,:
）,x）;
end
y=sum（w）;
扩展函数cplxcomp
functionI=cplxcomp（p1,p2）
I=[];
forj=1:
length（p2）
fori=1:
length（p1）
if（abs（p1（i）-p2（j））<0.0001）
I=[I,i];
end
end
end
I=I';
上机实验三
主要内容：
IIR滤波器设计。
设计Butterworth低通滤波器；设计Chebyshev低通滤波器
10-13设计Butterworth低通滤波器，指标为：
=0.2
，
；
；
wp=0.2*pi;%数字通带截止频率
ws=0.3*pi;%数字阻带截止频率
Rp=1;
As=15;
T=1;
OmegaP=wp*T;
OmegaS=ws*T;
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep*ep））;‘通带衰减值
Attn=1/（10^（As/20））;‘阻带衰减值
%计算
[cs,ds]=afd_butt（OmegaP,OmegaS,Rp,As）;
[b,a]=imp_invr（cs,ds,T）;
[C,B,A]=dir2par（b,a）;
%plot
figure
（1）;
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（221）,plot（w/pi,mag）;title（'magnituderesponse'）;
xlabel（'frequen'）;ylabel（'\H\'）;axis（[0,1,0,1.1]）
set（gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]）;
set（gca,'YTickMode','manual','YTick',[0,Attn,Ripple,1]）;grid;
subplot（223）,plot（w/pi,db）;title（'magnitudeindb'）;
xlabel（'frequen'）;ylabel（'decibels'）;axis（[0,1,-40,5]）
set（gca,'XTickMode','manual','XTick',[0,0.2,0.
升级会员 