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

第4章 数值运算基础 MATLAB的计算包括数值计算 和符号计算两类，本章将带大家学 习数值计算部分。其中将主要学习 与我们专业密切相关的多项式、方 程组求解、数据分析和数字信号处 理的快速傅里叶变换。 1物理与电气工程学院 第1节 多项式polynomial MATLAB用行向量表示多项式。将 多项式的系数按降幂次序存放在行向量 中。如： 的系数行向量为 注意：多项式中缺少的幂次系数一 定要用“0”补齐。 2物理与电气工程学院 一、创建多项式 1、系数矢量直接输入法 在命令窗口直接输入多式的系数向量 【例4-1】输入系数矢量，创建多项式 x3-4*x2+3*x+2。 p=1 -4 3 2 poly2sym(p) %将矢量P表示为多项式的手写形式 Polynomial coefficient vector to symbolic polynominal 3物理与电气工程学院 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) %由特征值函数求得的特征值 4 3、根矢量创建 p=poly(A) A为待求多项式的根 矢量，则返回值将 是对应多项式的系 数行矢量，该多项 式的根为矢量A。此 时p=poly(A)与 A=roots(p)互逆。系 统定义P0=1。 Ax1 x2 x3 p=1 p1 p2 p3 5 【例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对应的同维 向量或矩阵。 6 【例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) 7物理与电气工程学院 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) 8物理与电气工程学院 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)9物理与电气工程学院 多项式的除法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) 10 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) 11物理与电气工程学院 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。 12物理与电气工程学院 如果分母多项式a(s)不含重根，则两个多项式 可写成如下形式： 其中，pi称为极点，ri称为留数，k(s)是直项 。 如果b的次数低于a的次数，则直项为空。 如果分母多项式a(s)含m重根p(j)=p(j+m-1),则 这m项应展开为 13 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) %出现重根 笔算结果=？ 14物理与电气工程学院 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进行 多项式的曲线拟合，最后给出拟合的多项式 系数。 15物理与电气工程学院 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=(x-u1)/u2替代x，其中 mu=u1 u2,u1=mean(x),u2=std(x) 16物理与电气工程学院 【例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内拟合曲线拟近原 曲线，而在区间以外的曲线误差较大 17物理与电气工程学院 第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的列数。 18物理与电气工程学院 根据矩阵A的结构(m, n），可以将方 程分为以下3类： m=n 方阵系统，可偿试求精确解 mn 超定系统，可偿试求最小二乘解 mn,则方程没有精确解,此时方程组为超定 方程组。一般采用最小二乘法。 左除法：x=Ab 建立在奇异值分解基 础之上 广义逆法：x=pinv(A)*b 速度较快， 可靠性较差一些 实验中数据的曲线拟合就是就可以解超定方程组 的方法来解。一般情况下需要将非线性问题转换 为线性问题来解决 22物理与电气工程学院 一组实验数据，时间t和测量数据y，如下表所示： itexp(-ti)yi 10.010.82 20.30.7400.72 30.80.4490.63 41.10.3320.60 51.60.2020.55 62.20.1110.5 认为x和y有上式的关系式 ，则由6组实验数据就可形 成超定方程组，就可求解 出C1和C2。使得将t代入函 数得到的值与实际y值之间 差值的最小平方和最小 如何以矩阵的 形式表示？ 23物理与电气工程学院 【例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*) 24物理与电气工程学院 三 、欠定系统 程组Ax=b，A为mn矩阵,如果mc=0.612 c=norminv(0.612,3,2) p=normcdf(c,3,2) %用于验证 33物理与电气工程学院 二、协方差矩阵和相关阵 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就是该矩阵每个列向量的方 差 34 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 35物理与电气工程学院 【例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) 36物理与电气工程学院 三、微分、差分与梯度 1、微分及差分 计算数值向量或矩阵的数值差分。对于一个数值向 量或矩阵F，diff(F)就是计算F(2:m,:)-F(1:m-1,:)，其调 用格式为： Y=diff( F, n, dim) 其中，F是向量或数组，n是差分阶数，dim是指定沿着数组的 哪一维进行差分。 difference微分，差分 derivative导数 37物理与电气工程学院 【例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阶差分，空 38物理与电气工程学院 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指定点间隔 39物理与电气工程学院 【例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) 40物理与电气工程学院 第4节 插值 插值(Interpolation)是在原始数据点之间按照一 定的关系插入新的数据点，以便能较准确地分 析数据的变化规律。 一、一维插值 一维插值就是对一维函数y=f(x)的数据进行插 值: yi=interp1(x,y,xi,method) 它是在“基准”数据 的基础上，研究如何“平滑”地估算出“基准 ”数据间其他点的函数值。 原始数据点（x,y）,x为横坐标向量，y为纵坐标 向量；xi 待插值点的横坐标，yi 为插值函数计算出的 待插值点的纵坐标 41 X的数据必须按单调方式排列 如果y为矩阵，则插值将按照y的列 向量进行 xi为指定的拟插值点的横坐标，yi是 在xi指定位置计算出的插值结果 如果xi的某元素xi(i)超出x的定义范 围，那么线性插值和最近邻插值方 法将相应的yi(i)将取值为Nan。 42物理与电气工程学院 method用于指定插值的方法，缺省时默 认采用线性插值方法 Method含义特点和用途 Linear线线性插值值用作连连接图图上的数据点 Cubic 三次多项项式 插值值 精度高、速度慢，平滑性好 spline 三次样样条 插值值 精度最高、速度最慢，最平滑 nearst 最近领领域 插值值 精度低、速度快 43物理与电气工程学院 【例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 44物理与电气工程学院 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) 45物理与电气工程学院 二、二维插值 二维函数插值是基于一维函数插值的基 本思想，它是对两个变量的函数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参数同一维插值。 46 【例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) 47物理与电气工程学院 三、多维插值 多维插值包括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 48物理与电气工程学院 第5节 快速傅利叶变换 离散序列的傅立叶变换（Discrete Fourier Transform,DFT）和傅立叶逆变换（Invense Discrete Fourier Transform,IDFT）在数学上 定义和MATLAB软件中的定义基本相同。不 同在于：MATLAB中的序列或向量元素下标 从1开始，而不是数学上的从0开始。 49物理与电气工程学院 MATALB提供了fft(内置函数)、iff 、fft2、fftn、ifftn、fftshift、ifftshift等函 数，用来计算矩阵的离散快速傅立叶变换 。 在数据长度是2的幂次时，采用基-2算法 进行计算，计算速度会显著增加，因此在 使用时，尽量使用数据长度为2的幂次或 者用零来填补数据。 快速傅里叶变换的结果为复数。 50物理与电气工程学院 一、函数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为矩阵则计算每列的傅 立叶变换序列， 51物理与电气工程学院 Y=fft(X,N) 1 N N/2 N/2+1 点域：N/2+1 N 1 N/2 频域：-Fs/2 -Fs/N 0 Fs/2-Fs/N X的采样频率Fs 点频率=Fs/N 52物理与电气工程学院 【例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); %幅频特性 53物理与电气工程学院 54物理与电气工程学院 【例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相位曲线 55 56 二、快速傅里叶变换的长度与运算速度 fff函数的变换点数n决定变换的速度 n为2的整数次幂，运算速度最快 n为合数，fft采用质因数分解的算法，速度取决于质因数 的大小 n为质数，运算速度最慢 几个函数解释 isprime(x) 判断x是否质数 factor(x) 对x进行因式分解 cputime 计算机cpu在此时的时刻 57物理与电气工程学院 【例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-t t=cputime; y=fft(x,216); e=cputime-t t=cputime;