数字信号处理 实习大报告 22.docx
《数字信号处理 实习大报告 22.docx》由会员分享，可在线阅读，更多相关《数字信号处理 实习大报告 22.docx（24页珍藏版）》请在冰豆网上搜索。
数字信号处理实习大报告22
一、从给定的程序（文件包Friday.rar）中，选择一个源程序做详细标注。
（目的：
熟悉Matlab程序）
clear;
closeall;
Dimension=401;%赋值
X=zeros（Dimension,3）;%设置X矩阵,401行3列,所有元素为零
Y=zeros（Dimension,3）;%设置Y矩阵,401行3列,所有元素为零
f1=0;%赋值
f2=4*pi;%赋值
delta_f=（f2-f1）/（Dimension-1）;
ii=sqrt（-1）;
f=f1-delta_f;
fori=1:
Dimension%开始循环,从1到401
f=f+delta_f;%f以delta_f依次增加
X（i,1）=sin（f）;%由不同f值得出的sin（f）的值放入X矩阵的第一列;是频谱X
X（i,2）=abs（X（i,1））;%X频谱的绝对值即振幅
X（i,3）=angle（X（i,1））;%X频谱的相位谱
Y（i,1）=exp（f*ii）;%求另一个频谱Y
Y（i,2）=abs（Y（i,1））;%Y频谱的绝对值即振幅
Y（i,3）=angle（Y（i,1））;%Y频谱的相位谱
end
figure
（1）;%画图1
subplot（231）;%画图1的第一个子图
plot（f1:
delta_f:
f2,X（:
1）,'-k','LineWidth',2）;%用黑色实线画X（f）=sin（f）
holdon;%保留
plot（[f1f2],[00],':
r'）;%用红色虚线画一条直线
text（pi/2,-1.5,'X（f）=sin（f）','FontSize',12,'FontWeight','bold'）;%在（pi/2,-1.5）处标注”X（f）=sin（f）”
axis（[f1f2-22]）;%确定x、y的取值范围
title（'（1a）ContinuousSpectrum','FontSize',12,'FontWeight','bold'）;%写标题
subplot（232）;%画第二个子图
plot（f1:
delta_f:
f2,X（:
2）,'-k','LineWidth',2）;%X频谱的振幅图
holdon;%保留
plot（[f1f2],[00],':
r'）;%用红色虚线画一条直线
text（pi/2,-1.5,'|X（f）|=abs（sin（f））','FontSize',12,'FontWeight','bold'）;%在（pi/2,-1.5）处标注|X（f）|=abs（sin（f））
axis（[f1f2-22]）;%确定x、y的取值范围
title（'（1b）AmplitudeSpectrum','FontSize',12,'FontWeight','bold'）;%注标题
subplot（233）;%画第三个子图
plot（[f1f2],[00],':
r'）;%用红色虚线画一条直线
holdon;%保留
plot（f1:
delta_f:
f2,X（:
3）,'-k','LineWidth',2）;%画相位谱
text（pi/2,-2.5,'ArgX（f）=angle（sin（f））','FontSize',12,'FontWeight','bold'）;%在（pi/2,-2.5）处注“ArgX（f）=angle（sin（f））“
axis（[f1f2-3.53.5]）;%确定x、y的取值范围
title（'（1c）PhaseSpectrum','FontSize',12,'FontWeight','bold'）;%注标题
subplot（234）;%画第四个子图
plot（f1:
delta_f:
f2,real（Y（:
1））,'-k','LineWidth',2）;%画出Y频谱的实部部分的图
holdon;%保留
plot（f1:
delta_f:
f2,imag（Y（:
1））,'--b','LineWidth',2）;%画出Y频谱的虚部部分的图
holdon;%保留
plot（[f1f2],[00],':
r'）;%用红色虚线画一条直线
text（pi/2,-1.5,'Y（f）=exp（if）','FontSize',12,'FontWeight','bold'）;在（pi/2,-1.5）处标注“Y（f）=exp（if）”
axis（[f1f2-22]）;%确定x、y的取值范围
title（'（2a）ContinuousSpectrum','FontSize',12,'FontWeight','bold'）;%注标题
subplot（235）;%画第五个子图
plot（f1:
delta_f:
f2,Y（:
2）,'-k','LineWidth',2）;%画出Y频谱的振幅谱
holdon;%保留
plot（[f1f2],[00],':
r'）;
text（pi/2,-1.5,'|Y（f）|=abs（exp（if））','FontSize',12,'FontWeight','bold'）;在（pi/2,-1.5）标注|Y（f）|=abs（exp（if））
axis（[f1f2-22]）;
title（'（2b）AmplitudeSpectrum','FontSize',12,'FontWeight','bold'）;
subplot（236）;%画第五个子图
plot（f1:
delta_f:
f2,Y（:
3）,'-k','LineWidth',2）;%画出Y频谱的相位谱
holdon;
plot（[f1f2],[00],':
r'）;
text（pi/2,-2.5,'ArgY（f）=angle（exp（if））','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-3.53.5]）;
title（'（2c）PhaseSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT0.jpeg;
figure
（2）;%画图2，即单独画出图1的第一个子图
plot（f1:
delta_f:
f2,X（:
1）,'-k','LineWidth',2）;
holdon;
plot（[f1f2],[00],':
r'）;
text（1,-1,'X（f）=sin（f）','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-22]）;
title（'（1a）ContinuousSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT1.jpeg;
figure（3）;%画图3，即单独画出图1的第二个子图
plot（f1:
delta_f:
f2,X（:
2）,'-k','LineWidth',2）;
holdon;
plot（[f1f2],[00],':
r'）;
text（1,-1,'|X（f）|=abs（sin（f））','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-22]）;
title（'（1b）AmplitudeSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT2.jpeg;
figure（4）;%画图4，即单独画出图1的第三个子图
plot（[f1f2],[00],':
r'）;
holdon;
plot（f1:
delta_f:
f2,X（:
3）,'-k','LineWidth',2）;
text（1,-2,'ArgX（f）=angle（sin（f））','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-3.53.5]）;
title（'（1c）PhaseSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT3.jpeg;
figure（5）;%画图5，即单独画出图1的第四个子图
plot（f1:
delta_f:
f2,real（Y（:
1））,'-k','LineWidth',2）;
holdon;
plot（f1:
delta_f:
f2,imag（Y（:
1））,'--b','LineWidth',2）;
holdon;
plot（[f1f2],[00],':
r'）;
text（0.2,-1,'Y（f）=exp（if）','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-22]）;
title（'（2a）ContinuousSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT4.jpeg;
figure（6）;%画图6，即单独画出图1的第五个子图
plot（f1:
delta_f:
f2,Y（:
2）,'-k','LineWidth',2）;
holdon;
plot（[f1f2],[00],':
r'）;
text（1,-1,'|Y（f）|=abs（exp（if））','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-22]）;
title（'（2b）AmplitudeSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT5.jpeg;
figure（7）;%画图7，即单独画出图1的第六个子图
plot（f1:
delta_f:
f2,Y（:
3）,'-k','LineWidth',2）;
holdon;
plot（[f1f2],[00],':
r'）;
text（5.5,-2,'ArgY（f）=angle（exp（if））','FontSize',12,'FontWeight','bold'）;
axis（[f1f2-3.53.5]）;
title（'（2c）PhaseSpectrum','FontSize',12,'FontWeight','bold'）;
print-djpeg-r0Spectrum_Showing_CFT6.jpeg;
二、能够利用Matlab熟悉地画图，内容包括：
X、Y坐标轴上的label，每幅图上的title，绘画多条曲线时的legend，对图形进行适当的标注等。
（1）在一副图上画出多幅小图（程序：
problem2.1.m）
（2）画出一组二维图形（程序：
problem2.2.m）
（3）画出一组三维图形（程序：
problem2.3.m）
（4）画出复数的实部与虚部（程序：
problem2.4.m）
（5）对所做的一个源程序注释：
x=0:
0.1:
2*pi;%确定x的取值范围
plot（x,tan（x）,'bp'）;%画出tanx的函数
holdon;%保留
plot（x,cot（4*x）,'b-.'）;%画cot（4*x）的函数
holdon;%保留
plot（x,sec（x）,'ro'）;%画sec（x）的函数
title（'y1=tan（x）,y2=cot（x）,y3=sec（x）'）;%为此图命名
xlabel（'x'）;%标注x轴坐标
ylabel（'y'）;%标注x轴坐标
legend（'tan（x）','cot（x）','sec（x）',4）;
text（pi,sin（pi）,'x=pi'）;在（pi,sin（pi））处标注“x=pi”
三、计算普通褶积与循环褶积，分别使用时间域与频率域两种方法进行正、反演计算，指出循环褶积计算时所存在的边界效应现象；编写一个做相关分析的源程序。
一
（1）在时间域计算线性褶积（程序shixian1.m）
当两信号分别为x=[1357]'，y=[246]'时，利用matlab中conv命令求得的线性褶积z1=z2=[21028525842]'.
根据线性褶积公式：
得出的结果为z1=z2=[[21028525842]'
程序如下：
x=[1357]';
y=[246]';
%%%%%%%%%%%%%%%%%%%%%%%%利用conv求线性褶积
Z1=conv（x,y）;
Z2=conv（y,x）;
Z1-Z2;%利用matlab中conv命令验证线性褶积具有交换性
%%%%%%%%%%%%%%%%%%%线性褶积：
利用自己编写的程序进行计算。
N1=length（x）;
N2=length（y）;
N=N1+N2-1;
F=zeros（N,N2）;
forj=1:
N2
fori=j:
j+N1-1
F（i,j）=x（i-j+1）;
end
end
Z1=F*y;%利用公式编写程序计算出来的线性褶积结果
运行结果：
Z1=
2
10
28
52
58
42
（2）在时间域正演计算循环褶积（程序shixun2.m）
X=[1357]';
Y=[2460]';
N=length（X）;
F=zeros（N）;
AA=zeros（N）;
DD=zeros（N）;
w=exp（-sqrt（-1）*2*pi/N）;
fori=1:
N
forj=1:
N
AA（i,j）=X（mod（i-j,N）+1）;
DD（i,j）=Y（mod（j-i,N）+1）;
F（i,j）=power（w,（i-1）*（j-1））;
end
end
Z=AA*Y
运行结果为：
Z=
60
52
28
52
（3）在频率域反演计算循环褶积（程序pinxun3.m）
x=[1357]';
y=[2460]';
XX=fft（x）;%信号x的频谱
YY=fft（y）;%信号y的频谱
XY3=XX.*YY;%两信号频谱乘积
Z1=ifft（XY3）%由两频谱相乘再反变换得到的信号
运行程序如下：
Z1=
60
52
28
52
（4）验证循环褶积计算时所存在的边界效应现象（程序pinxun4.m）
程序如下：
clear;
closeall;
x=[1357]';
y=[246]';
z1=conv（x,y）;
N1=length（x）;
N2=length（y）;
N=N1+N2-1;
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）;
运行结果：
z1=
2
10
28
52
58
42
xy2=
2.0000
10.0000
28.0000
52.0000
58.0000
42.0000
由以上可知，当N》=N1+N2-1时，z1-xy2,即循环褶积等于线性褶积
由
（1）
（2）可知，N《N1+N2-1时，循环褶积不等于线性褶积
所以循环褶积的维度在length（x）+length（y）-1处存在边界效应，大于等于该值时，循环褶积等于线性褶积，小于该维度时二者不等。
（5）编写一个做相关分析的源程序（程序xiangguan.m）
x=[1357]';
y=[246]';
%对这两个向量作相关分析
N1=length（x）;
N2=length（y）;
N=N1+N2-1;
x0=zeros（N,1）;
y0=zeros（N,1）;
x0（1:
N1）=x;
y0（1:
N2）=y;
A1=zeros（N）;
forj=1:
N
fori=1:
N
if（i+j）>（N+1）
A1（i,j）=x0（i+j-N-1）;
else
A1（i,j）=x0（i+j-1）;
end
end
end
xy=A1*conj（y0）
运行结果：
xy=
44
68
38
14
6
四、设计一个病态（矩阵）系统，分析其病态程度；找出对应的解决方法（提示：
添加白噪因子）。
（程序biantai.m）
x=[2,0,0,7,0,1]';
y=[3,1,4,0,1,2]';
N=length（x）;
Fx=fft（x）;
Fx
（2）=0;
Fx（N）=0;
A=ifft（Fx）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W=exp（-2*pi*i/N）;
AA=zeros（N）;
F=zeros（N）;
fori=1:
N;
forj=1:
N;
F（i,j）=power（W,（i-1）*（j-1））;
AA（i,j）=A（mod（i-j,N）+1）;
end
end
Fx=fft（x）;
Fy=fft（y）;
FX=conj（Fx）.*Fy./（conj（Fx）.*Fx+0.001）;
X=zeros（N,1）;
X=ifft（FX）;
运行结果：
X=
-0.1473
0.0483
0.1769
0.4638
0.1038
0.4546
五、设计一个一维滤波处理程序（1、分别做低通、高通、带通、带阻等（程序process.m）
程序运行结果：
对于简单函数分别作四种滤波器后：
六、设计一个二维滤波处理程序（分别做低通、高通等处理）。
运行程序为：
原始信号频谱：
（程序problem6.1.m）
二维低通滤波：
（程序problem6.2.m）
二维高通（程序problem6.3.m）
七、验证时间域的循环褶积对应的是频率域的乘积；线性褶积则不然。
线性褶积的证明：
（程序xianxing.m）
x=[1357]';
y=[246]';
N1=length（x）;
N2=length（y）;
N=N1+N2-1;
F=zeros（N,N2）;
forj=1:
N2
fori=j:
j+N1-1
F（i,j）=x（i-j+1）;
end
end
Z1=F*y;
运行结果：
Z1=
2
10
28
52
58
42
M=max（N1,N2）;
x0=zeros（M,1）;
y0=zeros（M,1）;
x0（1:
N1）=x;
y0（1:
N2）=y;
X=fft（x0）;
Y=fft（y0）;
Z2=X.*Y
运行结果：
Z2=
192
32
-16
32
z1不等于z2
由此可知时间域的线性褶积对应的不是频率域的乘积。
循环褶积的证明（程序xunhuan.m）
x=[1357]';
y=[2460]';
N=length（x）;
F=zeros（N）;
AA=zeros（N）;
DD=zeros（N）;
w=exp（-sqrt（-1）*2*pi/N）;
fori=1:
N
forj=1:
N
AA（i,j）=x（mod（i-j,N）+1）;
DD（i,j）=y（mod（j-i,N）+1）;
F（i,j）=power（w,（i-1）*（j-1））;
end
end
Z1=AA*y
N1=length（x）;
N2=length（y）;
M=max（N1,N2）;
x0=zeros（N,1）;
y0=zeros（N,1）;
x0（1:
N1）=x;
y0（1:
N2）=y;
X=fft（x0）;
Y=fft（y0）;
Z2=ifft（X.*Y）
运行结果如下：
Z1=
60
52
28
52
Z2=
60
52
28
52
Z1等于z2
由此可知时间域的循环褶积对应的是频率域的乘积。
八、请用通俗、易懂的语言说明数字信号处理中的一种性质、一条定理或一个算例（顺便利用Matlab对其进行实现）。
（程序fanzhuan.m）
频谱的有些性质在我们的日常生活中也可以观察到，例如频谱的翻转定理，当我们看时间时，比如是8:
10，若我们通过镜子观察时，这时我们看到的时间是4:
50，假设6和12的连线为y轴，则时针和分针都关于y轴进行了翻转。
所得图形如下：
升级会员 