《数字信号处理》上机实习报告17.docx
- 文档编号:10765098
- 上传时间:2023-02-22
- 格式:DOCX
- 页数:30
- 大小:661.39KB
《数字信号处理》上机实习报告17.docx
《《数字信号处理》上机实习报告17.docx》由会员分享,可在线阅读,更多相关《《数字信号处理》上机实习报告17.docx(30页珍藏版)》请在冰豆网上搜索。
《数字信号处理》上机实习报告17
《计算机编程与数字信号处理实习》报告
实习目的:
熟练掌握matlab,能够使用该软件完成一定的数字信号处理的任务,并能从中有一定的领悟,达到应用和熟练的程度。
实习任务:
一、从给定的程序(文件包Friday.rar)中,选择一个源程序做详细标注。
(目的:
熟悉Matlab程序)
源程序见matlab文件夹“1”
二、能够利用Matlab熟悉地画图,内容包括:
X、Y坐标轴上的label,每幅图上的title,绘画多条曲线时的legend,对图形进行适当的标注等。
1.在一副图上画出多幅小图;
clearall;%%清除所有变量%%
closeall;%%关闭窗口%%
t=-10:
0.1:
10;%%产生变量序列%%
%%%%%%%%%画出正弦曲线并进行纵横坐标的标注&&&&&&&
subplot(2,2,1);%%图像一%%%
plot(t,sin(t));%%画出正弦曲线%%%
axis([-5,5,-1.5,1.5]);%%限定坐标范围%%
title('正弦曲线');%%添加图像标题%%
xlabel('x');%%%标注X轴%%%
ylabel('y');%%%标注Y轴%%
grid;%%产生网孔%%
%%%画出余弦曲线进行标注%%%
subplot(2,2,2);%%%图像二%%
plot(t,cos(t),'-r','LineWidth',2);%%%画出余弦曲线%%
axis([-5,5,-1.5,1.5]);%%限定坐标长度%%%
title('余弦曲线');%%添加图像标题%%
xlabel('x');%%标注X轴%%
ylabel('y');%%标注Y轴%%
grid;%%产生网孔%%
%%%画出正切曲线进行标注%%%
subplot(2,2,3);%图像三%%%
plot(t,tan(t),'-b','LineWidth',1);%画正切曲线%
axis([-5,5,-1.5,1.5]);%%限定图像范围%%%
title('正切曲线');%%添加图像标题%%
xlabel('x');%%标注X轴%%
ylabel('y');%%标注Y轴%%
%%%画出三角脉冲和矩形脉冲%%%%
subplot(2,2,4);
plot(t,rectpuls(t,2),'-r',t,sawtooth(t,1));%%%三角脉冲和矩形脉冲%%%
axis([-5,5,-1.5,1.5]);%%限定图像范围%%%
title('三角脉冲和矩形脉冲');
xlabel('x');%%标注X轴%%
ylabel('y');%%标注Y轴%%
Legend('rectpuls(t,2)','sawtooth(t,1)');%%%%%对图形进行标注%
grid;%%%添加网孔%%%
print-djpeg-r0Firstwork.jpeg;%%%%保存图片%%%
运行结果:
2.一个三维图像和复数的表示源程序.
igure
(2)
t=linspace(-6,6);
f1=t.^2;
f2=-t.^2;
plot(t,f1,'-k');
holdon;
plot(t,f2,':
r');
text(-2,20,'f1=t^2','FontSize',12,'FontWeight','bold');
text(-2,-20,'f2=-t^2','FontSize',12,'FontWeight','bold');
title('¶þάͼ')
print-djpeg-r02.jpeg
%%%%%%%%%%%%%%
figure(3)
[X,Y]=meshgrid(-3:
0.2:
3);
Z=81.4881+1.2877*X+2.9766*Y;
mesh(X,Y,Z)
holdon
plot3(X,Y,Z,'x','MarkerSize',3);
title('Èýάͼ')
print-djpeg-r03.jpeg
%%%%%%%%%%%%%%
figure(4)
b=3+4i;
compass(b);
title('¸´Êýͼ')
print-djpeg-r04.jpeg
三、计算普通褶积与循环褶积,分别使用时间域与频率域两种方法进行正、反演计算,指出循环褶积计算时所存在的边界效应现象;编写一个做相关分析的源程序。
普通褶积和循环褶积:
clearall;
closeall;
%普通褶积通过矩阵完成
aa=[2234]';
bb=[13553]';
%%%%%%%%%%%%%%%%
%利用conv函数直接计算,对计算结果进行验证
%%%%%%%%%%%%%%%%
aa_bb=conv(aa,bb);
bb_aa=conv(bb,aa);
%%%%%%%%%%%%
N1=length(aa);
N2=length(bb);
N=N1+N2-1;
A=zeros(N,N2);
B=zeros(N,N1);
forj=1:
N2;
A(j:
j+N1-1,j)=aa;
end
forj=1:
N1;
B(j:
j+N2-1,j)=bb;
end
%%%%%%%%%%%%%
%计算结果
AA=A*bb;
BB=B*aa;
figure
(1);
subplot(3,1,1);stem(aa,'.');title('x(n)');%画图
subplot(3,1,2);stem(bb,'.');title('y(n)');%画图
subplot(3,1,3);stem(AA'.');title('褶积');%画图
print-djpeg-r0secondwork1.jpeg%输出
%%%%%%%%%%%%%%%
%普通褶积通过循环褶积的利用MATLAB中的IFFT和FFT实现
aa1=zeros(N,1);
bb1=zeros(N,1);
aa1(1:
N1)=aa;
bb1(1:
N2)=bb;
AA1=fft(aa1);
BB1=fft(bb1);
CC=AA1.*BB1;
CC1=BB1.*AA1;
cc1=ifft(CC1);
cc=ifft(CC);
%循环褶积通过快速傅立叶变换
w=exp(-sqrt(-1)*2*pi/N);
aa2=zeros(N,1);
bb2=zeros(N,1);
aa2(1:
N1)=aa;
bb2(1:
N2)=bb;
F=zeros(N);
fori=1:
N
forj=1:
N
F(i,j)=power(w,(i-1)*(j-1));
end
end
X1=F*aa2;
Y1=F*bb2;
DD=1/N*conj(F)*(X1.*Y1);
DD1=1/N*conj(F)*(Y1.*X1);
%通过程序写出循环矩阵的系数矩阵,求取循环褶积DD=DD1说明褶积具有可交换性
N=5;
%%%%%%%可以改变N的值来验证N=N1+N2时,循环褶积和线性褶积相等;
C1=zeros(N);
C2=zeros(N);
aa3=zeros(N,1);
bb3=zeros(N,1);
aa3(1:
N1)=aa;
bb3(1:
N2)=bb;%%%%%%%%%%%%%%%%%%%双重循环
forj=1:
N
fori=1:
N
k=i+j-1;
ifk>N
k=k-N;
end
C1(k,j)=aa3(i);
C2(k,j)=bb3(i);
end
end
EE=C1*bb3;
EE1=C2*aa3;
figure
(2);
subplot(3,1,1);stem(aa,'.');title('x(n)');%画图
subplot(3,1,2);stem(bb,'.');title('y(n)');%画图
subplot(3,1,3);stem(EE,'.');title('褶积');%画图
print-djpeg-r0secondwork2.jpeg
运行结果:
边界效应:
两个离散的序列离散x(n)和y(n),他们的长度分别为N1和N2,如果循环褶积的长度N>N1+N2-1,则循环褶积和线性褶积的值相等。
否则,循环褶积和线性褶积的值不相等。
而在循环褶积的系数矩阵中,把第一列离散的序列依次折叠刀第二列,第三列等,这个时候,当循环褶积的N=N1+N2-1时,循环褶积的系数矩阵和线性褶积的系数矩阵不相同,循环褶积系数矩阵边沿是不是零,而线性褶积系数矩阵边沿是零,但是他们运算结果一致,这就是循环褶积在计算过程中的边界效应。
线性相关和循环相关源程序:
clearall;%%%%%%%%%%%%%%%%%%相关分析%%%%%%%%%%%
closeall;
aa=[1,2,3,4]';
bb=[1,2,5,3]';
%%%%%%%%%%%%%%线形相关分析
N1=length(aa);
N2=length(bb);
N=N1+N2-1;
N=4;
A=zeros(N);
aa3=zeros(N,1);
bb3=zeros(N,1);
aa3(1:
N1)=aa;
bb3(1:
N2)=bb;
forj=1:
N;
A(j:
N,N-j+1)=aa3(1:
N-j+1);
end
CC1=A*bb3;
figure
(1)
subplot(3,1,1);stem(aa,'.');title('x(n)');%画图
subplot(3,1,2);stem(bb,'.');title('y(n)');%画图
subplot(3,1,3);stem(CC1,'.');title('线性相关');%画图
print-djpeg-r0secondwork3.jpeg;
%%%%%%%%%%%%%%%循环相关分析
N1=length(aa);
N2=length(bb);
N=N1+N2-1;
N=4;
A1=zeros(N);
aa4=zeros(N,1);
bb4=zeros(N,1);
aa4(1:
N1)=aa;
bb4(1:
N2)=bb;
forj=1:
N;
A1(1:
N-j+1,j)=aa4(j:
N);
A1(N-j+2:
N,j)=aa(1:
j-1);
end
CC2=A1*bb3;
figure
(2)
subplot(3,1,1);stem(aa,'.');title('x(n)');%画图
subplot(3,1,2);stem(bb,'.');title('y(n)');%画图
subplot(3,1,3);stem(CC2,'.');title('循环相关');%画图
print-djpeg-r0secondwork4.jpeg;
四、设计一个病态(矩阵)系统,分析其病态程度;找出对应的解决方法(提示:
添加白噪因子)。
病态矩阵源程序:
clearall;
closeall;
a=[1.29690.8648;0.21610.1441];%转换前的病态矩阵,对它一个小扰动会造成线性方程组解发生巨大变化
N=length(a);
w=exp(-sqrt(-1)*2*pi/N);
fori=1:
N
forj=1:
N
F(i,j)=power(w,(i-1)*(j-1));end
end
k=F*a*inv(F)+0.001;%对A加白噪因子
m=inv(F)*k*F;%还原矩阵
A=real(m);%得到加白噪因子后的矩阵,利用它解AA*X=B
B=[0.86480.1441]';
X1=inv(A)*B%AA*X=B的解
AA=A+0.001;BB=B+0.001;%对AA,b进行小小的扰动
X2=inv(AA)*BB%扰动后的解
ans=X2-X1%这是一个接近于0的数组,说明线性方程组解变化不是很大,AA非病态的
运行结果:
这个结果的相差很小,说明病态矩阵经过添加白噪因子之后,已经不在是病态,给其一个很小的扰动,其线性方程的结果变化很小。
病态矩阵是求解方程组时对数据的小扰动很敏感的矩阵。
解线性方程组Ax=b时,若对于系数矩阵A及右端项b的小扰动,方程组的解x与原方程组Ax=b的解差别很大,则称矩阵A为病态矩阵。
方程组的近似解x一般都不可能恰好使剩余r=b-Ax为零,这时x亦可看作小扰动问题Ax=b-r(即)的解,所以当A为病态时,即使剩余r很小,仍可能得到一个与真解相差很大的近似解。
解决这个问题的方法就是给病态矩阵添加一个白噪因子,这样就可以解决其本身的病态问题。
五、设计一个一维滤波处理程序(1、分别做低通、高通、带通、带阻等理想滤波器进行处理;2、窗函数)。
%一维理想滤波器设计
closeall
clearall
f1=30;f2=45;%设定滤波器的两个频率,f1为低通滤波器的高通频率、带通滤波器的低通频率,f3为高通滤波器的高截频率、带通滤波器的高通频率
d=0.01;%设定采样间隔
t=0:
d:
1;%建立时间向量
x=sin(2*pi*10*t)+2*sin(2*pi*25*t)+8*cos(2*pi*45*t);%由时间向量构成的含有高、低、中频的原始信号
lx=length(x);%时间向量长度
fc=1/(2*d);%截止频率
f=0:
2*fc/lx:
2*fc;%求出频谱对应的频率序列
FX=fft(x);%求原始信号的频谱
%理想低通滤波器设计
H1=zeros(1,lx);
fori=1:
lx
iff(i)<=f1|f(i)>=(2*fc-f1)
H1(i)=1;
elseH1(i)=0;
end
end
DFX1=FX.*H1;
DX1=ifft(DFX1);
figure
(1);
subplot(2,2,1),plot(H1);axis([0100-0.51.5]);title('低通H1(f)');
subplot(2,2,2),plot(t,x);title('原始信号');
subplot(2,2,3),stem(abs(FX),'.');;title('原始信号频谱');
subplot(2,2,4),stem(abs(DFX1),'.');;title('经低通滤波后的低频信号的频谱');
print-djpeg-r0低通滤波器.jpeg
%理想高通滤波器设计
H2=zeros(1,lx);
fori=1:
lx
iff(i)<=2*f1&&f(i)>=f1
H2(i)=1;
elseH2(i)=0;
end
end
DFX2=FX.*H2;
DX2=ifft(DFX2);
figure
(2);
subplot(2,2,1),plot(H2);axis([0100-0.51.5]);title('高通H2(f)');
subplot(2,2,2),plot(t,x);;title('含有高低中频的原始信号');
subplot(2,2,3),stem(abs(FX),'.');title('原始信号频谱');
subplot(2,2,4),stem(abs(DFX2),'.');title('经高通滤波后的高频信号的频谱');
print-djpeg-r0高通滤波器.jpeg
%理想带通滤波器设计
H3=zeros(1,lx);
fori=1:
lx
iff(i)<=f2&&f(i)>=f1||(f(i)<=(fc+f2-f1)&&(f(i))>=(fc));
H3(i)=1;
elseH3(i)=0;
end
end
DFX3=FX.*H3;
DX3=ifft(DFX3);
figure(3);
subplot(2,2,1),plot(H3);axis([0100-0.51.5]);title('带通H3(f)');
subplot(2,2,2),plot(t,x);title('含有高低中频的原始信号');
subplot(2,2,3),stem(abs(FX),'.');title('原始信号频谱');
subplot(2,2,4),stem(abs(DFX3),'.');title('经带通滤波后的中频信号的频谱');
print-djpeg-r0带通滤波器.jpeg
%理想带阻滤波器设计
f1=20;f2=40;
H4=zeros(1,lx);
fori=1:
lx
if(f(i))<=f2&&(f(i))>=f1||((f(i))<=(f2+fc)&&(f(i))>=(f1+fc))
H4(i)=0;
elseH4(i)=1;
end
end
DFX4=FX.*H4;
DX4=ifft(DFX4);
figure(4);
subplot(2,2,1),plot(H4);axis([0100-0.51.5]);title('带阻H4(f)');
subplot(2,2,2),plot(t,x);title('含有高低中频的原始信号');
subplot(2,2,3),stem(abs(FX),'.');title('原始信号频谱');
subplot(2,2,4),stem(abs(DFX4),'.');title('经带阻滤波后的含有高低频信号的频谱');
print-djpeg-r0带阻滤波器.jpeg
窗函数:
运行结果:
六、设计一个二维滤波处理程序(分别做低通、高通等处理)。
%%%%%%二维滤波%%%%%%%%%
clearall;
closeall;
N=50;
fx1=15;fx2=35;fy1=15;fy2=35;
x=linspace(0,0.25,N);
y=linspace(0,0.25,N);
[X,Y]=meshgrid(x,y);
b=rand(N);
T=exp(2*sin(2*pi*X)+3*cos(2*pi*Y))+4*b;
fori=1:
N
ft1(i,1:
N)=fft2(T(i,1:
N));
end
forj=1:
N
ft(1:
N,j)=fft2(ft1(i,1:
N));
end
AA=size(ft);
fxc=2*N/0.25;
fyc=2*N/0.25;
%%%%%%%%%%%%%%%%%%%%%%%%%%低通滤波器
f1=zeros(AA
(1):
AA
(2));
fori=1:
fx1
forj=1:
fy1
f1(i,j)=1;
end
end
fori=1:
fx1
forj=fy1:
N
f1(i,j)=1;
end
end
fori=fx2:
N
forj=1:
15
f1(i,j)=1;
end
end
fori=fx2:
N
forj=fy2:
N
f1(i,j)=1;
end
end
tt1=ft*f1;
f2=zeros(AA
(1):
AA
(2));
fori=1:
fx1
forj=1:
fy1
f2(i,j)=1;
end
end
fori=1:
fx1
forj=fy2:
N
f2(i,j)=1;
end
end
fori=fx2:
N
forj=1:
fy1
f2(i,j)=1;
end
end
fori=fx2:
N
forj=fy2:
N
f2(i,j)=1;
end
end
tt=f2*tt1;
%%%%%%%%%%%%%%%高通%%%%
%%%%%%%%%%%%%%%%%%%%
f3=zeros(AA
(1):
AA
(2));
fori=1:
50
forj=1:
50
f3=1;
end
end
f4=f3-f2;
t3=f4*ft;
t4=t3*f4;
figure
(1);
surf(abs(T));title('原始信号');
figure
(2)
surf(abs(ft));title('原始信号频谱');
figure(3)
surf(abs(tt));title('低通滤波');
figure(4)
surf(abs(t4));title('高通滤波');
print-djpeg-r0二维滤波器.jpeg
运行结果:
原始信号:
原始信号频谱:
二维低通滤波:
二维高通滤波:
七、验证时间域的循环褶积对应的是频率域的乘积;线性褶积则不然。
closeall;
clearall;
aa=[2234]';
bb=[1355]';
N1=length(aa);
N2=length(bb);
N=N1+N2-1;
%%%%%%%%%%%系数矩阵计算循环褶积%%%%%
aa1=zeros(N,1);
bb1=zeros(N,1);
aa1(1:
N1)=aa;
bb1(1:
N2)=bb;
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)=aa1(i);
C2(k,j)=bb1(i);
end
end
EE=C1*bb1
%%%%%%%%利用FFt计算褶积%%%%
aa3=fft(aa1);
bb3=fft(bb1);
M=aa3.*bb3;
EE1=ifft(M)
max=abs(EE1-EE)
%%%%%max=0说明循环褶积时间域和频率域对应%%%
figure
(1);
subplot(4,1,1);stem(aa,'.');title('x(n)');%画图
subplot(4,1,2);stem(bb,'.');title('y(n)');%画图
subplot(4,1,3);stem(EE,'.');title('时间域计算循环褶积');%画图
subplot(4,1,4);stem(EE1,'.');title('频率域计算循环褶积');%画图
print-djpeg-r0sixwork1.jpeg
%%%%%%计算线形褶积%%%%
A=zeros(N,N2);
B=zeros(N,N1);
forj=1:
N2;
A(j:
j+N1-1,j)=aa;
end
forj=1:
N1;
B(j:
j+N2-1,j)=bb;
end
%%%%%%%%%%%%%%%%%%
%计算结果
AA=A*bb;
AA1=fft(aa);
BB1=fft(bb);
CC=AA1.*BB1
AA2=ifft(CC)
figure
(2);
subplot(4,1,1);stem(aa,'.');title('x(n)');%画图
subplot(4,1,2);stem(bb,'.');title('y(n)');%画图
subplot(4,1,3);stem(AA,'.');title('时间域计算普通褶积');%画图
subplot(4,1,4);stem(AA2,'.');title('频率域计算普通褶积');%画图
print-djpeg-r0
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号处理 数字信号 处理 上机 实习 报告 17