返回 相似
第 1 页 / 共 65 页
亲,该文档支持全文预览,点击此可以查看更多
资源描述:
第4章 数值运算基础 MATLAB的计算包括数值计算和符号计算两类,本章将带大家学习数值计算部分。其中将主要学习与我们专业密切相关的多项式、方程组求解、数据分析和数字信号处理的快速傅里叶变换。1物理与电气工程学院第1节 多项式polynomial  MATLAB用行向量表示多项式。将多项式的系数按降幂次序存放在行向量中。如:的系数行向量为 注意:多项式中缺少的幂次系数一定要用“0”补齐。 2物理与电气工程学院一、创建多项式1、系数矢量直接输入法 在命令窗口直接输入多式的系数向量 【例4-1】输入系数矢量,创建多项式 x^3-4*x^2+3*x+2。p=[1 -4 3 2]poly2sym(p) %将矢量P表示为多项式的手写形式Polynomial coefficient vector to symbolic polynominal3物理与电气工程学院2、方阵特征多项式输入法p=poly(A)   若A为n×n的矩阵,则返回值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。A=[x1 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】求多项式3x^2+2x+1在5、7和9处的值。p = [3 2 1];polyval(p,[5 7 9])y=polyvalm(p,x) 将矩阵整体(而不是矩阵元素)作为自变量进行计算。p是多项式的系数向量;相当于用矩阵x代替多项式的变量对矩阵而不是对数组进行计算,x必须是方阵【例4-5】求多项式3x^2+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】分别用两种方法求多项式    x^5-5x^4+3x^3-6x^2+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】计算两多项式x^4-5x^3+3x^2-4x+2和x^3+2x^2-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中求得的乘积被x^3+2x^2-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 的导数【例4-9】计算多项式3x^4-5x^3+2x^2-6x+10的微分。 p=[3 -5 2 -6 10]; dp=polyder(p) poly2sym(dp)polyint(p) 返回多项式系数p 的积分【例4-10】计算多项式12x^3-15x^2+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) 5x^3+3x^2-2x+7【例3-11】两多项式的比为—— = ————————, a(x) -4x^3+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的。省略部分。标从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是否为质数2^16 % 计算2的16次方factor(66000) % 计算66000的因数分解factor(65535) t=cputime; y=fft(x,65539); e=cputime-t t=cputime; y=fft(x,2^16); 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、在一个图形窗口的两个子窗口内分别绘制网格图和表面图。63物理与电气工程学院期中练习38、绘制曲线 ,并用“*”标出y最大值的点。9、用多项式函数 求在x=0.9时的值。10、x和 y的实验测试数据如下表,用多项式似合x和 y的二次关系式,并求出拟合值与原始值之间的差值平方均值。x0.10.20.30.40.50.6Y1.32.53.97.29.111.364物理与电气工程学院期中练习411、求下列方程组的解:65物理与电气工程学院
展开阅读全文

电脑版 |天天文库版权所有
备案: 粤ICP备19057495号