数字信号处理 实习大报告 22.docx

1、数字信号处理 实习大报告 22一、从给定的程序（文件包Friday.rar）中，选择一个源程序做详细标注。（目的：熟悉Matlab程序）clear;close all;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;for i=1:Dimension %开始循环,从1到401 f=f+de

2、lta_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频谱的相位谱endfigure(1); %画图1subplot(231); %画图1的第一个子图plot(f1:delta_f:f2,X(:,1),-k,LineWidth,2)

3、; % 用黑色实线画X(f)= sin(f)hold on; %保留 plot(f1 f2,0 0,:r); %用红色虚线画一条直线text(pi/2,-1.5,X(f)=sin(f),FontSize,12,FontWeight,bold); %在(pi/2,-1.5)处标注”X(f)= sin(f)” axis(f1 f2 -2 2); %确定x、y的取值范围title(1a) Continuous Spectrum,FontSize,12,FontWeight,bold); %写标题subplot(232); %画第二个子图plot(f1:delta_f:f2,X(:,2),-k,Lin

4、eWidth,2); %X 频谱的振幅图hold on; %保留plot(f1 f2,0 0,: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(f1 f2 -2 2); %确定x、y的取值范围title(1b) Amplitude Spectrum,FontSize,12,FontWeight,bold); %注标题subplot(233); %画第三个子图plot(f1 f2,0 0,:r); %用红色虚线画

5、一条直线hold on; %保留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(f1 f2 -3.5 3.5); %确定x、y的取值范围title(1c) Phase Spectrum,FontSize,12,FontWeight,bold); %注标题subplot(234); %画第四个子图plot(f1:delta_f:f2,re

6、al(Y(:,1),-k,LineWidth,2); %画出Y频谱的实部部分的图hold on; %保留plot(f1:delta_f:f2,imag(Y(:,1),-b,LineWidth,2); %画出Y频谱的虚部部分的图hold on; %保留plot(f1 f2,0 0,:r); %用红色虚线画一条直线text(pi/2,-1.5,Y(f)=exp(if),FontSize,12,FontWeight,bold); 在（pi/2,-1.5）处标注“Y(f)=exp(if)”axis(f1 f2 -2 2); %确定x、y的取值范围title(2a) Continuous Spectru

7、m,FontSize,12,FontWeight,bold); %注标题subplot(235); %画第五个子图plot(f1:delta_f:f2,Y(:,2),-k,LineWidth,2); %画出Y频谱的振幅谱hold on; %保留plot(f1 f2,0 0,: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(f1 f2 -2 2);title(2b) Amplitude Spectrum,FontSize,12,FontW

8、eight,bold);subplot(236); %画第五个子图plot(f1:delta_f:f2,Y(:,3),-k,LineWidth,2); %画出Y频谱的相位谱hold on;plot(f1 f2,0 0,:r);text(pi/2,-2.5,ArgY(f)=angle(exp(if),FontSize,12,FontWeight,bold);axis(f1 f2 -3.5 3.5);title(2c) Phase Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT0.jpeg;fi

9、gure(2); %画图2，即单独画出图1的第一个子图plot(f1:delta_f:f2,X(:,1),-k,LineWidth,2);hold on;plot(f1 f2,0 0,:r);text(1,-1,X(f)=sin(f),FontSize,12,FontWeight,bold);axis(f1 f2 -2 2);title(1a) Continuous Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT1.jpeg;figure(3); %画图3，即单独画出图1的第二个子图plot

10、(f1:delta_f:f2,X(:,2),-k,LineWidth,2);hold on;plot(f1 f2,0 0,:r);text(1,-1,|X(f)|=abs(sin(f),FontSize,12,FontWeight,bold);axis(f1 f2 -2 2);title(1b) Amplitude Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT2.jpeg;figure(4); %画图4，即单独画出图1的第三个子图plot(f1 f2,0 0,:r);hold on;plo

11、t(f1:delta_f:f2,X(:,3),-k,LineWidth,2);text(1,-2,ArgX(f)=angle(sin(f),FontSize,12,FontWeight,bold);axis(f1 f2 -3.5 3.5);title(1c) Phase Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT3.jpeg;figure(5); %画图5，即单独画出图1的第四个子图plot(f1:delta_f:f2,real(Y(:,1),-k,LineWidth,2);hold o

12、n;plot(f1:delta_f:f2,imag(Y(:,1),-b,LineWidth,2);hold on;plot(f1 f2,0 0,:r);text(0.2,-1,Y(f)=exp(if),FontSize,12,FontWeight,bold);axis(f1 f2 -2 2);title(2a) Continuous Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT4.jpeg;figure(6); %画图6，即单独画出图1的第五个子图plot(f1:delta_f:f2,Y(

13、:,2),-k,LineWidth,2);hold on;plot(f1 f2,0 0,:r);text(1,-1,|Y(f)|=abs(exp(if),FontSize,12,FontWeight,bold);axis(f1 f2 -2 2);title(2b) Amplitude Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT5.jpeg;figure(7); %画图7，即单独画出图1的第六个子图plot(f1:delta_f:f2,Y(:,3),-k,LineWidth,2);hold

14、 on;plot(f1 f2,0 0,:r);text(5.5,-2,ArgY(f)=angle(exp(if),FontSize,12,FontWeight,bold);axis(f1 f2 -3.5 3.5);title(2c) Phase Spectrum,FontSize,12,FontWeight,bold);print -djpeg -r0 Spectrum_Showing_CFT6.jpeg;二、能够利用Matlab熟悉地画图，内容包括：X、Y坐标轴上的label，每幅图上的title，绘画多条曲线时的legend，对图形进行适当的标注等。（1）在一副图上画出多幅小图(程序：pr

15、oblem2.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的函数hold on; %保留plot(x,cot(4*x),b-.); %画cot(4*x)的函数hold on; %保留plot(x,sec(x),ro); %画sec(x)的函数title(y1=tan(x),y2=cot(x),y3=sec(x); %为此图命名xlab

16、el(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=1 3 5 7，y=2 4 6 时，利用matlab中conv命令求得的线性褶积z1=z2=2 10 28 52 58 42.根据线性褶积公式： 得出的结果为

17、z1=z2=2 10 28 52 58 42程序如下：x=1 3 5 7;y=2 4 6 ;%利用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);for j=1:N2 for i=j:j+N1-1 F(i,j)=x(i-j+1); endendZ1=F*y; %利用公式编写程序计算出来的线性褶积结果 运行结果：Z1 = 2 10 28 52 58 42 (2

18、)在时间域正演计算循环褶积(程序shixun2.m)X=1 3 5 7;Y=2 4 6 0;N=length(X); F=zeros(N); AA=zeros(N);DD=zeros(N);w=exp(-sqrt(-1)*2*pi/N);for i=1:N for j=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); endendZ=AA*Y运行结果为：Z = 60 52 2852(3)在频率域反演计算循环褶积(程序pinxun3.m)x=1 3 5 7;y=2 4 6 0;XX=

19、fft(x); %信号x的频谱YY=fft(y); %信号y的频谱XY3=XX.*YY; %两信号频谱乘积Z1=ifft(XY3) %由两频谱相乘再反变换得到的信号运行程序如下：Z1 = 60 52 28 52(4)验证循环褶积计算时所存在的边界效应现象(程序pinxun4.m)程序如下：clear;close all;x=1 3 5 7; y=2 4 6; 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

20、); Y0=fft(y0); xy2=ifft(X0.*Y0); 运行结果： z1 = 2 10 28 52 5842xy2 = 2.0000 10.0000 28.0000 52.0000 58.0000 42.0000由以上可知，当N=N1+N2-1时，z1-xy2,即循环褶积等于线性褶积由（1）（2）可知，NN1+N2-1时，循环褶积不等于线性褶积所以循环褶积的维度在length(x)+length(y)-1处存在边界效应，大于等于该值时，循环褶积等于线性褶积，小于该维度时二者不等。（5）编写一个做相关分析的源程序(程序xiangguan.m)x=1 3 5 7; y=2 4 6;%对这

21、两个向量作相关分析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); for j=1:N for i=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=

22、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);for i=1:N; for j=1:N; F(i,j)=power(W,(i-1)*(j-1); AA(i,j)=A(mod(i-j,N)+1); endendFx=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.

23、0483 0.1769 0.4638 0.1038 0.4546五、设计一个一维滤波处理程序(1、分别做低通、高通、带通、带阻等(程序process.m)程序运行结果：对于简单函数分别作四种滤波器后：六、设计一个二维滤波处理程序（分别做低通、高通等处理）。运行程序为：原始信号频谱：（程序problem6.1.m） 二维低通滤波：（程序problem6.2.m）二维高通（程序problem6.3.m） 七、验证时间域的循环褶积对应的是频率域的乘积；线性褶积则不然。线性褶积的证明：(程序xianxing.m)x=1 3 5 7;y=2 4 6 ;N1=length(x);N2=length(y);

24、N=N1+N2-1;F=zeros(N,N2);for j=1:N2 for i=j:j+N1-1 F(i,j)=x(i-j+1); endendZ1=F*y;运行结果：Z1 = 2 10 28 52 58 42M =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 32z1不等于z2由此可知时间域的线性褶积对应的不是频率域的乘积。 循环褶积的证明（程序xunhuan.m）x=1 3 5 7;y=2 4 6 0;N=len

25、gth(x);F=zeros(N); AA=zeros(N);DD=zeros(N);w=exp(-sqrt(-1)*2*pi/N);for i=1:N for j=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); endendZ1=AA*yN1 =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 52Z2 = 60 52 2852Z1等于z2由此可知时间域的循环褶积对应的是频率域的乘积。八、请用通俗、易懂的语言说明数字信号处理中的一种性质、一条定理或一个算例（顺便利用Matlab对其进行实现）。(程序fanzhuan.m)频谱的有些性质在我们的日常生活中也可以观察到，例如频谱的翻转定理，当我们看时间时，比如是8:10，若我们通过镜子观察时，这时我们看到的时间是4:50，假设6和12的连线为y轴，则时针和分针都关于y轴进行了翻转。所得图形如下：