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

第4章 数值运算基础 MATLAB的的计计算算包包括括数数值值计计算算和和符符号号计计算算两两类类，本本章章将将带带大大家家学学习习数数值值计计算算部部分分。其其中中将将主主要要学学习习与与我我们们专专业业密密切切相相关关的的多多项项式式、方方程程组组求求解解、数数据据分分析析和和数数字字信信号号处处理的快速傅里叶变换。理的快速傅里叶变换。1物理与电气工程学院第第1 1节多项式节多项式polynomialpolynomial MATLAB用用行行向向量量表表示示多多项项式式。将将多多项项式式的的系系数数按按降降幂幂次次序序存存放放在在行行向向量量中。如：中。如：的系数行向量为的系数行向量为 注意：多项式中缺少的幂次系数一注意：多项式中缺少的幂次系数一定要用定要用“0”补齐。补齐。2物理与电气工程学院一、创建多项式一、创建多项式1、系数矢量直接输入法、系数矢量直接输入法在命令窗口直接输入多式的系数向量在命令窗口直接输入多式的系数向量 【例【例4-1】输入系数矢量，创建多项式】输入系数矢量，创建多项式 x3-4*x2+3*x+2。p=1 -4 3 2poly2sym(p)%将矢量将矢量P表示为多项式的手写形式表示为多项式的手写形式Polynomial coefficient vector to symbolic polynominal3物理与电气工程学院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)%由特征值函数求得的特征值由特征值函数求得的特征值43、根矢量创建、根矢量创建 p=poly(A)A为为待待求求多多项项式式的的根根矢矢量量，则则返返回回值值将将是是对对应应多多项项式式的的系系数数行行矢矢量量，该该多多项项式式的的根根为为矢矢量量A。此此时时p=poly(A)与与A=roots(p)互互逆逆。系系统定义统定义P0=1。Ax1 x2 x3 p=1 p1 p2 p35【例【例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)104、多项式微积分、多项式微积分polyder(p)返回多项式系数向量返回多项式系数向量p p 的导数的导数【例【例4-9】计算多项式】计算多项式3x4-5x3+2x2-6x+10的微分。的微分。p=3-5 2-6 10;dp=polyder(p)poly2sym(dp)polyint(p)返回多项式系数返回多项式系数p p 的积分的积分【例【例4-10】计算多项式】计算多项式12x3-15x2+4x-6的积分。的积分。p=12-15 4-6;ip=polyint(p)poly2sym(ip)11物理与电气工程学院5、多项式的部分分式展开、多项式的部分分式展开MATLABMATLAB提供了提供了residueresidue命令来执行部分分式命令来执行部分分式展开或多项式系数之间的转换。该命令通展开或多项式系数之间的转换。该命令通常用于信号与控制领域中。格式如下：常用于信号与控制领域中。格式如下：r,p,k=residue(b,a)该命令是求多项式之比该命令是求多项式之比b(s)/a(sb(s)/a(s)的部分分的部分分式展开，返回式展开，返回留数留数r r、极点极点p p和和直项向量直项向量k k。a a和和b b分别是分母和分子多项式的系数向量分别是分母和分子多项式的系数向量b,a=residue(r,p,k)上一条语句的上一条语句的逆过程逆过程，只是分母多项式系，只是分母多项式系数以归一化形式：最高次项系数数以归一化形式：最高次项系数a(1)a(1)为为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 1r2,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*x0p,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)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.50 2.5内内拟合曲线拟近原拟合曲线拟近原曲线，曲线，而在区间以而在区间以外外的曲线误差较大的曲线误差较大17物理与电气工程学院第第2 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)yi10.010.8220.30.7400.7230.80.4490.6341.10.3320.6051.60.2020.5562.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=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*)24物理与电气工程学院三三、欠定系统、欠定系统程组程组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=ab25物理与电气工程学院【例例4-18】使用两种方法求解欠定系统。使用两种方法求解欠定系统。a=1 1 1;1 1-1;b=1 0 6;p=abnorm(p)a*p-bq=pinv(a)*bnorm(q)a*q-b26物理与电气工程学院第3节数据分析MATLAB对数据分析函数有两条约定：对数据分析函数有两条约定：输入参数输入参数:向量，则不论是行向量还是列向量，向量，则不论是行向量还是列向量，运算都是对运算都是对向量整体向量整体进行的进行的输入参数输入参数:矩阵，则运算是按矩阵，则运算是按列列 进行的进行的一、基本统计命令一、基本统计命令mean(x)均值或期望均值或期望median(x)中位数，排序后中间的那个数中位数，排序后中间的那个数(或中间或中间两个数的均值两个数的均值)var(x)方差方差std(x)标准差标准差27sum(x)求和求和max(x)最大值最大值min(x)最小值最小值sum(x)求和求和prod(x)求积求积sx=cumsum(x)逐元素累和逐元素累和,sx结构同结构同xpx=cumprod(x)逐元素累积逐元素累积,px结构同结构同xxs,xn=sort(x)升序排序，升序排序，xs为为重新排重新排序序后数据，后数据，xn为为xs各元素在各元素在x中的中的原始位置原始位置y=sortrows(x,n)矩阵矩阵x以第以第n列为基准，按列为基准，按递递增增顺序排列顺序排列28物理与电气工程学院【例【例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)29物理与电气工程学院有关有关正态分布正态分布的几个常用函数的几个常用函数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 相反的相反的30物理与电气工程学院附加：附加：norm_1x=-30:0.1:30;ynorm=normpdf(x,0,1)%标准正态分布标准正态分布ynorm1=normpdf(x,3,2)%中心点由期望值确定，方差越大分布越分散中心点由期望值确定，方差越大分布越分散plot(x,ynorm,x,ynorm1)31物理与电气工程学院附加例附加例norm_3:设设X(3,22)，求，求PX3,P2Xc=0.612c=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就是该矩阵每个就是该矩阵每个列列向量的方差向量的方差34S=corrcoef(x)计算计算x的相关系数矩阵的相关系数矩阵S，将将x的每一的每一列列视为一个变量，视为一个变量，S(i,j)是是i列列与与j列之间的相关系数。列之间的相关系数。对于输入矩阵对于输入矩阵x，相关系数矩阵，相关系数矩阵S和协方和协方差矩阵差矩阵C之间的关系为：之间的关系为：S=corrcoef(x,y)相当于相当于corrcoef(x(:)y(:)，即先将即先将x和和y按展开成按展开成列列向量向量,再求这两个列再求这两个列向量的相关系向量的相关系S为对称矩阵，为对称矩阵，主对角线元素值为主对角线元素值为135物理与电气工程学院【例【例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);%计算计算zpx,py=gradient(z,.2,.2);%求二维梯度求二维梯度contour(x,y,z)%绘制等高线绘制等高线hold on%保持绘图保持绘图quiver(x,y,px,py)%绘制矢量图，绘制矢量图，抖动抖动hold offfiguremesh(x,y,z)40物理与电气工程学院第第4节插值节插值插值插值(Interpolation)是在原始数据点之间按是在原始数据点之间按照一定的关系插入新的数据点，以便能较准确照一定的关系插入新的数据点，以便能较准确地分析数据的变化规律。地分析数据的变化规律。一、一维插值一、一维插值一维插值就是对一维函数一维插值就是对一维函数y=f(x)的数据进行插的数据进行插值值:yi=interp1(x,y,xi,method)它是在它是在“基准基准”数据数据的基础上，研究如何的基础上，研究如何“平滑平滑”地估算出地估算出“基准基准”数据间其他点的函数值。数据间其他点的函数值。原始数据点（原始数据点（x,yx,y）,x,x为横坐标向量，为横坐标向量，y y为纵坐标为纵坐标向量；向量；xi 待插值点的横坐标，待插值点的横坐标，yi 为插值函数计算出的为插值函数计算出的待插值点的纵坐标待插值点的纵坐标41X的数据必须按的数据必须按单调单调方式排列方式排列如果如果y为为矩阵矩阵，则插值将按照，则插值将按照y的的列列向向量进行量进行 xi为指定的拟插值点的横坐标，为指定的拟插值点的横坐标，yi是是在在xi指定位置计算出的插值结果指定位置计算出的插值结果如果如果xi的某元素的某元素xi(i)超出超出x的定义范的定义范围围，那么线性插值和最近邻插值方，那么线性插值和最近邻插值方法将相应的法将相应的yi(i)将取值为将取值为Nan。42物理与电气工程学院method用于指定插值的方法，缺省时默用于指定插值的方法，缺省时默认采用线性插值方法认采用线性插值方法MethodMethod含义含义特点和用途特点和用途Linear线线性插性插值值用作用作连连接接图图上的数据点上的数据点Cubic三次多三次多项项式式插插值值精度高、速度慢，平滑性好精度高、速度慢，平滑性好spline三次三次样样条条插插值值精度最高、速度最慢，最平滑精度最高、速度最慢，最平滑nearst最近最近领领域域插插值值精度低、速度快精度低、速度快43物理与电气工程学院【例【例4-23】一维插值函数插值方法的对比。】一维插值函数插值方法的对比。x=0:10;y=sin(x);xi=0:.25:10;%将插值方法定义为单元数组将插值方法定义为单元数组strmod=nearest,linear,spline,cubic%将将X 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)end44物理与电气工程学院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数据间隔为数据间隔为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)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 flat48物理与电气工程学院第第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和和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为矩阵为矩阵则计算每则计算每列列的傅的傅立叶变换序列，立叶变换序列，51物理与电气工程学院Y=fft(X,N)1NN/2N/2+1点域：点域：N/2+1 N 1 N/2频域：频域：-Fs/2 -Fs/N 0 Fs/2-Fs/NX的的采样频率采样频率Fs点频率点频率=Fs/N52物理与电气工程学院【例【例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);%幅频特性53物理与电气工程学院54物理与电气工程学院【例【例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相位曲线相位曲线5556二、快速傅里叶变换的长度与运算速度二、快速傅里叶变换的长度与运算速度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-tt=cputime;y=fft(x,66000);e=cputime-tt=cputime;y=fft(x,65535);e=cputime-t58物理与电气工程学院三、三、fftshift和和ifftshift函数函数fftshift用于把傅立叶变换结果（频域数据）中用于把傅立叶变换结果（频域数据）中的直流分量（频率为的直流分量（频率为0Hz处的值）移到中间位置。处的值）移到中间位置。Y=fftshift(X)X是傅立叶变换结果是傅立叶变换结果如果如果X是向量，则交换是向量，则交换X的左右半边的左右半边如果如果X是矩阵，则交换其一、三象限和二、是矩阵，则交换其一、三象限和二、四象限四象限 ifftshift相当于相当于fftshift函数的操作逆转，用法函数的操作逆转，用法相同相同 59物理与电气工程学院【例【例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);60物理与电气工程学院期中练习期中练习11、生成、生成5行行6列的矩阵，矩阵元素为列的矩阵，矩阵元素为1到到61之间之间的偶数，要求元素数值按单下标顺序递增。的偶数，要求元素数值按单下标顺序递增。2、找出矩阵中大于、找出矩阵中大于9的元素的行列标的元素的行列标号，并有逻辑数组将其中元素值为号，并有逻辑数组将其中元素值为6的换为的换为8。3、将矩阵、将矩阵A、B、C组成一个新的矩阵，三个矩组成一个新的矩阵，三个矩阵分别为新矩阵的第阵分别为新矩阵的第1、2、3行。行。62物理与电气工程学院期中练习期中练习24、绘制一个周期为、绘制一个周期为3的三角波。的三角波。5、绘制分离饼图，其中、绘制分离饼图，其中5.6对应对应的扇区分离。的扇区分离。6、用阶梯图表现函数：、用阶梯图表现函数：7、在一个图形窗口的两个子窗口