matlab上机实验 数字信号处理实验课必备.docx
- 文档编号:10272124
- 上传时间:2023-05-24
- 格式:DOCX
- 页数:31
- 大小:71.73KB
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((n0 error('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.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab上机实验 数字信号处理实验课必备 matlab 上机 实验 数字信号 处理 必备