数字信号处理MATLAB上机实验教学讲义

1,Matlab基础介绍,、Matlab简介、Matlab的安装与启动、Matlab编程基础、Matlab在数字信号处理课程中的应用,2,一、Matlab简介,Matlab名字是由Matrix和Laboratory两个词的前三个字母组合而成的,是由MathWorks公司于1982年推出的一套高性能的数值计算和可视化数学软件。到20世纪90年代，MATLAB已经成为国际控制界公认的标准计算软件.Matlab编程运算与人进行科学计算的思路和表达方式完全一致。不象学习其它高级语言如Basic、Fortran和C等那样难于掌握。其基本元素是矩阵。具有强大的数值计算和图示能力。具有丰富的工具箱（ToolBox)。具有强大的与外部程序如C的交互接口。在高校中，Matlab成为数学，信息，控制等诸多学科有关课程的有效教学工具。,3,二、Matlab的安装,4,Matlab的启动,5,操作桌面的缺省外貌,6,通用操作界面指令窗（CommandWindow）历史指令窗（CommandHistory）当前目录浏览器（CurrentDirectoryBrowser）工作空间浏览器（WorkspaceBrowser）内存数组编辑器（ArrayEditor）M文件编辑器/调试器（Editor/Debugger）帮助导航/浏览器（HelpNavigator/Browser）,7,三、Matlab编程基础,1变量2数组，向量与矩阵的创建和访问3部分特殊变量和常数4部分常用运算符5Matlab的程序设计6m文件的创建，保存，运行7基本数学函数和基本绘图函数8编程举例,8,3.1变量,变量无需定义即可使用。变量名的大小写是敏感的。变量的第一个字符必须为英文字母，而且不超过31个字符。变量名可以包含下连字符、数字，但不能为空格符、标点。,9,局部（Local）变量：它存在于函数空间内部的中间变量，产生于该函数的运行过程中，其影响范围也仅限于该函数本身。全局（Global）变量：通过global指令，MATLAB也允许几个不同的函数空间以及基本工作空间共享同一个变量。,10,11,最简单的计算器使用法【例】求的算术运算结果。（1）用键盘在MATLAB指令窗中输入以下内容(12+2*(7-4)/32（2）在上述表达式输入完成后，按【Enter】键，该指令被执行。（3）在指令执行后，MATLAB指令窗中将显示以下结果。ans=2输入help*可查看该命令的帮助信息。,12,3.2数组，向量与矩阵的创建和访问,在Matlab中,这三个概念在创建和显示的时候没有任何区别。向量的创建【例】要创建一个向量，在命令窗口下输入：t=0:1:10屏幕显示：t=012345678910注意：向量的第一个元素的下标是1，而不是0。T=0:1:10产生了010步长为1的共11个数，保存在t(1),t(2),t(11)中。,13,定数线性采样法x=linspace(a,b,n),x=linspace(0,pi,11)x=Columns1through700.31420.62830.94251.25661.57081.8850Columns8through112.19912.51332.82743.1416,14,矩阵的创建【例】创建一个33的矩阵，输入：a=123;456;789或a=1,2,3;4,5,6;7,8,9或a=1,2,34,5,67,8,9屏幕输出a=123456789,15,整个输入数组必须以方括号“”为其首尾；数组的行与行之间必须用分号“;”或回车键【Enter】隔开；数组元素必须由逗号“,”或空格分隔。,a=2.7358;b=33/79;C=1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+iC=1.00005.4716+0.4177i0.69090.70714.82443.5000+1.0000i,16,标准数组生成函数,17,访问矩阵的某一个元素：a(3,2)ans=8访问矩阵某行(列）的所有元素：a(3,:)a(:,2)注意：这两个符号分别是1*3行向量和3*1列向量,18,【例】子数组的寻访。rand(state,0)%把均匀分布伪随机发生器置为0状态x=rand(1,5)%产生的均布随机数组x=0.95010.23110.60680.48600.8913x(3)%寻访数组x的第三个元素。ans=0.6068x(125)%寻访数组x的第一、二、五个元素组成的子数组。ans=0.95010.23110.8913x(1:3)%寻访前三个元素组成的子数组ans=0.95010.23110.6068,19,【例】创建二维数组C。a=2.7358;b=33/79;%这两条指令分别给变量a，b赋值。C=1,2*a+i*b,b*sqrt(a);sin(pi/4),a+5*b,3.5+i%这指令用于创建二维数组CC=1.00005.4716+0.4177i0.69090.70714.82443.5000+1.0000i,20,【例】二维数组的子数组寻访和赋值。A=zeros(2,4)%创建的全零数组A=00000000A(:)=1:8%全元素赋值方式A=13572468,21,3.3部分特殊变量和常数,ans最近生成的无名结果eps浮点数的相对误差pi3.14145926i虚数单位j虚数单位inf无穷大，如n/0NaN非数,真实记述0/0,/,0运算的后果；,22,a=0/0Warning:Dividebyzero.(TypewarningoffMATLAB:divideByZerotosuppressthiswarning.)a=NaNisnan(a)ans=1,23,t=-2*pi:pi/10:2*pi;%该自变量数组中，存在0值y=sin(t)./t;%在t=0处，计算将产生NaNtt=t+(t=0)*eps;%使0元素被一个“机器零”小数代替yy=sin(tt)./tt;%用sin(eps)/eps近似代替sin(0)/0极限subplot(1,2,1),plot(t,y),axis(-7,7,-0.5,1.2),xlabel(t),ylabel(y),title(残缺图形)subplot(1,2,2),plot(tt,yy),axis(-7,7,-0.5,1.2),24,3.4部分常用运算符,3.4.1算术运算符加减*乘(包括标量乘，矩阵乘，标量与矩阵乘)/除(包括标量除，矩阵除标量，数组除标量)矩阵求幂（矩阵必须为方阵）.*数组相乘./数组相除.数组求幂注意：.*./表示两个同维数组中的对应元素做乘，除.表示对数组的每个元素求幂,25,【例】输入A=1,3;2,4B=1,2;3,4C=A*B屏幕输出：A=1324B=1234C=10141420,26,C=A.*B屏幕输出：A=1324B=1234C=16616,27,【例】演示pow2的数组运算性质。A=1:4;5:8%生成数组A=12345678pow2(A)%计算的结果也是数组ans=248163264128256,28,【例】复数矩阵的生成及运算A=1,3;2,4-5,8;6,9*iB=1+5i,2+6i;3+8*i,4+9*iC=A*BA=1.0000-5.0000i3.0000-8.0000i2.0000-6.0000i4.0000-9.0000iB=1.0000+5.0000i2.0000+6.0000i3.0000+8.0000i4.0000+9.0000iC=1.0e+002*0.99001.1600-0.0900i1.1600+0.0900i1.3700,29,3.4.2关系运算符,大于=大于等于=等于=不等于运算法则：若关系式成立，结果为1；若关系式不成立，结果为0。,30,A=1:9,B=10-A,r0=(A4),r1=(A=B)A=123456789B=987654321r0=111000000r1=000010000,31,3.4.3逻辑运算符,n=0:(length(x)-1);w=(2*pi/length(x)*n;Am3,pha3=dft3(x);语句后不加；号可在command窗口看到结果。,35,第二种：不可单独运行，但能被其他M文件直接调用（须在MATLAB可寻的路径，如同一目录），相当于函数。%dft3.m:functionAm,pha=dft3(x)Xk=fft(x);Am=abs(Xk);pha=angle(Xk);,36,在MATLAB中，调用函数的常用形式是：输出参数1,输出参数2,=函数名(输入参数1,输入参数2,)函数调用可以嵌套，一个函数可以调用别的函数，甚至调用它自己（递归调用）。,37,程序流程控制,顺序结构：Matlab从上到下依次执行各语句，该结构最简单。循环结构(1)for-end循环用于循环次数事先确定的，格式为fori=n:s:m语句体ends为步长，可以为正数，负数或小数。,38,（2）while-end循环用于循环次数不能事先确定的，格式为while表达式语句体end只有表达式为真，就执行语句体，表达式为假，终止该循环。,39,分支结构（1）if语句格式为if表达式1语句体1elseif表达式2语句体2.else语句体elseend,40,（2）switch-case结构。,41,3.6m文件的创建，保存，运行,创建【File】菜单下【New】菜单选项的【M-File】命令打开matlab的m文件编辑器窗口。,42,保存单击M文件编辑器窗口工具栏中的【Save】图标，打开保存对话框,43,运行使M文件*.m所在目录成为当前目录，或让该目录处在MATLAB的搜索路径上，然后运行M文件名指令，便可。在command窗口可看到调试信息。command窗口可及时清理。,44,3.7基本数学函数和基本绘图函数,sin，cos，tan，asin(x)，acos(x)，atan(x)，atan2(x，y)，sinh(x)，cosh(x)，tanh(x)，asinh(x)，acosh(x)，atanh(x)：三角函数exp:指数sqrt：开平方log:自然对数log10:产用对数log2:以二为底的对数abs:取模Angle：取角度；conj:取共轭real:取实部；imag:取虚部round(x)，fix(x)，floor(x)，ceil(x)：取整,45,基本绘图指令,Matlab提供了强大的图形绘制功能。在大多数情况下，用户只需要指定绘图的方式，提供绘图数据，利用Matlab提供的丰富的二维，三维图形函数，就可以绘制出所需的图形。1、绘制二维连续函数Matlab中最常用的绘图函数是plot，plot的基本的命令格式有以下几种：,46,(1)plot(y),plot(y,S)当y为一向量时，以y的序号作为x轴，按向量y的值绘制曲线。S:用来指定线型、色彩、数据点形的选项字符串。(2)plot(x,y),plot(x,y,S)x,y均为向量时，以x向量作为X轴，向量y作为Y轴绘制曲线。注意：x和y种元素的个数必须相同！,47,2、绘制二维离散序列在Matlab用stem命令实现离散序列的绘制。stem命令的格式有以下几种:(1)stem(y);stem(y，S)以x=1,2,3为各点数据的x坐标，以y向量的各个对应元素为y坐标，在（X,Y）坐标面画一个空心小圆圈。(2)stem(x,y）；stem(x,y,S)以x向量的各个元素为x坐标，以y向量的各个对应元素为y坐标，在（X,Y）坐标面画一个实心小圆圈，并连接一条线段到X轴。,48,3、绘制子图：figure（n）;subplot(m,n,l);4、坐标控制：axis(0,1,-5,5);5、图上加标题:title(ActualPhaseResponseinradians);6、坐标标注:xlabel(frequencyinpiunits);ylabel(arg(Hexp(jw);7、在图面(xt,yt)坐标处书写字符注释text(xt,yt,S),49,8、gridon：画出分格线gridoff：不画分格线grid：是否画分格线的双向切换指令9、holdon：使当前轴及图形保持而不被刷新，准备接受此后将绘制的新图形holdoff：使当前轴及图形不再具备不给刷新的性质hold：当前图形是否具备刷新性质的双向切换开关,50,10、gtext：用鼠标把字符串放置到图形上11、legend：在指定位置建立图例,51,一般步骤,52,53,S：色彩和线型,54,S：数据，点形,55,3.8编程举例,例一：画出衰减震荡曲线及其包络线,t的取值范围是程序如下：,56,程序运行结果如下：,57,58,例二.用图形表示连续调制波形及其包络。程序如下：,59,60,例三.用图形表示连续调制波形及其包络。程序如下：,61,62,例四.x=linspace(0,2*pi,30);y=sin(x);z=cos(x);plot(x,y,b:p,x,z,r-,x,2*z,m+),63,例五.t=2*pi*(0:20)/20;y=cos(t).*exp(-0.4*t);stem(t,y,g);holdon;stairs(t,y,r);holdoff,64,65,四、Matlab在数字信号处理课程中的应用举例,常见连续信号离散信号离散卷积LTI系统,66,4.1常见连续信号,周期信号：正弦信号，周期方波【例】产生一个幅度为2，频率为4Hz,相位为的正弦信号,67,非周期信号：指数信号，阶跃信号，取样函数Sa(x)等【例】产生一个高度为1，宽度为3，延时为2秒的矩形脉冲信号,68,【例】取样函数：用sinc(x)命令,69,70,随机信号：【例】rand产生均匀分布的白嘈声，randn产生高斯分布的白噪声,71,4.2离散信号,周期序列【例】产生一个幅度为2，频率为4Hz,相位为的正弦信号,72,单位脉冲序列和单位阶跃序列u(n)可以用ones(1,n)和zeros(1,n)来生成单位脉冲序列和单位阶跃序列ones(1,n)命令产生1行n列的1值zeros(1,n)命令产生1行n列的0值,73,【例】产生单位脉冲序列和单位阶跃序列u(n-3),74,75,4.3离散卷积,matlab中卷积运算的指令是c=conv(a,b),其中a和b是有限长的序列。【例】计算下面两个序列的卷积a=-20113;b=120-1;两个序列的都从0开始的,76,77,78,【例】有序列和求这两个完整序列的卷积，并图示。,79,%完整序列卷积%conv1.ma=ones(1,10);n1=3;n2=12;%完整a(n)序列的非平凡值和区间端点b=ones(1,8);n3=2;n4=9;%完整b(n)序列的非平凡值和区间端点c=conv(a,b);nc1=n1+n3;nc2=n2+n4;%计算卷积和确定卷积非平凡区间端点kc=nc1:nc2;%构成非平凡区间的序号自变量stem(kc,c,b)title(convofxandh);ylabel(y);xlabel(Timeindexn);,80,81,4.4LTI系统,filter：因果LTI系统的的零状态响应，在matlab中可以用函数y=filter(num,den,x)实现。其中：num-系统传递函数分子系数组成的行向量den-系统传递函数分母系数组成的行向量x-输入的离散序列y-输出的离散序列，y的长度与x的长度一样,82,impz：在A(z)、B（z）已知情况下,用以下函数求系统的单位抽样响应h(n)，t是所记录的0到n-1的取样点数矢量：h,t=impz(num,den)，取样点数自动选取；h,t=impz(num,den,n)，取样点数选取为0到n-1或由矢量n指定，取样间隔为1；h,t=impz(num,den,n，Fs)，取样间隔为1/Fs；impz(num,den)，只在当前窗口用stem（t，h）绘出单位抽样响应h(n)。,83,freqz：因果LTI系统的的幅频响应和相频响应曲线，在matlab中可以用以下函数作出：H,w=freqz(num,den,N)，频率响应在单位圆的上半周取N个等分点；H,w=freqz(num,den,N,whole)，频率响应在单位圆取N个等分点；H=freqz(num,den,w)，频率响应在w设定的频率上。该函数可以同时作出幅频和相频响应图其中：num-系统传递函数分子系数组成的行向量；den-系统传递函数分母系数组成的行向量；N-是频率响应的点数，最好为2的幂缺省值512；w-是返回频率轴坐标向量。,84,zplane：离散系统的极零图，在matlab中可以用以下函数作出：zplane(z,p)，已知极点、零点；zplane(num,den)；hz，hp，ht=zplane(z,p)，返回极点、零点、坐标轴等的句柄。,85,【例】已知一个系统的传递函数为绘出极零图，求系统单位冲激响应hn，以及hn的幅频相频响应图。,86,87,figure(1)subplot(2,1,1);stem(n,h1,fill);gridon;title(单位冲激响应1);ylabel(h1(n);xlabel(n)subplot(2,1,2);h2,t=impz(num,den,n,Fs);%计算单位冲激响应2stem(t,h2,fill);gridon;title(单位冲激响应2);ylabel(h2(n);xlabel(n),88,figure(2)zplane(num,den);%作出极零图figure(3)freqz(num,den,N,Fs);gridon;%做出幅频和相频响应曲线title(频率响应);,89,90,91,92,【例】有一调幅信号用做频谱分析，要求能分辨xa(t)的所有频率分量，若用fs=3kHz频率抽样，抽样数据为点，做频谱分析，求X(k)=DFTx(n)，点，画出X(k)的幅频特性|X(k)|。,93,%cosfft0.mN=512;n=0:1:N-1;fs=3000%x=cos(2*pi/fs*100*n);x=cos(2*pi/fs*600*n)+0.5*cos(2*pi/fs*700*n)+0.5*cos(2*pi/fs*500*n);subplot(3,1,1);stem(n,x);ylabel(signalx(n),0=n=511);xlabel(n)axis(0,N-1,-2.5,2.5),94,X=fft(x,N);magX=abs(X(1:1:N)k=0:1:N-1;f=fs/N*k;axis(0,1500,0,250)subplot(3,1,2);plot(f,magX);ylabel(DFTMagnitude);xlabel(frequencyinHz)%图中横坐标值即为信号的实际频率（赫兹）。k=0:1:N-1;w=2*pi/N*k;subplot(3,1,3);plot(w,magX);ylabel(DFTMagnitude);xlabel(frequencyinrad)axis(0,pi,0,250)%图中横坐标值为信号的数字频率（弧度），乘以fs除以2*pi即为信号的实际频率（赫兹）。,95,96,【例】比较高密度谱和高分辨率谱之间的差异，考虑x（n）=cos(0.48*pi*n)+cos(0.52*pi*n)（）取x（n）（）时，求x（n）的DFT。（）将（）中的x（n）以补零方式使x（n）加长到，求x（n）的DFT。（）取x（n）（）时，求x（n）的DFT画出幅频特性。,97,%spectrum_com1.mxlabel(frequencyinpiunits)axis(0,1,0,10),98,%part2%Highdensityspectrum(100samples)basedonthefirst10samplesofx(n)figure(2)n3=0:1:99;x=cos(0.48*pi*n3)+cos(0.52*pi*n3);y3=x(1:1:10)zeros(1,90);subplot(2,1,1);stem(n3,y3);title(signalx(n),0=n=9+90zeros);xlabel(n)axis(0,100,-2.5,2.5)Y3=fft(y3);magY3=abs(Y3(1:1:51);k3=0:1:50;w3=2*pi/100*k3;subplot(2,1,2);plot(w3/pi,magY3);title(DTFTMagnitude);xlabel(frequencyinpiunits)axis(0,1,0,10),99,%part3%Highresolutionspectrumbasedon100samplesofthesignalx(n)figure(3)n=0:1:99;x=cos(0.48*pi*n)+cos(0.52*pi*n);subplot(2,1,1);stem(n,x);title(signalx(n),0=n=99);xlabel(n)axis(0,100,-2.5,2.5)X=fft(x);magX=abs(X(1:1:51);k=0:1:50;w=2*pi/100*k;subplot(2,1,2);plot(w/pi,magX);title(DTFTMagnitude);xlabel(frequencyinpiunits)axis(0,1,0,60),100,101,102,103,【例】有一连续信号频率成份分别为f1=4Hz，f2=8Hz，f3=10Hz。（1）、若用fs=32Hz频率抽样，抽样数据为N=16点做频谱分析，能分辨xa(t)的所有频率分量。X1（t）一个周期取8点，共取2个周期；X2（t）一个周期取4点，共取4个周期；X3（t）一个周期取3.2点，共取5个周期。或N取16的整数倍，但最少为16。,104,（2）、若用fs=64Hz频率抽样，抽样数据为N=32点做频谱分析，也能分辨xa(t)的所有频率分量。X1（t）一个周期取16点，共取2个周期；X2（t）一个周期取8点，共取4个周期；X3（t）一个周期取6.4点，共取5个周期。或N取32的整数倍，但最少为32。抽样频率fs可取32Hz的整数倍，但最少为32Hz。,105,若用fs=30Hz频率抽样，抽样数据为16点做频谱分析，不能分辨xa(t)的所有频率分量X(k)的幅频特性|X(k)|如下。,106,%cosdft30.m:N=16;n=0:1:N-1;k=0:1:N-1;fs=32;%x=1,1,1,1;x=cos(2*pi/fs*4*n)+2*cos(2*pi/fs*8*n)+3*