《数字信号处理》上机实习报告 5.docx
- 文档编号:8115541
- 上传时间:2023-01-28
- 格式:DOCX
- 页数:22
- 大小:649.93KB
《数字信号处理》上机实习报告 5.docx
《《数字信号处理》上机实习报告 5.docx》由会员分享,可在线阅读,更多相关《《数字信号处理》上机实习报告 5.docx(22页珍藏版)》请在冰豆网上搜索。
《数字信号处理》上机实习报告5
计算机编程与数字信号处理实习报告
一、实习目的
熟悉Matlab编程,并熟悉常用的信号处理手段,加深信号分析的课程知识。
二、实习内容
(一)从老师所给的源程序中文件中,选择合适的一个源程序,仔细学习程序的相关语句和做法,最后做详细标注。
以达到认识了解,从而熟练应用Matlab的目的,对今后的相关知识和学习做铺垫。
具体标注有M文件。
见文件Gibbs.m
(二)、能够利用Matlab熟悉地画图,内容包括:
X、Y坐标轴上的label,每幅图上的title,绘画多条曲线时的legend,对图形进行适当的标注等。
(1)在一副图上画出多幅小图;
(2)画出一组二维图形;
(3)画出一组三维图形;(4)画出复数的实部与虚部。
(5)完成对一个源程序进行详细注释。
(三)、计算普通褶积与循环褶积,分别使用时间域与频率域两种方法进行正、反演计算,指出循环褶积计算时所存在的边界效应现象;编写一个做相关分析的源程序。
(四)、设计一个病态(矩阵)系统,分析其病态程度;找出对应的解决方法(提示:
添加白噪因子)。
(五)、设计一个一维滤波处理程序(1、分别做低通、高通、带通、带阻等理想滤波器进行处理;2、窗函数)。
(六)、设计一个二维滤波处理程序(分别做低通、高通等处理)。
(七)、验证时间域的循环褶积对应的是频率域的乘积;线性褶积则不然。
(八)、请用通俗、易懂的语言说明数字信号处理中的一种性质、一条定理或一个算例(顺便利用Matlab对其进行实现)。
三、实习要求
(1)对每个问题进行编程计算,给出详细的注释;
(2)分析相关原理及编程思路;
(3)绘图显示每个问题的计算结果;
(4)编写总的实习报告(在本次实习最后一天的中午之前必须提交该报告)。
四、实习日记
1.6月21日,对matlab语言认识程度不深,只会一些最基本的用法,一天的时间都在研究matlab的编程,初步掌握了matlab编程过程。
2.6月22日,目标完成实习内容的第二题。
由于都是matlab的基础题目,完成起来相对简单,上午的时间就完成了所有小题。
下午再次去机房,紧接着进行第三题的研究。
3.6月23日,目标完成第三题。
通过昨天下午的初步研究,大致了解了编程的原理和方法,所以直接进行编程并验证操作结果。
经过一上午的努力,完成了第三题。
4.6月24日,开始进行第四题的研究。
通过之前数字信号处理课程上对病态矩阵的认识,结合matlab语言进行编程,上午基本完成病态矩阵的设计程序。
下午在宿舍继续进行操作,经过不断的验证和完善,第四道题也完成了。
晚上对第五题进行了初步的分析研究。
5.6月25日,目标完成第五道题。
上午在机房进行编程操作,在掌握了一位滤波器的编程原理之后,四种滤波器(低通、高通、带通、带阻)的编程就变得非常简单了。
下午回到宿舍将程序进行完善,完成了第四题。
6.6月26日和6月27日是星期六和星期天,没有去机房实习,但是自己在宿舍有对后面几道题的初步分析研究。
7.6月28日,目标完成第六题。
二维滤波器较一维滤波器,难度增加了不少,程序也比一维滤波器复杂很多,上午基本上了解了二维滤波器的原理及编程步骤,有了初步的编程思路。
下午回到宿舍继续进行此题的编程,由于已经掌握了编程的思路,在经过多次验证并对程序进行多次修改后,此题也被攻克。
8.6月29日,完成第六题的编程。
第七题相对简单,完成的也比较轻松。
9.7月2日,最后一天实习,完成了第八题,并对每天的实验记录进行汇总,编写实验报告。
五、具体编程内容
(二)、能够利用Matlab熟悉地画图,内容包括:
X、Y坐标轴上的label,每幅图上的title,绘画多条曲线时的legend,对图形进行适当的标注等。
(1)在一副图上画出多幅小图
源程序:
wd2_1.m
(2)画出一组二维图形
源程序:
wd2_2.m
(3)画出一组三维图形
源程序:
wd2_3.m
(4)画出复数的实部与虚部
源程序:
wd2_4.m
(5)完成对一个源程序进行详细注释
程序:
x=[0:
0.01:
5]*pi%建立数组x
y=sin(x)+i*cos(x)%建立复数方程y
subplot(2,1,1)%画两行一列组图中的第一幅图
plot(real(y))%画出复数的实部
title('复数y=sin(x)+i*cos(x)的实部')%标题
subplot(2,1,2)%画两行一列组图中的第二幅图
plot(imag(y))%画出复数的虚部
title('复数y=sin(x)+i*cos(x)的虚部')%标题
(三)、计算普通褶积与循环褶积,分别使用时间域与频率域两种方法进行正、反演计算,指出循环褶积计算时所存在的边界效应现象;编写一个做相关分析的源程序。
此题关键是对普通褶积(线性褶积)和循环褶积运算方式的掌握和对“for”循环语句的运用。
线性褶积1:
x=[123]';
y=[234]';
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;%证明褶积具有可交换性
运行结果:
xy1=yx1=2
7
16
17
12
线性褶积2:
x=[123]';
y=[234]';
N1=length(x);
N2=length(y);
N=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);
运行结果:
xy2=yx2=2.0000
7.0000
16.0000
17.0000
12.0000
循环褶积1:
x=[123]';
y=[234]';
N1=length(x);
N2=length(y);
N=N1+N2-1;
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);
运行结果:
xy3=yx3=2.0000-0.0000i
7.0000-0.0000i
16.0000-0.0000i
17.0000+0.0000i
12.0000+0.0000i
循环褶积2:
x=[123]';
y=[234]';
N1=length(x);
N2=length(y);
N=N1+N2-1;
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;
运行结果:
xy4=yx4=2
7
16
17
12
(四)、设计一个病态(矩阵)系统,分析其病态程度;找出对应的解决方法(提示:
添加白噪因子)。
程序:
A=[1988]';
B=[0427]';
N1=length(A);
N2=length(B);
N=N1+N2-1;
A0=zeros(N,1);%需要将原始数组转化为相同长度的两个数组(补零)
B0=zeros(N,1);
A0(1:
N1)=A;
B0(1:
N2)=B;
F=zeros(N);
w=exp(-sqrt(-1)*2*pi/N);%[w^(i*j)]
fori=1:
N
forj=1:
N
F(i,j)=power(w,(i-1)*(j-1));%DFT.w^((i-1)*(j-1))
end
end
FA=F*A0;
FB=F*B0;
FA
(2)=0;
FA(7)=0;
a=ifft(FA);
fori=1:
N
forj=1:
N
aa(i,j)=a(mod(i-j,N)+1);%aa为病态矩阵
end
end
x1=inv(aa)*B0;%x1为病态系数矩阵时的解
%添加白躁因子
C=F*aa*inv(F);
D=C*conj(C);
delta=0.01*eye(N);
aaa=D+delta;%添加白躁因子后的系数矩阵
X=inv(aaa)*conj(C)*FB%添加白躁因子后的解的谱
x2=ifft(X);%添加白躁因子后的解
运行结果:
X=
0.5000-0.0000i
-0.0000+0.0000i
-0.2714-0.8369i
1.1886+0.1326i
1.1886-0.1326i
-0.2714+0.8369i
0.0000+0.0000i
x1=1.0e+015*
-1.3208-0.0000i
-2.8831-0.0000i
-2.2743-0.0000i
0.0471+0.0000i
2.3330+0.0000i
2.8621+0.0000i
1.2360-0.0000i
x2=0.3335+0.0000i
-0.0006+0.0000i
0.2789-0.0000i
-0.2764-0.0000i
0.1714-0.0000i
0.4272-0.0000i
-0.4340+0.0000i
(五)、设计一个一维滤波处理程序(1、分别做低通、高通、带通、带阻等理想滤波器进行处理;2、窗函数)。
1.一维低通滤波器
程序:
a=pi/80;
x=(a*(1:
300))';
y=sin(x).*cos(x.^3);%要创建一个要被处理的信号;
b=length(x);
f=zeros(b,1);
fori=1:
b
ifi<0.2*b|i>0.8*b
f(i)=1.5;
end
end
fy=fft(y).*f;
y1=ifft(fy);
subplot(3,1,1);
plot(x,y);
text(2.5,1,'滤波前的信号');
title('一维低通滤波器')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,2)
plot(x,y1)
text(2.5,1,'滤波后的信号')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,3)
plot(f)
text(135,1.6,'频率域')
axis([0,300,0,2])
2.一维高通滤波器
程序:
a=pi/80;
x=(a*(1:
300))';
y=sin(x).*cos(x.^3);%要创建一个要被处理的信号;
b=length(x);
f=zeros(b,1);
fori=1:
b
ifi>0.2*b&i<0.8*b
f(i)=1.5;
end
end
fy=fft(y).*f;
y1=ifft(fy);
subplot(3,1,1);
plot(x,y);
text(2.5,1,'滤波前的信号');
title('一维高通滤波器')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,2)
plot(x,y1)
text(2.5,1,'滤波后的信号')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,3)
plot(f)
text(135,1.6,'频率域')
axis([0,300,0,2])
3.一维带通滤波器
程序:
a=pi/80;
x=(a*(1:
300))';
y=sin(x).*cos(x.^3);%要创建一个要被处理的信号;
b=length(x);
f=zeros(b,1);
fori=1:
b
if(i>0.15*b&i<0.35*b)|(i>0.65*b&i<0.85*b)
f(i)=1.5;
end
end
fy=fft(y).*f;
y1=ifft(fy);
subplot(3,1,1);
plot(x,y);
text(2.5,1,'滤波前的信号');
title('一维带通滤波器')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,2)
plot(x,y1)
text(2.5,1,'滤波后的信号')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,3)
plot(f)
text(135,1.6,'频率域')
axis([0,300,0,2])
4.一维带阻滤波器
程序:
d=pi/100;
x=(d*(1:
200))';
y=sin(x).*cos(x.^3);%要创建一个要被处理的信号
n=length(x);
f=zeros(n,1);
fori=1:
n
if(i>0.15*n&i<0.35*n)|(i>0.65*n&i<0.85*n)
f(i)=-1;
end
end
f=f+1
fy=fft(y).*f;
y1=ifft(fy);
subplot(3,1,1);
plot(x,y);
text(2.5,1,'滤波前的信号');
title('一维带阻滤波器')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,2)
plot(x,y1)
text(2.5,1,'滤波后的信号')
axis([0,2*pi,-1.5,1.5]);
subplot(3,1,3)
plot(f)
text(90,1.5,'频率域')
axis([0,200,0,2])
(六)、设计一个二维滤波处理程序(分别做低通、高通等处理)。
程序:
%低通滤波器
x=(peaks(85)*rand(85));
subplot(2,2,1);
surf(x)
title('滤波前信号');
X=fft(x);
subplot(2,2,2);
surf(abs(X))
title('滤波前频谱');
X(30:
70,:
)=0;
X(:
30:
70)=0;
subplot(2,2,3);
surf(abs(X))
title('滤波后频谱');
subplot(2,2,4);
Y=ifft(X);
surf(abs(Y))
title('滤波后信号');
二维低通
%高通滤波器
x=(peaks(85)*rand(85));
subplot(2,2,1);
surf(x)
title('滤波前信号');
subplot(2,2,2);
X=fft(x);
surf(abs(X))
title('滤波前频谱');
X(1:
30,:
)=0;
X(70:
100,:
)=0;
X(:
1:
30)=0;
X(:
70:
100)=0;
subplot(2,2,3);
surf(abs(X))
title('滤波后频谱');
Y=ifft(X);
subplot(2,2,4);
surf(abs(Y))
title('滤波后信号');
二维高通
(七)、验证时间域的循环褶积对应的是频率域的乘积;线性褶积则不然。
1.时间域的循环褶积对应的是频率域的乘积
程序:
%时间域求a和b的循环褶积
m=zeros(5);
a=[1,2,2,0,1];
fori=1:
5
forj=1:
5
m(i,j)=a(mod(i-j,5)+1);
end
end
b=[0,1,3,4,7]';
x=m*b
%求a和b频率的乘积对应的循环褶积
a=[1,2,2,0,1];
b=[0,1,3,4,7]';
N1=length(a);
N2=length(b);
N=5;
w2=exp(-sqrt(-1)*2*pi/N);
a0=zeros(N,1);
b0=zeros(N,1);
a0(1:
N1)=a;
b0(1:
N2)=b;
F=zeros(N);
fori=1:
N
forj=1:
N
F(i,j)=power(w2,(i-1)*(j-1));
end
end
A1=F*a0;
B1=F*b0;
y=1/N*conj(F)*(A1.*B1)%x=y,所以对于时间域的循环褶积,对应的是频率域的乘积
%时间域求a和b的线性褶积
a=[0,1,2,3]';
b=[1,2,3,4,5]';
a1=zeros(5,1);
a1(1:
4)=a;
A=a1;
N1=length(a);
N2=length(b);
N=N1+N2-1;
B=zeros(N,N2);
forj=1:
N2
fori=j:
j+N1-1
B(i,j)=a(i-j+1);
end
end
Z=B*b
%%求a和b频率的乘积对应的线性褶积
M=fft(A);
K=fft(b);
L=M.*K;
m=zeros(8,1);
m(1:
5)=L;
R=m;
T=ifft(R)%Z!
=T,所以对于时间域的线性褶积,对应的不是频率域的乘积
运行结果:
x=23
18
9
19
21
y=23.0000+0.0000i
18.0000+0.0000i
9.0000-0.0000i
19.0000-0.0000i
21.0000-0.0000i
Z=0
1
4
10
16
22
22
15
T=13.7500
11.8077-1.3464i
13.5853+2.3353i
10.9425+0.1274i
11.2500+4.0981i
6.9335-1.7880i
11.7097-0.4597i
10.0212-2.9667i
(八)、请用通俗、易懂的语言说明数字信号处理中的一种性质、一条定理或一个算例(顺便利用Matlab对其进行实现)。
解:
褶积的可交换性:
设有两个矩阵A和B,则有
证明程序:
x=[123]';
y=[234]';
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
xy=A1*y;
yx=A2*x;
运行结果:
xy1=2xy2=2
77
1616
1717
1212
xy=yx,所以得证。
六、结束语
此次信号实习基本上完成了任务,对matlab有了更深的认识,更好的掌握了matlab的使用方法。
通过此次实习,加深了信号的感性认识。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号处理 数字信号处理上机实习报告 数字信号 处理 上机 实习 报告