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

物理与电气工程学院,1,第4章数值运算基础,MATLAB的计算包括数值计算和符号计算两类，本章将带大家学习数值计算部分。其中将主要学习与我们专业密切相关的多项式、方程组求解、数据分析和数字信号处理的快速傅里叶变换。,物理与电气工程学院,2,第1节多项式polynomialMATLAB用行向量表示多项式。将多项式的系数按降幂次序存放在行向量中。如：,的系数行向量为,注意：多项式中缺少的幂次系数一定要用“0”补齐。,物理与电气工程学院,3,一、创建多项式1、系数矢量直接输入法在命令窗口直接输入多式的系数向量,【例4-1】输入系数矢量，创建多项式x3-4*x2+3*x+2。p=1-432poly2sym(p)%将矢量P表示为多项式的手写形式,Polynomialcoefficientvectortosymbolicpolynominal,4,2、方阵特征多项式输入法p=poly(A)若A为nn的矩阵，则返回值P将是一个含有n+1个元素的行向量，也就是该矩阵特征多项式的系数,【例4-2】求矩阵123;456;780的特征多项式系数，并转换为多项式形式。a=123;456;780;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。,Ax1x2x3,p=1p1p2p3,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=321;polyval(p,579),y=polyvalm(p,x)将矩阵整体（而不是矩阵元素）作为自变量进行计算。p是多项式的系数向量；相当于用矩阵x代替多项式的变量对矩阵而不是对数组进行计算，x必须是方阵,【例4-5】求多项式3x2+2x+1对于矩阵25;79的值p=321;polyvalm(p,25;79),物理与电气工程学院,8,2、求多项式的根格式：C=roots(p)p为多项式的系数矢量，C为函数返回多项式的根矢量如果C为复数，则必成对出现。,【例4-6】分别用两种方法求多项式x5-5x4+3x3-6x2+4x-10的根。a=1-53-64-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-53-42;b=12-53;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-1230-3633-226;b=12-53;q,r=deconv(c,b),物理与电气工程学院,11,4、多项式微积分polyder(p)返回多项式系数向量p的导数,【例4-9】计算多项式3x4-5x3+2x2-6x+10的微分。p=3-52-610;dp=polyder(p)poly2sym(dp),polyint(p)返回多项式系数p的积分,【例4-10】计算多项式12x3-15x2+4x-6的积分。p=12-154-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=-4083;b=53-27;r,p,k=residue(b,a)b1,a1=residue(r,p,k)%分母最高次项归1r2,p2,k2=residue(11,1-21)%出现重根,笔算结果=？,物理与电气工程学院,15,6、多项式拟合对于给定的一组数据(xi,yi)，i=1,2,n,如果要采用多项式模型对数据组进行描述，形成如多项式y(x)=f(x,p)=p1xn+p2xn-1+p3xn-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=10.512,则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=u1u2,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),可以看出，拟合区间02.5内拟合曲线拟近原曲线，,而在区间以外的曲线误差较大,第2节线性代数给定两个矩阵A和B,求X的解,使得:AX=BXA=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,则方程没有精确解,此时方程组为超定方程组。一般采用最小二乘法。,左除法：x=Ab建立在奇异值分解基础之上,广义逆法：x=pinv(A)*b速度较快，可靠性较差一些,实验中数据的曲线拟合就是就可以解超定方程组的方法来解。一般情况下需要将非线性问题转换为线性问题来解决,物理与电气工程学院,23,一组实验数据，时间t和测量数据y，如下表所示：,认为x和y有上式的关系式，则由6组实验数据就可形成超定方程组，就可求解出C1和C2。使得将t代入函数得到的值与实际y值之间差值的最小平方和最小,如何以矩阵的形式表示？,物理与电气工程学院,24,【例4-16】求表中数据的最小二乘解。t=00.30.81.11.62.3;y=0.820.720.630.600.550.50;e=ones(size(t)exp(-t)%6组实验数据变换得到的c=eycp=pinv(e)*yt1=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=111;11-1;b=106;p=abnorm(p)a*p-bq=pinv(a)*bnorm(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是互逆的,probabilitydensity概率密度cumulativedistribution累计分布inverse相反的,物理与电气工程学院,31,附加：norm_1x=-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,P2Xc=0.612c=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);%计算zpx,py=gradient(z,.2,.2);%求二维梯度contour(x,y,z)%绘制等高线holdon%保持绘图quiver(x,y,px,py)%绘制矢量图，抖动holdofffiguremesh(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;fori=1:4yi=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数据间隔为1mesh(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,69.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,69.5,2,-2.2);shadingflat,物理与电气工程学院,49,第5节快速傅利叶变换离散序列的傅立叶变换（DiscreteFourierTransform,DFT）和傅立叶逆变换（InvenseDiscreteFourierTransform,IDFT）在数学上定义和MATLAB软件中的定义基本相同。不同在于：MATLAB中的序列或向量元素下标从1开始，而不是数学上的从0开始。,物理与电气工程学院,50,MATALB提供了fft(内置函数)、iff、fft2、fftn、ifftn、fftshift、ifftshift等函数，用来计算矩阵的离散快速傅立叶变换。,在数据长度是2的幂次时，采用基-2算法进行计算，计算速度会显著增加，因此在使用时，尽量使用数据长度为2的幂次或者用零来填补数据。,快速傅里叶变换的结果为复数。,物理与电气工程学院,51,一、函数fft和ifftY=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+1N1N/2,频域：-Fs/2-Fs/N0Fs/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);%加入高斯噪声正胘信号，频率为50HzN=512;%变换序列点数，最好是2的幂y=fft(x,N);Pf=fs/N;%点频率fs/Nf=(0:N/2-1)*Pf;plot(f,abs(y(1:N/2);%幅频特性,物理与电气工程学院,54,55,【例4-27】产生一个正弦衰减曲线，进行快速傅里叶变换，并画出幅值图、相位图、实部图和虚部图。tp=0:1/6:2048;%时域采样时间间隔ts=1/6fs=6;%抽样频率%根据奈奎斯特定理，抽样频率大于信号最高频率的2倍yt=sin(4*pi*tp).*exp(-tp/80);%信号最高频率2Hzyf=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-28】比较快速傅里叶变换的长度与运算速度的关系x=rand(70000,1);isprime(65539)%判断65539是否为质数216%计算2的16次方factor(66000)%计算66000的因数分解factor(65535)t=cputime;y=fft(x,65539);e=cputime-tt=cputime;y=fft(x,216);e=cputime-tt=cputime;y=fft(x,66000);e=cputime-tt=cputime;y=fft(x,65535);e=cputime-t,物理与电气工程学院,59,三、fftshift和ifftshift函数fftshift用于把傅立叶变换结果（频域数据）中的直流分量（频率为0Hz处的值）移到中间位置。Y=fftshift(X),X是傅立叶变换结果如果X是向量，则交换X的左右半边如果X是矩阵，则交换其一、三象限和二、四象限ifftshift相当于fftshift函数的操作逆转，用法相同,物理与电气工程学院,60,【例fft_shift】变换结果频谱直流分量移位fs=102;ts=1/fs;%采样间隔;t=0:ts:5;x=100*sin(2*pi*50*t)+randn(size(t);N=1024;%最好是2的幂y=fft(x,N);f=(0:N/2-1)*fs/N;%点频率fs/Nfigure(1);plot(f,abs(y(1:N/2);figure(2);F=(-N/2:(N/2-1)*fs/N;%频率轴坐标值plot(F,abs(y);%请判断频谱的正误figure(3);Y1=fftshift(y);%变换结果频谱直流分量移位F1=(-N/2:(N/2-1)*fs/N;plot(F1,abs(Y1);,物理与电气工程学院,62,期中练习1,1、生成5行6列的矩阵，矩阵元素为1到61之间的偶数，要求元素数值按单下标顺序递增。,2、找出矩阵中大于9的元素的行列标号，并有逻辑数组将其中元素值为6的换为8。,3、将矩阵A、B、C组成一个新的矩阵，三个矩阵分别为新矩阵的第1、2、3行。,物理与电气工程学院,63,期中练习2,4、绘制一个周期为3的三角波。,5、绘制分离饼图，其中5.6对应的扇区分离。,6、用阶梯图表现函数：,7、在一个图形窗口的两个子窗口内分别绘制网格图和表面图。,物理与电气工程学院,64,期中练习3,8、绘制曲线，并用“*”标出y最大值的点。,9、用多项式函数求在x=0.9时的值。,10、x和y的实验测试数据如下表，用多项式似合x和y的二次关系式，并求出拟合值与原始值之间的差值平方均值。,物理与电气工程学院,65,期中练习4,11、求下列方程组的解：,