计算机编程与数字信号处理实习.docx
- 文档编号:29309255
- 上传时间:2023-07-22
- 格式:DOCX
- 页数:62
- 大小:1.20MB
计算机编程与数字信号处理实习.docx
《计算机编程与数字信号处理实习.docx》由会员分享,可在线阅读,更多相关《计算机编程与数字信号处理实习.docx(62页珍藏版)》请在冰豆网上搜索。
计算机编程与数字信号处理实习
计算机编程与数字信号处理实习
能够利用Matlab熟悉地画图,内容包括:
X、Y坐标轴上的label,每幅图上的title,绘画多条曲线时的legend,对图形进行适当的标注等。
(1)在一副图上画出多幅小图;
clc;clear;close;
x=linspace(0,5,200);
y1=sin(3*x);
y2=sin(x.^3);
y3=exp(x-2)
y4=sin(2*x)+cos(3*x);
fork=1:
4;
ifk==1
subplot(2,2,k)
plot(x,y1)
holdon;
title('x与y1')
elseifk==2
subplot(2,2,k)
plot(x,y2)
holdon;
title('x与y2')
elseifk==3
subplot(2,2,k)
plot(x,y3)
holdon;
title('x与y3')
elseifk==4
subplot(2,2,k)
plot(x,y4)
holdon;
title('x与y4')
end
end
print-djpeg-r0xinhao1.jpeg
(2)画出一组二维图形;(以分段函数为例)
clc;clear;close;
x=-3:
0.01:
3;
y1=zeros(size(x));
y2=zeros(size(x));
y3=zeros(size(x));
N=length(x);
fork=1:
N
ifx(k)<-1&x(k)>=-3;
y1(k)=(-x(k).^2-4*x(k)-3)/2;
elseifx(k)>=-1&x(k)<1;
y2(k)=-x(k).^2+1;
elsex(k)<=3&x(k)>=1;
y3(k)=(-x(k).^2+4*x(k)-3)/2;
end
end
y=y1+y2+y3;
plot(x,y)
xlabel('x');ylabel('y');
print-djpeg-r0xinhao2.jpeg
(3)画出一组三维图形;
clear
x=0:
0.2:
3*pi;
y=0:
0.25:
5*pi;
[xx,yy]=meshgrid(x,y);
z1=sin(xx).*sin(yy);
x=-3:
0.25:
3;
y=x;
[xx,yy]=meshgrid(x,y);
z2=xx-0.5*xx.^3+0.2*yy.^2+1;
x=-8:
0.5:
8;
y=x;
[xx,yy]=meshgrid(x,y);
r=sqrt(xx.^2+yy.^2)+eps;
z3=sin(r)./r;
subplot(2,2,1),mesh(z1);
title('z1')
subplot(2,2,2),mesh(z2);
title('z2')
subplot(2,2,3),waterfall(z2);
title('z2')
subplot(2,2,4),mesh(z3);
title('z3')
print-djpeg-r0xinhao3.jpeg
说明:
第一个图和第二个图中对x1赋值方式不同,结果相同。
可以看出,无论是第一段中对X1的定义,还是第二段中对X1的定义。
只要X1的维度为10,就符合surf的输入条件。
且第一段中X1的范围就是常数9,X2中同样不变,所以画出的图像当然是一样的。
(4)画出复数的实部与虚部。
clc;clear;close;
x=20;
y=rand(x,1)+rand(x,1)*i;
m=real(y);
n=imag(y);
subplot(2,2,1)
plot(y,'.');title('复数');
holdon;gridon;
subplot(2,2,2)
plot(m,n,'.');title('复数(实部为横,虚部为纵)');
xlabel('实部');ylabel('虚部');
holdon;gridon;
subplot(2,2,3)
plot(m,'.');title('复数实部');
holdon;gridon;
subplot(2,2,4)
plot(n,'.');title('复数虚部');
holdon;gridon;
print-djpeg-r0xinhao4.jpeg
(5)完成对一个源程序进行详细注释。
(仍以分段函数为例)
clc;clear
x=-3:
0.01:
3;%从-3到+3以0.01为步长给x赋值
y1=zeros(size(x));%初始化y1的值
y2=zeros(size(x));%初始化y1的值
y3=zeros(size(x));%初始化y1的值
N=length(x);%以x的维度给N赋值
fork=1:
N%k从1循环到N
ifx(k)<-1&x(k)>=-3;%分段赋值
y1(k)=(-x(k).^2-4*x(k)-3)/2;
elseifx(k)>=-1&x(k)<1;
y2(k)=-x(k).^2+1;
elsex(k)<=3&x(k)>=1;
y3(k)=(-x(k).^2+4*x(k)-3)/2;
end
end
y=y1+y2+y3;%各段叠加成分段函数y
plot(x,y)%以x为横坐标,y为纵坐标作图
二、计算普通褶积与循环褶积,分别使用时间域与频率域两种方法进行正、反演计算,指出循环褶积计算时所存在的边界效应现象;编写一个做相关分析的源程序。
(1)作线性褶积,时间域正演,频率域反演(无)
clc;clear;close;
x=[012]';
y=[0123]';
xy0=conv(x,y);%利用公式计算线性褶积
yx0=conv(y,x);%利用公式计算线性褶积
%%%%%%编程计算%%%%%%
N1=length(x);%
N2=length(y);%
NN=N1+N2-1;%
%xy=zeros(1,NN);%
%yx=zeros(1,NN);%
A1=zeros(NN,N2);%
A2=zeros(NN,N1);%
forj=1:
N2%
fori=j:
j+N1-1%
A1(i,j)=x(i-j+1);%
end%
end%
forj=1:
N1%
fori=j:
j+N2-1%
A2(i,j)=y(i-j+1);%
end%
end%
xy1=A1*y;%
yx1=A2*x;%
%%%%%%结果验证%%%%%%
xy0
yx0
xy1
yx1
max(abs(xy0-yx0))%说明褶积的可交换性
max(abs(xy1-yx1))%说明褶积的可交换性
max(abs(xy0-xy1))%说明公式计算和编程计算的结果一致性
%%%作图%%%
subplot(2,2,1)
stem(x)
gridon;holdon;
subplot(2,2,2)
stem(y)
gridon;holdon;
subplot(2,2,3)
stem(xy0)
gridon;holdon;
subplot(2,2,4)
stem(yx0)
gridon;holdon;
print-djpeg-r0xinhao.jpeg
运行结果:
xy0=
0
0
1
4
7
6
yx0=
0
0
1
4
7
6
xy1=
0
0
1
4
7
6
yx1=
0
0
1
4
7
6
ans=
0
ans=
0
ans=
0
(2)作循环褶积,时间域正演与频率域反演
clc;clear;close;
x=[012]';
y=[0123]';
xy0=conv(x,y);%利用公式计算线性褶积
yx0=conv(y,x);%利用公式计算线性褶积
%%%%%%自带公式频率域反演计算%%%%%%
N1=length(x);%
N2=length(y);%
N=N1+N2-1;%
NN=N1+N2-1;%
w=exp(-sqrt(-1)*2*pi/N);%
x0=zeros(N,1);%需要将原始数组转化为相同长度的两个数组(补零)。
y0=zeros(N,1);%需要将原始数组转化为相同长度的两个数组(补零)。
x0(1:
N1)=x;%
y0(1:
N2)=y;%
X0=fft(x0);%
Y0=fft(y0);%
xy2=ifft(X0.*Y0);%
yx2=ifft(Y0.*X0);%
%%%%%%频率域反演%%%%%%
F=zeros(N);%
fori=1:
N%
forj=1:
N%
F(i,j)=power(w,(i-1)*(j-1));%
end%
end%
X1=F*x0;%
Y1=F*y0;%
xy3=1/N*conj(F)*(X1.*Y1);%inv(F)=conj(F)/N
yx3=1/N*conj(F)*(Y1.*X1);%
%%%%%%时间域直接计算的方法%%%%%%
C1=zeros(N);%
C2=zeros(N);%
forj=1:
N%
fori=1:
N%
k=i+j-1;%
ifk>N%
k=k-N;%
end%
C1(k,j)=x0(i);%
C2(k,j)=y0(i);%
end%
end%
xy4=C1*y0;%
yx4=C2*x0;%
%%%%%%结果验证%%%%%%
xy2
yx2
xy3
yx3
xy4
yx4
max(abs(xy2-yx2))
max(abs(xy3-yx3))
max(abs(xy4-yx4))
max(abs(xy2-xy3))
max(abs(xy2-xy4))
max(abs(xy2(1:
NN)-xy0))%验证结果:
可以利用循环褶积的方法计算线性褶积。
%%%作图%%%
subplot(2,2,1)
stem(x)
gridon;holdon;
subplot(2,2,2)
stem(y)
gridon;holdon;
subplot(2,2,3)
stem(xy2)
gridon;holdon;
subplot(2,2,4)
stem(yx2)
gridon;holdon;
print-djpeg-r0xinhao.jpeg
运行结果:
xy2=
0
0
1
4
7
6
yx2=
0
0
1
4
7
6
xy3=
0.0000-0.0000i
-0.0000-0.0000i
1.0000-0.0000i
4.0000+0.0000i
7.0000+0.0000i
6.0000-0.0000i
yx3=
0.0000-0.0000i
-0.0000-0.0000i
1.0000-0.0000i
4.0000+0.0000i
7.0000+0.0000i
6.0000-0.0000i
xy4=
0
0
1
4
7
6
yx4=
0
0
1
4
7
6
ans=
0
ans=
0
ans=
0
ans=
7.8221e-015
ans=
0
ans=
0
这些ans=0说明了各种方法计算得到的褶积结果是相等的,并且说明了褶积的可交换性。
(3)循环褶积计算时所存在的边界效应现象
变量赋值如上所示,程序不再累述
分别以xy1,xy4代表线性褶积和循环褶积
当N=4时
结果如下:
xy1=
0
0
1
4
7
6
xy4=
7
6
1
4
当N=5时
结果如下:
xy1=
0
0
1
4
7
6
xy4=
6
0
1
4
7
当N=6时
结果如下:
xy1=
0
0
1
4
7
6
xy4=
0
0
1
4
7
6
当N=7时
结果如下:
xy1=
0
0
1
4
7
6
xy4=
0
0
1
4
7
6
0
由此可见,当N≥N1+N2-1时,可以利用循环褶积计算线性褶积,N1+N2-1可以看做一个边界值,这种现象称之为循环褶积计算时所存在的边界效应现象。
(4)相关分析的示例
clc;clear;close;
x=[012]';
y=[0123]';
N1=length(x);%
N2=length(y);%
NN=N1+N2-1;%
y0=zeros(NN,1);%构造相关乘积用的列向量
y0(N1:
NN)=y;%把y的值赋给y0
A1=zeros(NN,N2);%线性褶积矩阵
A2=zeros(NN,N1);%线性相关矩阵
forj=1:
N2%
fori=j:
j+N1-1%
A1(i,j)=x(i-j+1);%
end%
end%
forj=1:
NN
fori=NN-j+1:
min(NN-j+1+2,NN)%
A2(i,j)=x(i-(NN-j));%
end
end%
xyZhe=A1*y;%线性褶积
xyXiang=A2*y0;%线性相关
xyXiangN=xcorr(x,y)%由matlab自带的函数计算相关
%%%作图%%%
subplot(2,2,1)
stem(x,'r','.');
xlabel('n');ylabel('x');
title('x向量');
gridon;
subplot(2,2,2)
stem(y,'r','.');
xlabel('n');ylabel('y');
title('y向量');
gridon;
subplot(2,2,3)
stem(xyXiang,'r','.');
xlabel('n');ylabel('xyXiang');
title('自己计算的线性相关');
gridon;
subplot(2,2,4)
stem(xyXiangN,'r','.');
xlabel('n');ylabel('xyXiangN');
title('xcorr计算得到的相关');
gridon;
print-djpeg-r0xiangguang.jpeg
运行结果:
由此可见,编程用循环得到的相关结果和Matlab提供的xcorr函数计算的结果一致,只不过,xcorr计算得到的结果是奇数项表示的。
三、设计一个病态(矩阵)系统,分析其病态程度;找出对应的解决方法(提示:
添加白噪因子)。
PartOne
病态矩阵和及其病态程度的判断:
组别的变量有冗余的矩阵,或者说特征值有0的矩阵就是病态矩阵。
但是在实际工作中,判断一个矩阵是否为病态矩阵是相对的。
通常通过用小扰动后解的改变程度来判断。
因此定义病态矩阵是求解方程组时对数据的小扰动很敏感的矩阵。
解线性方程组Ax=b时,若对于系数矩阵A及右端项b的小扰动
,方程组
的解x与原方程组Ax=b的解差别很大,则称矩阵A为病态矩阵。
方程组的近似解x一般都不可能恰好使剩余r=b-Ax为零,这时x亦可看作小扰动问题Ax=b-r(即
)的解,所以当A为病态时,即使剩余r很小,仍可能得到一个与真解相差很大的近似解。
例如,取
如果取
作为近似解,剩余
已经很小,但x与真解
相差还是很大。
用
来衡量矩阵的病态程度,K(A)称为A的条件数,它很大时,称A为病态,否则称良态;K(A)愈大,A的病态程度就愈严重。
计算得
,这说明A的病态程度很大。
我们对A,b都做一个小小的扰动,
那么
的解
与原来的真解
差别很大,由此说明A是病态矩阵。
如果用单位矩阵来做比较的话:
它的解x就等于b
我们再对E,b做一个小小的扰动,
那么
的解
与真解
完全相等,这说明单位矩阵不是病态的,再用cond函数计算得
,结论就是E不是病态矩阵。
因此解决病态问题的关键是在其扰动的条件下,仍然保持解的不变,这就需要添加白噪因子(可以看做第一次扰动)。
PartTwo
解决病态矩阵方程组AA*X1=C0问题
clc;clear;close;
A=[1221]';
D=[1357]';
B=[1234]';
N=length(A);
F=zeros(N);%[w^(i*j)]
AA=zeros(N);
DD=zeros(N);
w=exp(-sqrt(-1)*2*pi/N);
fori=1:
N
forj=1:
N
AA(i,j)=A(mod(i-j,N)+1);%
DD(i,j)=D(mod(i-j,N)+1);
F(i,j)=power(w,(i-1)*(j-1));%DFT.%w^((i-1)*(j-1))
end
end
FA=F*A;
FB=F*B;
FAA=F*AA;%定义FAA
FDD=F*DD;%定义FDD
C0=AA*B;%=ifft((fft(A)).*fft(B));
%C1=real(conj(F)/N*(FA.*FB));%inv(F)=conj(F)/N;
%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算病态矩阵方程式AA*X1=C0
%以及良态矩阵方程式DD*X2=C0
%%%%%%%%%%%%%%%%%%%%%%%%%%
%AA*X1=C0
%F*AA*inv(F)*F*X1=F*C0
%F*AA*((1/N)*conj(F))*F*X1=F*C0
%(1/N)*conj(F)*AA*F*AA*conj(F)*F*X1=conj(F)*AA*F*C0
%(1/N)*C_FA*F*AA*conj(F)*F*X1=C_FA*FC
%(1/N)*C_FA*FAA*conj(F)*F*X1=FC2
%(1/N)*(C_FA*FAA+0.001*eye(N))*conj(F)*F*X1=FC2%添加白噪因子
%(1/N)*FA2*conj(F)*F*X1=FC2
%FA2*X1=FC2
%X=inv(FA2)*FC2
FC=F*C0;
C_FA=conj(FAA);
FC2=C_FA*FC;
FA2=C_FA*FAA+0.001*eye(N);
X1=inv(FA2)*FC2
REALX1=inv(AA)*C0
%BX=inv(AA+0.001*eye(N))*C0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DD*X2=C0
%F*DD*inv(F)*F*X2=F*C0
%F*DD*((1/N)*conj(F))*F*X2=F*C0
%(1/N)*conj(F)*DD*F*DD*conj(F)*F*X2=conj(F)*DD*F*C0
%(1/N)*C_FD*F*DD*conj(F)*F*X2=C_FD*FC
%(1/N)*C_FD*FDD*conj(F)*F*X2=FC2
%(1/N)*(C_FD*FDD+0.001*eye(N))*conj(F)*F*X2=FC2%添加白噪因子
%(1/N)*FD2*conj(F)*F*X2=FC2
%FD2*X2=FC2
%X=inv(FD2)*FC2
FC=F*C0;
C_FD=conj(FDD);
FC2=C_FD*FC;
FD2=C_FD*FDD+0.001*eye(N);
X2=inv(FD2)*FC2
REALX2=inv(DD)*C0
运行结果:
X1=
-1.9995+0.0006i
5.0004-0.0008i
0.0000+0.0000i
6.9996+0.0001i
REALX1=
Inf
Inf
Inf
Inf
X2=
0.6875+0.0000i
1.1875+0.0000i
1.1875-0.0000i
0.6875-0.0000i
REALX2=
0.6875
1.1875
1.1875
0.6875
结果分析:
病态矩阵AA添加白噪因子前解为Inf,添加白噪因子后方可求解。
良态矩阵DD添加白噪因子前解与添加白噪因子后的解是相等的。
四、设计一个一维滤波处理程序(1、分别做低通、高通、带通、带阻等理想滤波器进行处理;2、窗函数法或者镶边法设计滤波器,以消除Gibbs现象)。
直接设计的理想滤波器(使用时会出现Gibbs现象):
clc;clear;close;
f1=15,f2=35;%设定滤波器的两个频率
d=0.01;%采样间隔
t=0:
d:
1;%建立时间向量
N=length(t);%时间向量长度
x=5*rand(size(t))-2;%随机信号,含中高低频
fc=1/(2*d);%截止频率
f=0:
(1/d)/N:
(1/d);%求出频谱对应的频率序列
FX=fft(x);%求原始信号的频谱
%理想低通滤波器设计
H1=zeros(1,N);
fori=1:
N
iff(i)<=f1|f(i)>=(1/d)-f1
H1(i)=1;
elseH1(i)=0;
end
end
FX1=FX.*H1;%频率滤波
X1=ifft(FX1);%频率滤波后的反演信号
%理想高通滤波器设计
H2=zeros(1,N);
fori=1:
N
if(f(i)>=f2)&(f(i)<(1/d)-f2)
H2(i)=1;
elseH2(i)=0;
end
end
FX2=FX.*H2;
X2=ifft(FX2);
%理想带通滤波器设计
H3=zeros(1,N);
fori=1:
N
if(f(i)<=f2&f(i)>=f1)|(f(i)>2*fc-f2&f(i)<(1/d)-f1);
H3(i)=1;
elseH3(i)=0;
end
end
FX3=FX.*H3;
X3=ifft(FX3);
%理想带阻滤波器设计
H4=zeros(1,N);
fori=1:
N
if(f(i
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 计算机 编程 数字信号 处理 实习