第2章-MATLAB的数值计算.ppt

第2章MATLAB的数值计算 所谓数值计算 就是指计算的表达式 变量中不得包含未经定义的自由变量 本章主要内容 2 1变量及其赋值2 2矩阵和数组的算术运算2 3线性方程组的求解2 4矩阵函数2 5多项式运算2 6练习 2 1变量及其赋值 在命令窗口可直接输入变量并赋值 无须声明其类型和维数 MATLAB将自动处理 MATLAB的变量名称由英文大小写52个字母 数字和下划线等符号组成 并且第一个字符必须是英文字母 变量和常量最长允许31个字符 须注意的是 在默认状态下 MATLAB对字母大小写敏感 MATLAB数据格式与显示形式 MATLAB既可用传统的十进制数表达数值 也可以用科学计数表达数值 用e代表10的指数形式 用i和j来代表虚数 MATLAB内部数据格式只有一种 是IEEE浮点标准的双精度二进制 64位 相应于十进制的16位有效数 范围为10 308 10 308 为了人机交互的友好性 数据输出显示格式有8种 可用菜单选项或format命令选择 MATLAB定义的特别变量及其意义 2 2矩阵和数组的算术运算 2 2 1矩阵的定义 MATLAB中的变量或常量都代表矩阵 矩阵变量中的元素用 中的数字来注明其下标 维数不同 数字个数也不同 数字之间用 隔开使用下标 用户可以单独给元素赋值 如x 2 1 7321 a 2 3 6等 实例 访问矩阵中的元素 x 0 pi 4 2 pi y sin x A x y 先构造矩阵A上面构造的A由列向量x和在x基础上生成的列向量y合并组成 由定义 x是一长度为9的列向量 因此这样构造的A是一个9 2的矩阵 b A 2 1 此命令访问了A的第二行第一列元素 创建特殊矩阵的函数 zeros 矩阵元素都是零ones 矩阵元素都是1注意 中的数字用于定义矩阵维数 内不是数字的情况参看Help rand 矩阵元素是0到1之间均匀分布的随机数randn 矩阵元素是正态分布的随机数 2 2 2矩阵运算 MATLAB提供矩阵算术运算有 加 减 乘 除 幂次方 转置 注意 对于矩阵来说 乘法有左乘右乘之分 A B时要求A的列数等于B的行数 即满足矩阵相乘的条件 MATLAB中 为求解线性方程组 引入了矩阵除法 MATLAB中的矩阵除法有左除 和右除 两种 以下是矩阵运算的合法命令形式 A B A B A B A B A B A p 此处p只能是标量 A必须为方阵 A 矩阵的左除和右除 Case1 方程为A x B若A非奇异 用左除解得x A B如A为方阵 则A 形式上相当于A的逆阵 但须注意 MATLAB中 求解线性方程组并不采用求逆的方法 Case2 方程为x A B若A非奇异 用右除解得x B A如A为方阵 A形式上相当于A的逆矩阵 2 2 3数组运算 和矩阵运算符比较 除了加减符号外 其余的数组运算符号均须加上一个 符号 即 加 减 乘 左除 右除 幂次方 数组运算是矩阵中元素对元素的操作 因而进行数组运算时 加 减 乘 除均要求A B有相同的维数 对加减运算而言 A和B中有一个是标量也可以 数组相除也有左除和右除之分 其中A B表示A的元素为被除数 而A B则表示A的元素是除数 数组运算A p中 如p是标量 则此运算是指对A中每个元素按幂次p进行乘方运算 但若数组运算A B中 A为标量 B为矩阵 则A B的含义为AB i 范数 范数 信号处理问题中 信号常表示为线性空间中的一个向量 而向量的长度是一个非常重要的量 利用上述的数组运算A p 向量X的长度的平方很容易表示为sum X 2 范数 norm 概念可容易地用来度量向量长度 MATLAB中设置了专用的norm函数norm X p sum X p 1 p 注意 这里用到了上面的sum 函数 常用的p值为1 2 p 2时返回X的二阶范 即通常意义下向量的欧氏长度 2 3线性方程组的求解 实际情况下 通常遇到的线性方程组的形式为A X B其中 系数矩阵A不必为方阵设其为m n维矩阵 1 m n 系统为超定 X A B 给出系统最小二乘解 2 m n 如矩阵A非奇异 X A B 给出系统的精确解 3 m n 系统为欠定 X A B 给出系统的一组通解 不具惟一性 例 均衡滤波器 图中是一个5抽头的横向滤波器 用作均衡滤波器 C1表示抽头因子构成的列向量 其输入为h n 设h n 长为5点h 0 1 0 2 0 9 0 3 0 15 理想情况下 均衡滤波器输出y n 的理想值应为y 000010000 长度为9点 我们的任务是确定抽头因子的数值 例 均衡滤波器的最小均方误差 MMSE 解与迫零 Zeroforcing 解 MMSE解 1 求出h的卷积矩阵A convmtx h 5 2 令y1 y 左除求解超定方程组A C1 y1 即得MMSE解 即最小二乘解C1 Zeroforcing解 1 去掉A中的上下各两行 构成方阵B A 3 7 1 5 2 令y2 00100 左除求解B C2 y2 即得Zeroforcing解C2 运用范数概念比较两个解 对于得到的C1 C2 均衡滤波器输出分别为y11 conv h C1 y22 conv h C2 长度相同 均为9点 y11 y22与理想输出y之间的误差分别为y y11和y y22 这两个误差向量的二阶范数就是方均根误差 也可用均方误差即二阶范的平方来衡量 MinimumMSE解优于Zeroforcing解 Why 2 4矩阵函数 本节给出了一些重要的矩阵函数 在信号处理中 这些函数常会用到 对初学者而言 本节内容如感困难 可略过不学 不影响后续内容的学习 矩阵范数 秩与条件数 行列式与逆矩阵 三种矩阵分解 特征值与特征向量 奇值分解 SVD 2 5多项式运算 多项式在工程中有着广泛的应用 如实验数据的多项式拟合表达等 2 5 1多项式的定义 1 直接定义令a x 代表一个如下的多项式在MATLAB中 这个多项式可以用它的系数向量A a 1 a 2 a n a n 1 来表示 2 用命令poly创建如果A是矩阵 则poly A 将创建A矩阵的特征多项式 如果A是向量 则poly A 将创建以A中各元素为根的多项式 2 5 2多项式的四则运算 1 多项式相乘 以两个多项式为例 conv函数本用于两个有限长序列的卷积 但卷积和的运算也符合多项式相乘规则 因此在命令窗口键入 a 2 4 6 8 b 3 6 9 c conv a b 后 可得c 624609610272而 2 多项式相加 MATLAB规定 只有长度相同即维数相同的向量才能相加 因此 必须把短的向量前面补以若干个零元素 也即对多项式高次项系数以零表示 才能用MATLAB的矩阵加法运算符 可以编写一个子程序polyadd m 以使计算机自动完成两个多项式符合相加的条件 关于子程序的编写 下一章将会学习到 functiony polyadd x1 x2 n1 length x1 n2 length x2 ifn1 n2x2 zeros 1 n1 n2 x2 elseifn1 n2x1 zeros 1 n2 n1 x1 end y x1 x2 3 多项式相除 相除是相乘的逆运算 不能除尽时会有余子式产生 故在MATLAB中 多项式相除后会返回一项余子式 能够除尽时 余子式全为零 函数形式为 q r deconv c a 其中q是商 r是余子式 注意deconv的原意是反卷积 工程中也称之为逆滤波 但由于deconv允许出现余子式 这个函数不能保证得到满意的逆滤波解 2 5 3多项式求导 求根和求值 1 多项式求导数 polyder 函数形式 e polyder c 对上面的c x 使用这一函数得到e 3096180192102表示 2 多项式求根 函数形式 ra roots a 其中a表示多项式a x 注意 前面已经学习过根据多项式的根创建多项式系数的函数poly 实际上 这是roots的逆运算 即a poly ra 这里ra是一个由根构成的向量 3 多项式求值 polyval 将多项式a中的自变量x赋予值xv时 该多项式的值可用F polyval a xv 求得 其中xv可以是复数 也可以是矩阵或数组 对于后者 polyval对输入变元的运算结果也是矩阵或数组 这一性质可以用于求线性系统的频率特性 例2 1 LTI系统频率响应 求此系统的频率响应并画出幅频与相频特性 解 设a为该系统分母系数向量 b为系统分子系数向量 即 a 1234 b 123 该系统的频率特性为程序中使用polyval函数求出了系统幅频特性和相频特性在频率数组w各个取值处的数值 MATLAB程序见教材p32 33 频率特性图见图2 2 例2 1设有LTI系统的系统函数为 2 5 4多项式拟合 函数 p polyfit x y n 用于在最小均方误差准则下对数据进行n阶多项式曲线拟合 其中x y是已知的N个数据点坐标向量 其长度均为N 现欲用一个n阶多项式来表示x y之间的函数关系 待定的是该n阶多项式的系数 多项式阶数n越大 拟合程度越好 设p为用于拟合的多项式系数个数 p n 1 例2 2 测量放大器3dB通频带的测定 参见pp33 35 2 6练习 命令窗口中 键入 demo 出现MATLAB的demo演示窗口 然后选择 1 MATLAB下的Matric