第9次课-信号处理算法Matlab仿真

第9讲MATLAB知识回顾数字信号处理与Matlb仿真题型一、填空：10个题，每空2分，共20分；二、选择：5个题，每题2分，共10分；三、分析与计算题（共4题,每题6分，共24分）三、综合编程题（共4题,共46分）引言MATLAB软件介绍：MATLAB是矩阵实验室（MatrixLaboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。版本更新4Matlab的特点与功能

Matlab具有很强的数值计算功能

Matlab以矩阵作为数据操作的基本单位，但无需预先指定

矩阵维数（动态定维）按照IEEE的数值计算标准进行计算提供十分丰富的数值计算函数，方便计算，提高效率

Matlab命令与数学中的符号、公式非常接近，可读性强，

容易掌握

Matlab是一个交互式软件系统输入一条命令，立即就可以得出该命令的结果5

Matlab符号计算功能Matlab和著名的符号计算语言Maple相结合

Matlab的编程功能Matlab具有程序结构控制、函数调用、数据结构、输入输出、面向对象等程序语言特征，而且简单易学、编程效率高。

Matlab的绘图功能Matlab提供丰富的绘图命令，

很方便实现数据的可视化Matlab的特点与功能6

Matlab丰富的工具箱（toolbox）根据专门领域中的特殊需要而设计的各种可选工具箱

Matlab的Simulink动态仿真集成环境提供建立系统模型、选择仿真参数和数值算法、启动仿真程序对该系统进行仿真、设置不同的输出方式来观察仿真结果等功能。SymbolicMathPDEOptimizationSignalprocessImageProcessStatisticsControlSystemSystemIdentification…

…Matlab的特点与功能7主要应用数值分析数值和符号计算工程与科学绘图控制系统的设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程89界面介绍1011一、命令窗的常用控制指令指令含义指令含义cd设置当前工作目录。exit关闭/退出MATLABclf清除图形窗quit关闭/退出MATLABclc清除指令窗中显示内容more使其后的显示内容分页进行clear清除MATLAB工作空间中保存的变量return返回到上层调用程序；结束键盘模式dir列出指定目录下的文件和子目录清单type显示指定M文件的内容edit打开M文件编辑器which指出其后文件所在的目录11二、MATLAB基本操作变量与矩阵矩阵及元素的赋值矩阵运算1、变量---命名规则变量名是对字母大小写敏感的。如:变量myvar和MyVar是不同变量。sin是MATLAB定义的正弦名，但SIN、Sin等都不是。变量名的第一个字母必须是英文字母，最多可包含31个字符（英文、数字和下划线）。如：myvar201合法。变量名中不得含空格、标点，但可以有下连符号如:变量名my_var_201合法的，且读起来更方便。

而my,var201不合法。132.矩阵及元素的赋值Matlab中的变量或常量都代表矩阵，标量看作1×1阶的矩阵赋值语句：变量=表达式（或数）

a=[123;456;789]元素的输入可以用表达式代替

x=[-1.3,sqrt(3),(1+2+3)/5*4]矩阵的值放在方括号中，同一行中各元素之间用逗号或空格分开，不同的行则以分号隔开，此时会立即显示运算结果。如果不希望显示结果，以分号结尾，此时结果在workspace可以看到。14变量的元素用()中的数字标明，一维矩阵中的元素用一个下标表示，二维矩阵可有两个下标数，以逗号隔开。三维或更高维的矩阵可有三个或更多下标。可单独给元素赋值

x(2)=1.7321a(2,3)=6linspace(a,b,n)%在a和b之间均匀153.矩阵运算（一）预定义变量（二）运算符+

-

*

/

\

^

sqrt

~=“点乘”与普通乘法的区别；‘.*’and‘*’2.运算优先级：指数运算>乘除>加减；3.“=”和其他运算符两侧，允许有空格，以增加程序可读性。（三）矩阵（数组）生成函数zeros(m,n)生成一个m

行n

列的零矩阵，m=n

时可简写为zeros(n)ones(m,n)生成一个m行n列的元素全为1的矩阵,

m=n

时可写为ones(n)eye(m,n)生成一个主对角线全为1的m

行n

列矩阵,

m=n

时可简写为eye(n)，即为n

维单位矩阵diag(X)若X是矩阵，则diag(X)为X的主对角线向量若X是向量，diag(X)产生以X为主对角线的对角矩阵tril(A)提取一个矩阵的下三角部分triu(A)提取一个矩阵的上三角部分rand(m,n)产生0～1间均匀分布的随机矩阵，m=n

时简写为rand(n)randn(m,n)产生均值为0，方差为1的标准正态分布随机矩阵，

m=n

时简写为randn(n)随机矩阵为了进行信号分析与处理，经常需要对接收信号进行仿真，而在信号仿真中离不开随机数的产生，MATLAB提供的rand和randn函数可分别产生均匀分布和正态分布的随机数。例如，要产生[0，1]之间均匀分布的随机向量r1(2×3)，可输入：>>r1=rand(2,3) ％产生[0，1]之间均匀分布的随机矩阵

r1=0.45650.82140.61540.01850.44470.791919>>r2=5-10*rand(2,3)

％产生[-5，5]之间均匀分布的随机矩阵

r2=0.89734.4211-3.1317

-3.93651.47134.9014>>r3=randn(2,3)％产生均值为0、方差为1的标准正态分布的随机矩阵r3=1.18920.3273-0.1867

-0.03760.17460.725820通过下标引用矩阵的元素，例如A(3,2)=200矩阵元素按列存储，例如

A=[1,2,3;4,5,6];

A(3)

ans=

2显然，序号(Index)与下标(Subscript)是一一对应的，以m×n矩阵A为例，矩阵元素A(i,j)的序号为(j-1)*m+i。其相互转换关系也可利用sub2ind和ind2sub函数求得。矩阵元素（四）矩阵操作21A(r,c) A第r行第c列元素A(r,:) A第r行所有列元素A(:,c) A所有行第c列元素A(:,k:m)A的第k到第m列组成的子矩阵A(:) 矩阵所有元素排列为一维向量（从左到右 按列连接）A(s) 按序号寻访元素A(r,c)=K r行c列元素赋值为KA(:)=D(:) A、D两矩阵总元素数相同（行列不一 定），将D每个元素赋给A中每个元素A(:)与A(:,:)的区别?如何获得由A的第一、三行和第一、二列组成的子矩阵？自己动手2323查看矩阵的大小：size与length

size(A)

列出矩阵A的行数和列数

size(A,1)

返回矩阵A的行数

size(A,2)

返回矩阵A的列数例：>>

A=[123;456]>>

size(A)>>

size(A,1)>>

size(A,2)

length(x)返回向量

X的长度

length(A)等价于max(size(A))24改变矩阵的形状：reshapereshape(A,m,n):将矩阵元素按列方向进行重组重组后得到的新矩阵的元素个数必须与原矩阵元素个数相等！25矩阵基本运算

矩阵的加减：对应分量进行运算要求参与加减运算的矩阵具有相同的维数例：>>

A=[123;456];B=[321;654]>>

C=A+B;D=A-B;

矩阵的普通乘法要求参与运算的矩阵满足线性代数中矩阵相乘的原则例：>>

A=[123;456];B=[21;34;56]>>

C=A*B26矩阵的数组运算

数组运算：对应元素进行运算点与算术运算符之间不能有空格！

数组运算包括：点乘、点除、点幂

相应的数组运算符为：“.*”，“./”，“.\”和“

.^”参与运算的对象必须具有相同的形状！例：>>

A=[123;456];B=[321;654];>>

C=A.*B;D=A./B;E=A.\B;F=A.^B;27线性代数运算的MATLAB命令MATLAB是矩阵化程序设计语言，所以处理矩阵和向量运算特别方便。关于矩阵和向量的一些基本运算命令已在前面有所介绍，常用的命令和函数还有zeros生成0矩阵eig特征值、特征向量ones生成1矩阵diag对角矩阵eye生成单位矩阵trace方阵的迹linspace生成等距行向量rank矩阵的秩rand生成随机矩阵rref行最简形det方阵的行列式orth正交规范inv方阵的逆null求基础解系norm范数jordanJordan分解cond方阵的条件数28Matlab中常见数学函数sin、cos、tan、cot、sec、csc、…asin、acos、atan、acot、asec、acsc、…exp、log、log2、log10、sqrtabs、conj、real、imag、signfix、floor、ceil、round、mod、remmax、min、sum、mean、sort、fftnorm、rank、det、inv、eig、lu、qr、svd……log

是自然对数，即以

e

为底数mod(x,y)

结果与y

同号，rem(x,y)

则与

x

同号max

等函数的参数是矩阵时，是作用在矩阵各列上29二、Matlab绘图作为一个功能强大的工具软件，Matlab具有很强的图形处理功能，提供了大量的二维、三维图形函数。由于系统采用面向对象的技术和丰富的矩阵运算，所以在图形处理方面既方便又高效。1概述及图形窗口的控制【1】figure命令格式：figure或figure(N)功能：创建一个图形窗口。使编号为N的图形窗口成为当前图形窗口，即图形窗口处于可视状态。如果窗口N不存在，则将创建一个句柄为N的图形窗口。31311概述及图形窗口的控制【2】clf命令格式：clf功能：清除当前图形窗口中的所有的内容例>>x=-pi:pi/10:pi;>>y=sin(x);>>plot(x,y)>>clf32321概述及图形窗口的控制【3】close命令格式一：close功能：

关闭当前图形窗口格式二：close(N)功能：

关闭指定编号N的图形窗口格式三：closeall功能：

关闭所有图形窗口33二维图形plot函数函数格式：plot(x,y)

其中x和y为坐标向量函数功能：以向量x、y为轴，绘制曲线。例：在区间0≤x≤2内，绘制正弦曲线y=sin（x），其程序为：x=0:pi/100:2*pi;y=sin(x);plot(x,y)34线型与颜色格式：plot(x,y1,’cs’,...)其中c表示颜色，s表示线型。例：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,'g-',x,y2,'b-.')其中参数‘g-’和‘b-.’表示图形的颜色和线型。g表示绿色，-表示图形线型为实线；b表示蓝色，-.表示图形线型为点划线。35设置曲线的样式属性色彩（c）说明线型（l）说明点型/标记（m）说明y黄色-实线.点m品红--长划线o圆圈c青色：虚线x十字r红色-.点划线+加号g绿色*星号b蓝色^上三角w白色v下三角k黑色>右三角<左三角s矩形d菱形36>>figure,>>t=0:pi/20:2*pi;>>y=sin(t);>>y1=sin(t-pi/2);>>y2=sin(t-pi);>>plot(t,y,'c-v',t,y1,'k--s',t,y2,'m:o')37若使用plot函数时不指定线型，则绘制的曲线将不被连接起来>>figure,>>t=0:pi/20:2*pi;>>y=sin(t);>>y1=sin(t-pi/2);>>y2=sin(t-pi);>>plot(t,y,'cv',t,y1,'ks',t,y2,'mo')3838使用命令grid，可以将图形窗体的轴在有网格线和无网格线之间切换gridon——图形窗体的轴显示坐标网格线gridoff——图形窗体的轴不显示坐标网格线3939>>t=0:pi/20:2*pi;>>y=sin(t);>>y1=sin(t-pi/2);>>y2=sin(t-pi);>>plot(t,y,'cv',t,y1,'ks',t,y2,'mo')>>gridon>>gridoff40图形标记在绘制图形的同时，可以对图形加上一些说明，如图形名称、图形某一部分的含义、坐标说明等，将这些操作称为添加图形标记。title('加图形标题');xlabel('加X轴标记');ylabel('加Y轴标记');

41设定坐标轴用户若对坐标系统不满意，可利用axis命令对其重新设定。axis([xminxmaxyminymax])设定最大和最小值axis(‘auto’)将坐标系统返回到自动缺省状态axis(‘square’)将当前图形设置为方形axis(‘equal’)两个坐标因子设成相等axis(‘off’)关闭坐标系统axis(‘on’)显示坐标系统42例：在坐标范围0≤X≤2π,-2≤Y≤2内重新绘制正弦曲线，其程序为：x=linspace(0,2*pi,60);%生成含有60个数据元素的向量Xy=sin(x);plot(x,y);axis([02*pi-2

2]);%设定坐标轴范围43加图例给图形加图例命令为legend。该命令把图例放置在图形空白处，用户还可以通过鼠标移动图例，将其放到希望的位置。格式:legend('图例说明','图例说明');

例：为正弦、余弦曲线增加图例，其程序为：x=0:pi/100:2*pi;y1=sin(x);y2=cos(x);plot(x,y1,x,y2,'--');legend('sin(x)','cos(x)');44subplot（m,n,p）该命令将当前图形窗口分成m×n个绘图区，即每行n个，共m行，区号按行优先编号，且选定第p个区为当前活动区。subplot函数45x=linspace(0,2*pi,60);y=sin(x);z=cos(x);t=sin(x)./(cos(x)+eps);%eps为系统内部常数ct=cos(x)./(sin(x)+eps);subplot(2,2,1);plot(x,y);title('sin(x)');axis([02*pi-11]);subplot(2,2,2);plot(x,z);title('cos(x)');axis([02*pi-11]);subplot(2,2,3);

plot(x,t);title('tangent(x)');axis([02*pi-4040]);subplot(2,2,4);plot(x,ct);title('cotangent(x)');axis([02*pi-4040]);例：在一个图形窗口中同时绘制正弦、余弦、正切、余切曲线46多图形窗口需要建立多个图形窗口，绘制并保持每一个窗口的图形，可以使用figure命令。每执行一次figure命令，就创建一个新的图形窗口，该窗口自动为活动窗口，若需要还可以返回该窗口的识别号码，称该号码为句柄。句柄显示在图形窗口的标题栏中，即图形窗口标题。用户可通过句柄激活或关闭某图形窗口，而axis、xlabel、title等许多命令也只对活动窗口有效。47重新绘制上例4个图形，程序变动后如下：x=linspace(0,2*pi,60);y=sin(x);z=cos(x);t=sin(x)./(cos(x)+eps);ct=cos(x)./(sin(x)+eps);figure;plot(x,y);title('sin(x)');axis([02*pi-11]);figure;plot(x,z);title('cos(x)');axis([02*pi-11]);figure;

plot(x,t);title('tangent(x)');axis([02*pi-4040]);figure;

plot(x,ct);title('cotangent(x)');axis([02*pi-4040]);48hold命令若在已存在图形窗口中用plot命令继续添加新的图形内容，可使用图形保持命令hold。发出命令holdon后，再执行plot命令，在保持原有图形或曲线的基础上，添加新绘制的图形。49阅读如下程序：x=linspace(0,2*pi,60);y=sin(x);z=cos(x);plot(x,y,'b');%绘制正弦曲线holdon;%设置图形保持状态plot(x,z,'r');

%保持正弦曲线同时绘制余弦曲线axis([02*pi-11]);legend('sin','cos');holdoff%关闭图形保持5051二维绘图函数小结plot二维图形基本函数fplotf(x)函数曲线绘制fill填充二维多边图形polar极坐标图bar条形图loglog双对数坐标图semilogxX轴为对数的坐标图semilogyY轴为对数的坐标图stairs阶梯形图axis设置坐标轴clf清除图形窗口内容close关闭图形窗口figure创建图形窗口grid放置坐标网格线gtext用鼠标放置文本hold保持当前图形窗口内容subplot创建子图text放置文本title放置图形标题xlabel放置X轴坐标标记ylabel放置Y轴坐标标记52三维图形plot3函数最基本的三维图形函数为plot3，它是将二维函数plot的有关功能扩展到三维空间，用来绘制三维图形。函数格式：plot3(x1,y1,z1,'c1',x2,y2,z2,'c2',…)其中x1,y1,z1…表示三维坐标向量，c1，c2…表示线形或颜色。函数功能：以向量x，y，z为坐标，绘制三维曲线。53绘制三维网线图mesh(z)：z为m×n的矩阵Z坐标——矩阵中元素的值X与Y坐标——矩阵中元素的下标X——1：n，对应矩阵的列Y——1：m，对应矩阵的行54图像的读/写用imread函数读取图像，调用格式为：

x=imread(‘filename.fmt’)filename是图像的文件名，fmt指定图像文件的格式。图像读、写与显示图像的读/写图像读、写与显示图像的读/写55a=imread('jlu.jpg');size(a)运行结果ans=2082083如：56在MATLAB中，imshow函数用于图像显示调用格式为：imshow(x)例子：x=imread(lena.jpg');figure,imshow(x);图像读、写与显示图像的显示57三、Matlab程序设计三、Matlab程序设计常用指令流程控制语句四种流程语句：if语句switch语句while语句for语句选择结构循环结构60（一）选择结构选择结构的语句有if语句和switch语句。1．if语句格式一：

if条件语句组

end%必不可少，若没有则表达式为0时找不到

%继续执行的程序入口格式二：

if条件语句组1else

语句组2

end61格式三：

if条件1

语句组1elseif条件2

语句组2……elseif条件m

语句组melse

语句组m+1end62例：输入三角形的三条边，求面积。

A=input('请输入三角形的三条边：');ifA(1)+A(2)>A(3)&A(1)+A(3)>A(2)&A(2)+A(3)>A(1)p=(A(1)+A(2)+A(3))/2;s=sqrt(p*(p-A(1))*(p-A(2))*(p-A(3)));disp(s);elsedisp('不能构成一个三角形。')end运行：请输入三角形的三条边：[456]9.921663（二）循环结构实现循环结构的语句：for语句和while语句。1.for语句：格式：

for循环变量=表达式1:表达式2:表达式3

循环体语句

end表达式1的值为循环变量的初值;表达式2的值为步长;表达式3的值为循环变量的终值。步长为1时，表达式2可以省略。64例：计算1+3+5+…+99=?sum=0;forn=1:2:99sum=sum+n;endsum运行结果：sum=250065例：当n=100时，求的值。程序如下:ticy=0;n=100;fori=1:ny=y+1/(2*i-1);endtocyy=3.284366在实际MATLAB编程中，采用循环语句会降低其执行速度，所以前面的程序通常由下面的程序来代替：ticn=100000;k=1:2:(2*n-1);y=sum(1./k);tocyy=3.284367672.while语句格式为：

while(条件)

循环体语句

end例：计算1~10的和i=1;sum=0;while(i<=10)sum=sum+i;i=i+1;endsum运行结果：5568M文件及程序调试在入门阶段，通常在行命令模式下工作，键入一条命令后，让系统立即执行该命令。该方法程序可读性差且难以存储。对于复杂的问题应编成可存储的程序文件，让Matlab执行该程序文件，这种工作模式成为程序文件模式。由Matlab语句构成的程序文件称为m文件，用m作为文件的扩展名，可直接阅读并用任何编辑器来建立。m文件分主程序和子程序主程序：由用户为解决问题而编写的子程序：函数文件，必须由其他m文件来调用69主程序：用clear、closeall等语句开始，清除工作空间中原有的变量和图形，以避免其他已执行的程序残留数据对本程序的影响。主程序中如有全局变量，即在子程序中与主程序共用的变量，在程序起始部分注明

global变量名1变量名2……整个程序按Matlab标识符的要求起名，加上后缀.m70函数文件用来定义子程序由function起头，后跟的函数名必须与文件名相同有输入输出变量，可进行变量传递除非用global声明，子程序中的变量均为局部变量，不保存在工作空间中71functiony=mean_my(x)%定义函数名，输入和输出变量[m,n]=size(x);ifm==1m=n;endy=sum(x)/mm，n，y为局部变量在主程序中调用：x=…a=mean_my(x)…72变量的存储和下载存储：save：把工作空间的变量存入磁盘，.matsaveaa%把工作空间的全部变量存入保存在名为aa.mat文件中;saveaaabc%把a，b，c三个变量保存在名为aa.mat文件中;下载：load，把mat数据文件传回到工作空间中，但不能选择变量loadaa%把aa.mat文件中的全部数据和变量名一起下载到工作空间中73四、离散傅立叶变换(DFT)74定义：有限长序列x(n)(0≤n≤N-1)有N个样本值。它的傅立叶变换X(jω)在频率区间(0≤ω<2π)的N个等间隔分布的点ωk=2πk/N(0≤k≤N-1)上也有N个取样值。这两组同样长的序列之间可以用简单的关系联系起来：其中：1、DFT定义75它的反变换称为离散傅立叶反变换(IDFT)，可以表为：1、DFT定义76若用矩阵式X=W*x及x=W-1*X求变换矩阵：W=dftmtx(N)，

反变换矩阵：W-1=conj(dftmtx(N))/N，但通常可直接调用FFT程序：正变换：X=fft(x,N) X=

fft(x)N默认为x的长度反变换：X=ifft(x)或X=

ifft(x,N)

4、用工具箱函数计算DFT77FFT程序求频谱的范围规范化为[0,2π)，对应的频谱下标为k=0,1,…N-1。设频率分辨率 dw=2π/N，则频率向量为

w=k*dw若要画出乃奎斯特频率范围[-π,π)内的频谱，需要把[π,2π)范围的频谱平移到[-π,0)内，对于这一运算，有专门的工具箱函数fftshift来实现。

X1=fftshift(X) X1就是在[-π,π)内的频谱向量。4、用工具箱函数计算DFT78傅里叶变换函数常用的FFT及反变换函数如下表所示。函数说明fft计算快速离散傅立叶变换fftshift调整fft函数的输出顺序，将零频位置移到频谱的中心abs求幅值angle求相角ifft计算离散傅立叶反变换7980FFT计算频谱举例例：选择合适的变换区间长度N，用DFT对下列信号进行谱分析，画出幅频、相频特性。解：(1)为周期序列，取一个周期用DFT谱分析%(1)周期序列,截取一个周期作谱分析N1=100;n1=0:N1-1;x1n=2*cos(0.2*pi*n1);%生成序列x1nX1k=fft(x1n,N1);%x1n的N1点DFTsubplot(1,2,1);stem(n1/N1*pi,abs(X1k),'.');title('X1k-Amplitude');su