MATLAB第4章数值运算基础.ppt

物理与电气工程学院,1,第4章 数值运算基础,MATLAB的计算包括数值计算和符号计算两类，本章将带大家学习数值计算部分。其中将主要学习与我们专业密切相关的多项式、方程组求解、数据分析和数字信号处理的快速傅里叶变换。,物理与电气工程学院,2,第1节 多项式polynomial MATLAB用行向量表示多项式。将多项式的系数按降幂次序存放在行向量中。如：,的系数行向量为,注意：多项式中缺少的幂次系数一定要用“0”补齐。,物理与电气工程学院,3,一、创建多项式 1、系数矢量直接输入法 在命令窗口直接输入多式的系数向量,【例4-1】输入系数矢量，创建多项式 x3-4*x2+3*x+2。 p=1 -4 3 2 poly2sym(p) %将矢量P表示为多项式的手写形式,Polynomial coefficient vector to symbolic polynominal,4,2、方阵特征多项式输入法 p=poly(A) 若A为nn的矩阵，则返回值P将是一个含有n+1个元素的行向量，也就是该矩阵特征多项式的系数,【例4-2】求矩阵1 2 3;4 5 6;7 8 0的特征多项式系数，并转换为多项式形式。 a=1 2 3;4 5 6;7 8 0; p=poly(a); poly2sym(p) %将矢量P表示为多项式的手写形式 d1=roots(p) %由特征多项式求得的特征值 d2=eig(a) %由特征值函数求得的特征值,5,3、根矢量创建 p=poly(A) A为待求多项式的根矢量，则返回值将是对应多项式的系数行矢量，该多项式的根为矢量A。此时p=poly(A)与A=roots(p)互逆。系统定义P0=1。,Ax1 x2 x3, p=1 p1 p2 p3,6,【例4-3】根据根矢量-0.5 -0.3+0.4i -0.3-0.4i，创建多项式 r=-0.5 -0.3+0.4i -0.3-0.4i; p=poly(r) pr=real(p) ppr=poly2sym(pr),二、多项式运算 1、求多项式的值 MATLAB的多项式求值方式有两种，按数组运算规则和按矩阵运算规则计算多项式的值。,y=polyval(p,x) 按数组规则运算，计算多项式在x处的值，p是多项式的系数矢量；x是指定自变量的值，可以是标量、向量或矩阵。如果x是向量或矩阵，则函数的返回值是与x对应的同维向量或矩阵。,物理与电气工程学院,7,【例4-4】求多项式3x2+2x+1在5、7和9处的值。 p = 3 2 1; polyval(p,5 7 9),y=polyvalm(p,x) 将矩阵整体（而不是矩阵元素）作为自变量进行计算。p是多项式的系数向量；相当于用矩阵x代替多项式的变量对矩阵而不是对数组进行计算，x必须是方阵,【例4-5】求多项式3x2+2x+1对于矩阵2 5;7 9的值 p= 3 2 1; polyvalm(p,2 5;7 9),物理与电气工程学院,8,2、求多项式的根 格式：C=roots(p) p为多项式的系数矢量，C为函数返回多项式的根矢量 如果C为复数，则必成对出现。,【例4-6】分别用两种方法求多项式 x5-5x4+3x3-6x2+4x-10的根。 a=1 -5 3 -6 4 -10; r=roots(a),物理与电气工程学院,9,3、多项式的乘除运算 多项式的乘法conv 格式： c=conv(a,b) 多项式的乘法运算,也是矢量的卷积运算 向量a长度为m，向量b长度为n，a和b的卷积定义为：,运算结果矢量c为长度k=m+n-1,【例4-7】计算两多项式x4-5x3+3x2-4x+2和x3+2x2-5x+3的乘法 a=1 -5 3 -4 2;b=1 2 -5 3; c=conv(a,b) poly2sym(c),10,多项式的除法deconv 格式： q, r=deconv( c,a) 多项式的除法运算，也是矢量的解卷积运算过程。 向量a对向量c进行解卷积得到“商”向量q和“余”向量r。 q和r分别是商多项式和余多项式；c和a分别是被除多项式和除多项式，使得：c=conv(a,q)+r,【例4-8】计算例4-7中求得的乘积被x3+2x2-5x+3除所得结果 c=1 -3 -12 30 -36 33 -22 6; b=1 2 -5 3; q ,r=deconv(c,b),物理与电气工程学院,11,4、多项式微积分 polyder(p) 返回多项式系数向量p 的导数,【例4-9】计算多项式3x4-5x3+2x2-6x+10的微分。 p=3 -5 2 -6 10; dp=polyder(p) poly2sym(dp),polyint(p) 返回多项式系数p 的积分,【例4-10】计算多项式12x3-15x2+4x-6的积分。 p=12 -15 4 -6; ip=polyint(p) poly2sym(ip),物理与电气工程学院,12,5、多项式的部分分式展开 MATLAB提供了residue命令来执行部分分式展开或多项式系数之间的转换。该命令通常用于信号与控制领域中。格式如下： r,p,k=residue(b,a) 该命令是求多项式之比b(s)/a(s)的部分分式展开，返回留数r、极点p和直项向量k。a和b分别是分母和分子多项式的系数向量,b,a=residue(r,p,k) 上一条语句的逆过程，只是分母多项式系数以归一化形式：最高次项系数a(1)为1。,13,如果分母多项式a(s)不含重根，则两个多项式可写成如下形式：,其中，pi称为极点，ri称为留数，k(s)是直项。 如果b的次数低于a的次数，则直项为空。,如果分母多项式a(s)含m重根p(j)=p(j+m-1),则这m项应展开为,物理与电气工程学院,14,b(x) 5x3+3x2-2x+7 【例3-11】两多项式的比为 = , a(x) -4x3+8x+3 求部分分式展开。 a = -4 0 8 3; b = 5 3 -2 7; r, p, k = residue(b,a) b1,a1 = residue(r,p,k) %分母最高次项归1 r2,p2,k2 = residue(1 1,1 -2 1) %出现重根,笔算结果=？,物理与电气工程学院,15,6、多项式拟合 对于给定的一组数据(xi ,yi)，i=1,2,n,如果要采用多项式模型对数据组进行描述，形成如多项式y(x)=f(x,p)=p1 xn+ p2 xn-1+ p3 xn-2+ pn+1的形式，求取参数p使得量值2(p)的值最小的过程，称为对数据组进行多项式拟合，其中,MATLAB系统设计polyfit函数采用最小二乘法原理对给定的数组(xi,yi)，i=1,2,n进行多项式的曲线拟合，最后给出拟合的多项式系数。,物理与电气工程学院,16,p=polyfit(x,y,n) 其中，x,y分别表示横、纵坐标向量；n是给定的拟合多项式的最高阶数，返回一个多项式系数向量p。 如n=3,若p=1 0.5 1 2,则 y=1*x3+0.5*x2+1*x1+2*x0,p,S=polyfit(x,y,n) 返回n阶多项式系数p, S为误差估计结构,p,S,mu=polyfit(x,y,n) 对归一化以后的数据进行多项式拟合，用x=(x-u1)/u2替代x，其中mu=u1 u2,u1=mean(x),u2=std(x),物理与电气工程学院,17,【例3-13】 求误差函数的6阶拟合多项式。 x = (0: 0.1: 2.5); % 生成0至2.5间隔为0.1的自变量 y = erf(x); % 计算误差函数 p = polyfit(x,y,6) % 求6阶拟合多项式 x = (0: 0.1: 5); % 生成0至5间隔为0.1的自变量 y1 = erf(x); % 计算误差函数 f = polyval(p,x); % 计算拟合函数的值 plot(x,y1,o,x,f,-) % 绘图函数 p0,s0,mu0 = polyfit(x,y,6) %x=(x-mean(x)/std(x) p1,s1,mu1 = polyfit(x-mu(1)/mu(2),y,6),可以看出，拟合区间0 2.5内拟合曲线拟近原曲线，,而在区间以外的曲线误差较大,第2节 线性代数 给定两个矩阵A和B,求X的解,使得: AX=B XA=B,在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“”和“/”求解： X=AB 是AX=B 的解 X=B/A 是XA=B 的解 对于方程AX=B，要求矩阵A和B有相同的行数；X和B有相同的列数，且它们的行数等于A的列数。,物理与电气工程学院,19,根据矩阵A的结构(m, n），可以将方程分为以下3类： m=n 方阵系统，可偿试求精确解 mn 超定系统，可偿试求最小二乘解 mn 欠定系统，可偿试求少m个非零解,“反斜线运算因子”对于不同结构矩阵A，将采用不同的算法处理。MATLAB会自动检测参数特性,选择运用不同的算法解方程组。,20,一、方阵系统 由n个未知数构成n个方程，方程有唯一的一组解，其一般形式写成如下形式：Ax=b 其中A是方阵，b是列向量 格式：x=Ab,A为非奇异，且条件数不大：x精确,A为非奇异，且条件数很大：注意解的可靠性,A为奇异：解不存在或存在但不唯一,A为接近奇阵或病态矩阵 Matrix is close to singular or badly scaled. Results may be inaccurate.,尽量使用格式：x=Ab,而不用格式：x=inv(A)*b,物理与电气工程学院,21,【例4-14】 a为方阵，b为矢量，求方阵系统的根。 a=fix(15*rand(3,3) b=fix(15*rand(3,1) det(a) cond(a) x=ab,【例4-15】 a和b均为方阵，求方阵系统的根。 a=fix(15*rand(3,3) b=fix(15*rand(3,3) x=ab,物理与电气工程学院,22,二、超定系统 方程组Ax=b，A为mn矩阵,如果A列满秩,且mn,则方程没有精确解,此时方程组为超定方程组。一般采用最小二乘法。,左除法：x=Ab 建立在奇异值分解基础之上,广义逆法：x=pinv(A)*b 速度较快，可靠性较差一些,实验中数据的曲线拟合就是就可以解超定方程组的方法来解。一般情况下需要将非线性问题转换为线性问题来解决,物理与电气工程学院,23,一组实验数据，时间t和测量数据y，如下表所示：,认为x和y有上式的关系式，则由6组实验数据就可形成超定方程组，就可求解出C1和C2。使得将t代入函数得到的值与实际y值之间差值的最小平方和最小,如何以矩阵的形式表示？,物理与电气工程学院,24,【例4-16】求表中数据的最小二乘解。 t=0 0.3 0.8 1.1 1.6 2.3; y=0.82 0.72 0.63 0.60 0.55 0.50; e=ones(size(t) exp(-t) %6组实验数据变换得到的 c=ey cp=pinv(e)*y t1=0:0.1:2.5; y1=ones(size(t1),exp(-t1)*c; yp=ones(size(t1),exp(-t1)*cp; plot(t,y,t1,y1,ro,t1,yp,b*),物理与电气工程学院,25,三 、欠定系统 程组Ax=b，A为mn矩阵,如果mn,即未知数的个数多于方程的个数，理论上有无穷个解，MATLAB寻求一个基本特解,左除法：x=Ab 最少元素解 广义逆法：x=pinv(A)*b 最小范数解,【例3-17】求解欠定系统。 a=fix(15*rand(2,3); b=fix(15*rand(2,1) p=ab,物理与电气工程学院,26,【例4-18】使用两种方法求解欠定系统。 a=1 1 1;1 1 -1; b=1 0 6; p=ab norm(p) a*p-b q=pinv(a)*b norm(q) a*q-b,27,第3节 数据分析 MATLAB对数据分析函数有两条约定： 输入参数:向量，则不论是行向量还是列向量，运算都是对向量整体进行的 输入参数:矩阵，则运算是按列 进行的,一、基本统计命令,mean(x) 均值或期望,median(x) 中位数，排序后中间的那个数(或中间两个数的均值),var(x) 方差,std(x) 标准差,物理与电气工程学院,28,sum(x) 求和,max(x) 最大值,min(x) 最小值,sum(x) 求和,prod(x) 求积,sx=cumsum(x) 逐元素累和,sx结构同x,px=cumprod(x) 逐元素累积,px结构同x,xs,xn=sort(x) 升序排序，xs为重新排序后数据，xn为xs各元素在x中的原始位置,y=sortrows(x, n) 矩阵x以第n列为基准，按递增顺序排列,物理与电气工程学院,29,【例4-19】基本统计函数的应用。 A=fix(10*randn(4,4); Amed=median(A) ; Amean=mean(A) Avar=var(A); Astd=std(A) %等同于std(A,0); Astd1=std(A,1) %std(A)2是对A方差的最佳无偏估计 Amax=max(A); Amin=min(A) Asum=sum(A); Aprod=prod(A) Acumsum=cumsum(A) Acumprod=cumprod(A) As,An=sort(A),物理与电气工程学院,30,有关正态分布的几个常用函数 fp=normpdf(x, mu, sigma) 正态分布密度函数在x处的值, mu=mean(X), sigma=std(X),P=normcdf(x0, mu, sigma) 正态分布随机变量小于等于x0的累积概率值,x0=norminv(P, mu, sigma) 正态分布逆累积函数，p=PXx0。与normcdf是互逆的,probability density 概率密度 cumulative distribution累计分布 inverse 相反的,物理与电气工程学院,31,附加：norm_1 x=-30:0.1:30; ynorm=normpdf(x,0,1) %标准正态分布 ynorm1=normpdf(x,3,2) %中心点由期望值确定，方差越大分布越分散 plot(x,ynorm,x,ynorm1),物理与电气工程学院,32,附加例norm_3:设X(3,22)，求PX3,P2X5 p1=normcdf(3,3,2) p2=normcdf(5,3,2) - normcdf(2,3,2),物理与电气工程学院,33,附加例norm_3:设X(3,22),确定c,使得PXc=0.612 c=norminv(0.612,3,2) p=normcdf(c,3,2) %用于验证,34,二、协方差矩阵和相关阵 C=cov(x) 如果x是向量，则返回C是该向量的方差,C是一个标量；,C=cov(x)如果x是矩阵，则每一列相当于一个变量，返回C是该矩阵的列与列之间的协方差矩阵：C是一个列数与x相同的对称方阵，主对角线元素是x对应列矢量的方差 c(i,j)=meanx(:,i)-mean(x(:,i)x(:,j)-mean(x(:,j),C=cov(x,y) x和y是同维向量或矩阵，相当于cov(x(:) y(:),既先将x和y按展开成列向量,再求这两个列向量的协方差,diag(cov(x) 如果x是矩阵，则C就是该矩阵每个列向量的方差,物理与电气工程学院,35,S=corrcoef(x) 计算x的相关系数矩阵S，将x的每一列视为一个变量，S(i , j)是i列与j列之间的相关系数。,对于输入矩阵x，相关系数矩阵S和协方差矩阵C之间的关系为：,S=corrcoef(x,y) 相当于corrcoef(x(:) y(:)，即先将x和y按展开成列向量,再求这两个列向量的相关系,S为对称矩阵，主对角线元素值为1,物理与电气工程学院,36,【例4-20】计算协方差和相关系数矩阵。 x=rand(10,3); y=rand(10,3); z=rand(1,10) cx=cov(x) cy=cov(y) cz=cov(z) vz=var(z) cxy=cov(x,y),物理与电气工程学院,37,三、微分、差分与梯度 1、微分及差分 计算数值向量或矩阵的数值差分。对于一个数值向量或矩阵F，diff(F)就是计算F(2:m,:)-F(1:m-1,:)，其调用格式为： Y=diff( F, n, dim),其中，F是向量或数组，n是差分阶数，dim是指定沿着数组的哪一维进行差分。,difference微分，差分 derivative导数,物理与电气工程学院,38,【例4-21】计算三维数组的差分。 a=rand(3,3,2) diff(a) %1阶差分，列内 diff(a,2)%第2阶差分，列内 diff(a,3) %第3阶差分，行内 diff(a,4) %第4阶差分，行内 diff(a,5)%第5阶差分，页间 diff(a,6)%第6阶差分，空,物理与电气工程学院,39,2、近似梯度 gradient函数用来进行数值梯度的计算，一般情况下，函数F(x,y,z,)的梯度定义为：,FX,FY=gradient(F,h1,h2) FX,FY,FZ=gradient(F,h1,h2,h3),FX、FY、FZ分别按矩阵的列、行、页方向；,h1,h2,h3指定点间隔,物理与电气工程学院,40,【例4-22】计算表达式xe(-x.2 - y.2)的梯度。 v = -2:0.2:2; % 生成-2到2间隔为0.2的自变量 x,y = meshgrid(v); % 产生数据网格 z = x .* exp(-x.2 - y.2); % 计算z px,py = gradient(z,.2,.2); % 求二维梯度 contour(x,y,z) % 绘制等高线 hold on % 保持绘图 quiver(x,y,px,py) % 绘制矢量图，抖动 hold off figure mesh(x,y,z),41,第4节 插值 插值(Interpolation)是在原始数据点之间按照一定的关系插入新的数据点，以便能较准确地分析数据的变化规律。,一、一维插值 一维插值就是对一维函数y=f(x)的数据进行插值: yi=interp1(x,y,xi,method),它是在“基准”数据的基础上，研究如何“平滑”地估算出“基准”数据间其他点的函数值。,原始数据点（x,y）,x为横坐标向量，y为纵坐标向量；xi 待插值点的横坐标，yi 为插值函数计算出的待插值点的纵坐标,物理与电气工程学院,42,X的数据必须按单调方式排列,如果y为矩阵，则插值将按照y的列向量进行,xi为指定的拟插值点的横坐标，yi是在xi指定位置计算出的插值结果,如果xi的某元素xi(i)超出x的定义范围，那么线性插值和最近邻插值方法将相应的yi(i)将取值为Nan。,物理与电气工程学院,43,method用于指定插值的方法，缺省时默认采用线性插值方法,物理与电气工程学院,44,【例4-23】一维插值函数插值方法的对比。 x=0:10; y=sin(x); xi=0:.25:10; % 将插值方法定义为单元数组 strmod=nearest,linear,spline,cubic % 将X轴标识定义为单元数组 strlb=(a) method=nearest, (b) method=linear,. (c) method=spline, (d) method=cubic; for i=1:4 yi=interp1(x,y,xi,strmodi); subplot(2,2,i),plot(x,y,ro,xi,yi,b) xlabel(strlb(i) end,物理与电气工程学院,45,spline三次样条插值函数 yi=spline(x,y,xi) 参数使用方法与interp1相似,并指定用三次样条方法,【例4-24】三次样条插值。 x0=0:10; y0=sin(x0); x=0:.25:10; y=spline(x0,y0,x); plot(x0,y0,or,x,y,k),46,二、二维插值 二维函数插值是基于一维函数插值的基本思想，它是对两个变量的函数z=f(x,y)进行插值，z为矩阵，由x,y确定的点上的值。 zi=interp2(x,y,z,xi,yi,method),输入参数为原始数据点(x,y,z)；x,y为两个独立的矩阵，它们也必须是按照单调方式排列的，z为函数矩阵，是由x,y确定的点上的值,xi为指定插值点横坐标的数值数组，yi是指定插值点纵坐标的数值数组。,method参数同一维插值。,物理与电气工程学院,47,【例4-25】二维插值四种方法的对比 x,y,z=peaks(7); % 生成双峰函数值,x,y数据间隔为1 mesh(x,y,z) % 绘制网格图 xi,yi=meshgrid(-3:0.2:3,-3:0.2:3); % 生成待插值的数据网格 z1=interp2(x,y,z,xi,yi,nearest); z2=interp2(x,y,z,xi,yi,linear); z3=interp2(x,y,z,xi,yi,spline); z4=interp2(x,y,z,xi,yi,cubic); figure,mesh(xi,yi,z1) figure,mesh(xi,yi,z2) figure,mesh(xi,yi,z3) figure,mesh(xi,yi,z4),物理与电气工程学院,48,三、多维插值 多维插值包括interp3和interpn函数,使用方法与interp2类似。 vi=interp3(x,y,z,v,xi,yi,zi,method),【例4-26】三维插值实例。 x,y,z,v = flow(10); slice(x,y,z,v,6 9.5,2,-2 .2) xi,yi,zi = meshgrid(0:.25:10, -3:.25:3, -3:.25:3); vi = interp3(x,y,z,v,xi,yi,zi); figure,slice(xi,yi,zi,vi,6 9.5,2,-2 .2); shading flat,物理与电气工程学院,49,第5节 快速傅利叶变换 离散序列的傅立叶变换（Discrete Fourier Transform,DFT）和傅立叶逆变换（Invense Discrete Fourier Transform,IDFT）在数学上定义和MATLAB软件中的定义基本相同。不同在于：MATLAB中的序列或向量元素下标从1开始，而不是数学上的从0开始。,物理与电气工程学院,50,MATALB提供了fft(内置函数)、iff、fft2、fftn、ifftn、fftshift、ifftshift等函数，用来计算矩阵的离散快速傅立叶变换。,在数据长度是2的幂次时，采用基-2算法进行计算，计算速度会显著增加，因此在使用时，尽量使用数据长度为2的幂次或者用零来填补数据。,快速傅里叶变换的结果为复数。,物理与电气工程学院,51,一、函数fft和ifft Y=fft(X) X为矢量傅里叶变换序列，变换序列长度与X等长;,Y=fft(X,n) 指定变换序列长度n，如果X的长度小于n，则用0来补足，如果X的长度大于n，则去掉X多出的部分。n可缺省时，视变换序列长度与X等长。,Y=fft(X,n,dim) 对指定维进行离散傅里叶变换,Y=fft(X,dim) n可缺省时，视变换序列长度与X等长,若X为矩阵则计算每列的傅立叶变换序列，,物理与电气工程学院,52,Y=fft(X,N),点域：N/2+1 N 1 N/2,频域：-Fs/2 -Fs/N 0 Fs/2-Fs/N,X的采样频率Fs,点频率=Fs/N,物理与电气工程学院,53,【例fft_test】fft实验 fs=102; %采样频率fs=1/ts,采样频率要大于信号最高频率的2倍 ts=1/fs; %采样间隔 t=0:ts:5; x=100*sin(2*pi*50*t)+randn(size(t); %加入高斯噪声正胘信号，频率为50Hz N=512; %变换序列点数，最好是2的幂 y=fft(x,N); Pf=fs/N; %点频率fs/N f=(0:N/2-1)*Pf; plot(f,abs(y(1:N/2); %幅频特性,物理与电气工程学院,54,55,【例4-27】产生一个正弦衰减曲线，进行快速傅里叶变换，并画出幅值图、 相位图、实部图和虚部图。 tp=0:1/6:2048; % 时域采样时间间隔 ts=1/6 fs=6; %抽样频率 % 根据奈奎斯特定理，抽样频率大于信号最高频 率的2倍 yt=sin( 4*pi*tp).*exp(-tp/80); % 信号最高频率2Hz yf=fft( yt ); % 快速傅里叶变换序列点数与yt长度相同 pf=fs / length( yt ); %变换序列点频率 lee=length(yt)/2 ; %频谱取段的长度 f=(0:lee-1).*pf; %频域轴坐标 ya=abs(yf(1:lee); % 幅值 yp=angle(yf(1:lee)*180/pi; % 相位(角度) plot(tp,yt) % 绘制正弦衰减曲线 figure, plot(f,ya) % 绘制FFT幅值曲线 figure, plot(f,yp) % 绘制FFT相位曲线,56,物理与电气工程学院,57,二、快速傅里叶变换的长度与运算速度 fff函数的变换点数n决定变换的速度,n为2的整数次幂，运算速度最快,n为合数，fft采用质因数分解的算法，速度取决于质因数的大小,n为质数，运算速度最慢,几个函数解释 isprime(x) 判断x是否质数 factor(x) 对x进行因式分解 cputime 计算机cpu在此时的时刻,物理与电气工程学院,58,【例3-