数字信号处理实验报告.docx

1、数字信号处理实验报告数字信号处理实验讲义实验一 MATLAB简介实验目的 1熟悉MATLAB软件的使用方法； 2MATLAB的绘图功能； 3用MATLAB语句实现信号的描述及变换。实验原理 1在MATLAB下编辑和运行程序在MATLAB中，对于简单问题可以在命令窗（command windows）直接输入命令，得到结果；对于比较复杂的问题则可以将多个命令放在一个脚本文件中，这个脚本文件是以m为扩展名的，所以称之为M文件。用M文件进行程序的编辑和运行步骤如下：（1）打开MATLAB，进入其基本界面；（2）在菜单栏的File项中选择新建一个M文件；（3）在M文件编辑窗口编写程序；（4）完成之后，可

2、以在编辑窗口利用Debug工具调试运行程序，在命令窗口查看输出结果；也可以将此文件保存在某个目录中，在MATLAB的基本窗口中的File项中选择Run The Script，然后选择你所要运行的脚本文件及其路径，即可得出结果；也可以将此文件保存在当前目录中，在MATLAB命令窗口，“”提示符后直接输入文件名。 2MATLAB的绘图功能 plot(x,y) 基本绘图函数，绘制 x和y之间的坐标图。 figure(n ) 开设一个图形窗口n subplot(m,n,N) 分割图形窗口的MATLAB函数，用于在一个窗口中显示多个图形，将图形窗口分为m行n列，在第N个窗口内绘制图形。axis(a0,b

3、0,a1,b1 ) 调整坐标轴状态 title( ) 给图形加题注xlabel ( ) 给x轴加标注ylabel ( ) 给y轴加标注grid 给图形加网格线 3信号描述及变换 信号描述及变换包括连续时间信号和离散时间信号内容，详细内容请见课本第1章、第2章。实验内容1试用MATLAB绘制出下列信号的波形：(Signal 1.6)（1） ；（2） （3） ；（4） ；（5） 【程序代码】clear all;close all;clc;syms t;x1=exp(-1.5*t)x2=3*sin(0.5*pi*t)x3=0.5+0.5*sym(sign(t)x4=sym(heaviside(t)+

4、sym(heaviside(t-1)-sym(2*heaviside(t-2)x5=.5*t*(sym(heaviside(t)-sym(heaviside(t-4)subplot(2,3,1);ezplot(x1);axis(-6 3 -500 7000);title(x1(t)=exp(-1.5t);grid onsubplot(2,3,2);ezplot(x2);title(x2(t)=3sin(0.5t);grid onsubplot(2,3,3);fplot(sign(t)/2 + 1/2,-10 10,1e-8);ezplot(x3,-10 10);axis(-10 10 -.2

5、1.2);xlabel(t);title(x3(t)=0.5+0.5sgn(t);grid onsubplot(2,3,4);ezplot(x4,-1 3);title(x4(t)=u(t)+u(t-1)-2u(t-2);grid onsubplot(2,3,5);ezplot(x5,-2 6);title(x5(t)=0.5tu(t)-u(t-4);grid onsubplot(2,3,6);axis off2 已知连续时间信号(Signal 1.7)，试用MATLAB绘制出下列信号的波形：（1） ；（2） ；（3） ；（4） ；（5） 。【程序代码】clear all;close all;

6、clc;figure(2)syms t;x1=(4-t)*(sym(heaviside(t)-sym(heaviside(t-4)x2=exp(-2*t)*sym(heaviside(t)x3=sin(2*pi*t)x4=subs(x1,t,t/2)x5=subs(x4,t,t-2)x6=subs(x2,t,-t)x7=x2+x6x8=x7*x3subplot(2,3,1)text(0,0.9,x1(t)=(4-t)u(t)-u(t-4) );text(0,0.7,x2(t)=exp(-2t)u(t);text(0,0.5,x3(t)=sin(2t);axis off;box off;subp

7、lot(2,3,2);ezplot(x5,-3,10);gridtitle(x4(t)=x1(t/2);subplot(2,3,3);ezplot(x5,-1,12);gridtitle(x5(t)=x4(t-2);subplot(2,3,4);ezplot(x6,-10,5);axis(-5 5 -0.1 1.1);gridtitle( x6(t)=x2(-t);subplot(2,3,5);ezplot(x7,-5 5);axis(-5 5 -0.1 1.1);grid;title(x7(t)=x2(t)+x6(t);subplot(2,3,6);tv8=-2.5:0.05:2.5;xv8

8、=subs(x8,tv8);plot(tv8,xv8)axis(-2.5 2.5 -1 1);gridxlabel(t);title(x8(t)=x7(t)x3(t);clear tv8 xv83列出单位冲激信号、单位阶跃信号、正弦信号的MATLAB表达式，并绘出信号波形。【程序代码】clear all;close all;clcsyms t;x1=sym(dirac(t);x2=sym(Heaviside(t);x3=sin(t);tn=-6.3:0.1:6.3;xn1=subs(x1,t,tn);xn2=subs(x2,t,tn);xn3=subs(x3,t,tn);plot(tn,xn1

9、,k,tn,xn2,r,tn,xn3,m);gridxlabel(t);ylabel(x(t);legend(dirac,Heaviside,sin)hold onplot(tn,xn1,k.,tn,xn2,r.,tn,xn3,m.)实验二 用FFT实现信号的谱分析实验目的 1了解FFT在信号谱分析中的作用； 2了解谱分析的一般步骤和方法。实验原理关于信号谱分析的步骤和方法参见教材第3章相关内容。为了解信号的特点，了解信号频谱分布情况，应该对信号进行谱分析，计算出信号的幅度谱、相位谱和功率谱。信号的谱分析可以用FFT实现，讨论如下：1 谱分析中的参数选择；A 若已知信号的最高频率，为防止混叠，

10、选定采样频率： (1)B 根据实际需要，选定频率分辨，一但选定后，即可确定FFT所需的点数N (2)我们希望越小越好，但越小，N越大，计算量、存储量也随之增大。一般取N为2的整次幂，以便用FFT计算，若已给定N，可用补零方法便N为2的整次幂。C 和N确定后，即可确定所需相应模拟信号的长度 (3)分辨率反比于T，而不是N，在给定的T的情况下，靠减小来增加N是不能提高分辨率的，因为为常数 2谱分析步骤； A 数据准备 (4) B 使用FFT计算信号的频谱 (5) (6) C 由频谱计算幅度谱、相位谱和功率谱 (7) (8) (9)3实验中用到的一些基本函数简介y=fft(x,n) ； 计算n点的F

11、FT。abs(x) ； 取绝对值。angle(z) ； 取相角。Pxx, f= periodogram (xn, nfft, fs, window) ；%周期图谱估计Pxx, f=pwelch (xn, nfft, fs, window, noverlap)；%平均周期图法 Pxx=psd (xn) ；功率谱密度 实验内容 1. 已知序列x(n)=2sin(0.48n)+cos(0.52n) 0n100,试绘制x(n)及它的频谱图。若x(n)=sin(0.56n)+2cos(0.25n),结果又如何？【程序代码】clear all;close all;clcN=100;n=0:N-1;xn=2

12、*sin(0.48*pi*n)+cos(0.52*pi*n);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);subplot(1,2,1)plot(n,xn)xlabel(n);ylabel(x(n);title(x(n)N=100);subplot(1,2,2)k=0:length(magXK)-1;stem(k,magXK,.);xlabel(k);ylabel(|X(K)|);title(X(K)N=100); 2. 对下面信号进行频谱分析，求幅度谱和相位谱。 (1)， (2)，【程序代码】clear all;close all;clcfs1=500

13、0;dt1=1/fs1;N1=0.004/dt1;n1=0:N1-1;xn1=0.8.(n1*dt1);Xk1=fft(xn1,N1);mag1=2*abs(Xk1)/N1;pha1=angle(Xk1);f1=n1*fs1/N1;subplot(221)stem(f1,mag1,fill);title(Magnitude 1)xlabel(f);ylabel(|Xk(1)|);subplot(222)stem(f1,pha1,fill);title(Phase 1)xlabel(f);ylabel(|Phk(1)|);fs2=128;dt2=1/fs2;N2=18;n2=0:N2-1;xn2

14、=sin(n2*dt2)./(n2*dt2+eps);Xk2=fft(xn2,N2);mag2=2*abs(Xk2)/N2;pha2=angle(Xk2);f2=n2*fs2/N2;subplot(223)stem(f2,mag2,fill);title(Magnitude 2)xlabel(f);ylabel(|Xk(2)|);subplot(224)stem(f2,pha2,fill);title(Phase 2)xlabel(f);ylabel(|Phk(2)|) 3. 给定信号，现在对采样，采样点数，采样频率=50Hz，设采样序列为，编写程序计算的频谱，并绘图；改变采样频率，得到序列，

15、计算的频谱，并绘图；增大采样点数，得到序列，计算的频谱，并绘图；采样点数N=64，采样频率=300Hz，在采样点后补零得到新序列，计算的频谱，并绘图。(Signal 3.18)【程序代码】clear all;close all; clcN=16;n=0:N-1;f1=15;f2=18;fs=50;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(1);subplot(121)plot(f,mag);gridtitle(Magnit

16、ude);xlabel(f);ylabel(|Xk|);subplot(122)plot(f,pha);gridtitle(Phase);xlabel(f);ylabel(?Xk?);N=16;n=0:N-1;f1=15;f2=18;fs=36; dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(2);subplot(221)plot(f,mag,r);gridtitle(Magnitude(36Hz Sampling reque

17、ncy);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,r);gridtitle(Phase(36Hz Sampling requency);xlabel(f);ylabel(?Xk?);fs=100; dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,r);gridtitle(Magnitude(100Hz Sampling Frequency);xl

18、abel(f);ylabel(|Xk|);subplot(224)plot(f,pha,r);gridtitle(Phase(100Hz Sampling Frequency);xlabel(f);ylabel(?Xk?);N=64;n=0:N-1;f1=15;f2=18;fs=50;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(3);subplot(221)plot(f,mag,m);gridtitle(Magnitude(

19、64 Sampling Points);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,m);gridtitle(Phase(64 Sampling Points);xlabel(f);ylabel(?Xk?);N=512; n=0:N-1;dt=1/fs;xn=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,m);gridtitle(Magnitude(512 Sampli

20、ng Points);xlabel(f);ylabel(|Xk|);subplot(224)plot(f,pha,m);gridtitle(Phase(512 Sampling Points);xlabel(f);ylabel(?Xk?);N=64;n=0:N-1;f1=15;f2=18;fs=300;dt=1/fs;xn0=sin(2*pi*f1*n*dt)+2*sin(2*pi*f2*n*dt);xn=xn0,zeros(1,64);Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;figure(4);subplot(221)plot(

21、f,mag,g);gridtitle(Magnitude(filled with 64 zeros);xlabel(f);ylabel(|Xk|);subplot(222)plot(f,pha,g);gridtitle(Phase (filled with 64 zeros);xlabel(f);ylabel(?Xk?);xn=xn0,zeros(1,448); Xk=fft(xn,N);mag=2*abs(Xk)/N;pha=angle(Xk);f=fs/N*n;subplot(223)plot(f,mag,g);gridtitle(Magnitude(filled with 448 zer

22、os);xlabel(f);ylabel(|Xk|);4. 试求下列差分方程所描述的输出序列的功率谱并作图。(a)(b)(c)式中，是方差为(例如，=1/12)的白噪声。【程序代码】clear all;close all;clcnfft=512;window=hanning(256);noverlap=128;w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=-0.81*x(i-2)+w(i)-w(i-1);endendPxx f=psd(x,nfft,2,window,noverlap,mean);figure(1)plot

23、(f,10*log10(Pxx);gridxlabel(f);ylabel(Power Spectrum,dB);title(a): x(n)=-0.81x(n-2)+w(n)-w(n-1)w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=w(i)-w(i-2);endendPxx f=psd(x,nfft,2,window,noverlap,mean);figure(2)plot(f,10*log10(Pxx),m);gridxlabel(f);ylabel(Power Spectrum,dB);title(b): x(n)

24、=w(n)-w(n-2)w=randn(1,1000);for i=1:1000if i=2x(i)=randn(1,1);elsex(i)=-0.81*x(i-2)+w(2);endendPxx f=psd(x,nfft,2,window,noverlap);figure(3)plot(f,10*log10(Pxx),r);gridxlabel(f);ylabel(Power Spectrum,dB);title(c): x(n)=-0.81x(n-2)+w(n)5. 一序列是由两个频率相距为的模拟信号采样得来的，即 n=0,1,15已知序列长度N=16，试采用周期图法，应用DFT分别计算当

25、=0.06及=0.01时的功率谱估计，并通过作图说明从功率谱估计的分布是否能分辨出这两个正弦信号的真实频谱？若N=64又有什么变化？【程序代码】clear all;close all;clcnfft=512;window=hanning(256);noverlap=128;df=0.06 0.01;N=16;n=0:N-1;x1=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(1)*n);x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2)*n);Pxx1 f1=psd(x1,nfft,1,window,noverlap,mean);P

26、xx2 f2=psd(x2,nfft,1,window,noverlap,mean);figure(1)plot(f1,Pxx1);gridxlabel(f);ylabel(Power Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.06)n N=16)figure(2)plot(f2,Pxx2);gridxlabel(f);ylabel(Power Spectrum );title(sin2?(0.135)n+cos2?(0.135+0.01)n N=16)N=64;n=0:N-1;x1=sin(2*pi*0.135*n)+cos(2*pi*(0.13

27、5+df(1)*n);x2=sin(2*pi*0.135*n)+cos(2*pi*(0.135+df(2)*n);Pxx1 f1=psd(x1,nfft,1,window,noverlap,mean);Pxx2 f2=psd(x2,nfft,1,window,noverlap,mean);figure(3)plot(f1,Pxx1,m);gridxlabel(f);ylabel(Power Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.06)n N=64)figure(4);plot(f2,Pxx2,m);gridxlabel(f);ylabel(Po

28、wer Spectrum);title(sin2?(0.135)n+cos2?(0.135+0.01)n N=64)6. 用MATLAB产生256点白噪声序列，应用Welch法估计其功率谱，每段长64点，重叠32点，输出平均后的功率谱图以及对256点一次求周期图的功率谱图。【程序代码】N=256;x=rand(1,N);window=hanning(64);noverlap=32;nfft=128;Pxx1 f1=periodogram(x);figure(1)plot(f1/pi,10*log10(Pxx1);gridxlabel(f);ylabel(Power Spectrum,dB);title(256)Pxx2 f2=psd(x,nfft,2,window,noverlap,line);figure(2)plot(f2,10*log10(Pxx2),