《数字信号处理》上机实习报告 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的使用方法。
通过此次实习，加深了信号的感性认识。
升级会员 