《实用大众线性代数第七章.pdf》由会员分享，可在线阅读，更多相关《实用大众线性代数第七章.pdf（20页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、 第第 7 章章 在科技及工程中的应用实例在科技及工程中的应用实例 . 1 7.1 由拉压杆组成的桁架结构 . 1 7.2 格型梯形滤波器系统函数的推导 . 1 7.3 计算频谱用的 DFT 矩阵 . 2 7.4 显示器色彩制式转换问题 . 4 7.5 人员流动问题 . 5 7.6 二氧化碳分子结构的振动频率 . 5 7.7 二自由度机械振动 . 7 7.8 FIR 数字滤波器最优化设计12 . 8 7.9 弹性梁的柔度矩阵 . 10 7.10 用二次样条函数插值 5 个点 . 11 7.11 飞行器三维空间运动的矩阵描述 . 12 7.12 金融公司支付基金的流动 . 14 7.13 质谱图
2、实验结果分析 . 15 7.14 用特征方程解 Fibonacci 数列问题 . 16 7.15 简单线性规划问题. 18 第 2 章 时域中的离散信号和系统 1 1 第第 7 章章 在科技及工程中的应用实例在科技及工程中的应用实例 7.1 由拉压杆组成的桁架结构 由 13 根拉压杆件组成的桁架结构， 如图 7-1 所示，13 个平衡方程已给出，它们来自 6 个中 间节点，每个节点有 x,y 两个方向的平衡方程， 还有一个整体结构的 y 方向平衡方程。 现求其各 杆所受的力。 解解：按照题给方程组改写成矩阵形式，令 11 22 11 cos14/ 162 1420.6585 cos16/ 16
3、2 1620.7071 sin16/ 162 1420.7526 k k k =+= =+= =+= 列 方程时假设各杆的受力均为拉力，其相应的方程组及化为矩阵后的形式为： 22 1 2 26 3 41 52 1 2 133 5 71 84 3 89 101 56 93 5 2 117 2 1112 123 8 132 11 F +k F=0 k1 0 0 0 0 0 0 0 0 0 0 0 -F +F =0 0 - F =2000 F +k F -k F=0 k F+F +k F =-1000 F +k F -F =0 k F +F = -500 F -k F -F =0 F +k F =
4、4000 k F -F =0, k F +F =-500 F +k F = 2000 F +k F =0 21 23 1 3 1 3 2 2 3 2 1 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 -k 0 0 1 k0 0 0 0 0 0 0 0 k0 1 0 k0 0 0 0 0 0 0 0 00 0 -1 0 0 1 k0 0 0 0 0 00 0 0 0 0 0 k1 0 0 0 0 0 0 0 0 -k -1 0 0 0 1 0 0 0 0 0 0 0 k0 0 0 1 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 k0
5、 0 0 0 0 0 0 0 0 0 0 0 k1 0 0 0 0 0 0 0 0 k0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 k0 1 1 2 3 4 5 6 7 8 9 10 11 12 13 F 0 F 0 F 2000 F0 F-1000 F0 F-500(7.1.1) 0F 4000F 0F -500 F 2000 F 0 F = 将它看作 A*F=B，编成的程序为 pla701，核心语句为给 A,B 赋值，再求 F=AB，结果为： F= -7236; 5117; 2000; -6969; 2812; 5117; -4883; -3167; 1883; 6969;
6、-6906; 4383; 4883 其中负号表示杆受的是压力。 7.2 格型梯形 滤波器系统函数的 推导 使用计算机解题 后，用矩阵模型几乎是 最简便的数学方法了。 这将给后续课的建模 图 7-1 由拉压杆件组成的桁架结构 图 7-2 三阶格型梯形滤波器的结构图 第 2 章 时域中的离散信号和系统 2 2 和计算带来革命性的好处。例如要求出图 7-2 所示的滤波器的系统函数： 先列出方程，令 q=z-1，得到 x1= u k3x4; x2=x1; x3=k3x2+x4; x4=qx7; x5= x2 k2x8; x6=x5; x7=k2x6+x8; x8=qx11; x9= x6 k1x12;
7、 x10=x9; x11=k1x10+x12; x12=qx10; x13=y= C0 x12+ C1x11+ C2x7+ C3x3 这是一组含有 13 个变量的 13 个联立方程，用过去的手工方法一个一个消元，理论上是可 行的，但它运算极其繁琐，可以预期，95%以上的师生恐怕一个小时也解不出来，而且做对的 概率极低。 用矩阵的思路和方法来解就完全不同，它不是通过消元来减少变量，而是想办法补上所有 的零元素，把方程扩充为完整的矩阵形式： 134 21 3324 47 5228 65 72 1 2 3 4 6 5 6 7 8 9 8 811 961 12 109 111 1012 1210 13
8、 10 11 12 13 - - - xu k x xx xk xx xqx xxk x xx xk xx xqx xxk x xx xk xx xqx x x x x x x x x x x x x x x = = =+ = = = = + = = = =+ = = 3 3 0, 0, 0, - , 0, 0, 0, 0, 0, 0, 0, 0, 0 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 0, , 0, 1, 0, 0, 0, k k = 2 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, q, 0, 0, 0, 0, 0, 0 0
9、, 1, 0, 0, 0, 0, 0, -, 0, 0, 0, 0, 0 0, 0, 0, k 2 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, , 0, 1, 0, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, q, 0, k 1 0 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, - , 0 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0 0, 0, 0, 0, 0, 0, 0, 0, k 1 3210 0, , 0, 1, 0 0, 0, 0, 0, 0, 0, 0
10、, 0, 0, q, 0, 0, 0 0, 0, C , 0, 0, 0, C , 0, 0, 0, C , C , 0 k 1 2 3 4 5 6 7 8 9 10 11 12 13 1 0 0 0 0 0 0 0 0 0 0 0 0 x x x x x x xu x x x x x x + ,/X = QX+PUX U = inv(I-Q)*P 看似把模型搞复杂了，其实计算却非常容易。程序 pla703 先对 P，Q 矩阵赋值，键入 W = inv(I-Q)*P，马上就得出了系统函数。 编程时要注意，本例虽然是数值计算，但计算的内容中带有 z 变换算子 q=z-1，所以 P，Q 矩阵仍然必
11、须用符号属性，对 P，Q 赋值时第一个元素必须取含 q 的算式。熟练后不必列出 Q 和 P 的矩阵形式，可以按其下标规律直接进行元素赋值。 用以下参数：k0=1, k1=1/4, k2=1/2, k3=1/3, C0=-0.2, C1=0.8, C2=1.5, C3=1，编成了程序 pla703。运行此程序就得到： 32321 32321 (13)244942.930.8244942.930.8 (13) 81513248151324 xqqqzzz W uqqqzzz = + = + 用矩阵模型解信号流图的最大优点是一步到位，依靠计算机，既快速，又极易查错。 7.3 计算频谱用的 DFT 矩
12、阵 有限长序列x(n)(0nN-1)有N个样本值。 它的傅里叶变换( )X k在频率区间 （0 2）的N个等间隔分布的点k=2k/N(0kN-1)上也有 N 个样本值。这两组有限长的 序列之间可以用简单的关系联系起来： 第 2 章 时域中的离散信号和系统 3 3 1 ( )( )01 N kn N n kx n WkN = X 其中()exp 2/ N Wj N=是一个相角2/ N为的单位向量，也称为旋转因子。X(k)就 称为x(n)的离散傅里叶变换(DFT)，也就是 x(n)的频谱。用矩阵来表示,可写成: 所以信号频谱的计算，可以简单地用一个矩阵乘法来完成。信号是 N1 维数组，矩阵 W W
13、 称为 DFT 矩阵，它是 NN 维的复数阵。把矩阵乘法看作一个变换，我们就可以把频谱计 算看作信号从时域到频域的线性变换。 这个矩阵运算可用下列几条 MATLAB 语句来实现，程序名为 pla704。 xn=input(xn=(长度为N 的数组) ), N=length(xn); % 输入数据 n = 0:1:N-1; k = 0:1:N-1; % 设定n 和k 的行向量 WN = exp(-j*2*pi/N); % 设定Wn 因子 nk = n*k; % 产生一个含nk 值的N 乘N 维的整数矩阵 WNnk = WN . nk; % 求出W 矩阵 Xk = xn * WNnk % 求出离散
14、傅里叶级数系数 plot(k,Xk) % 画出幅频特性 前两条语句无须解释。第三条语句把n转为列向量n，再与行向量k做矩阵乘，它产 生的是一个整数矩阵： 这个整数方阵恰好是算式中 W W 矩阵的指数部分。 第四条语句就产生了 W W 矩阵， 而最后一条语句就完成了 DFS 的全部运算。由于实际上离散傅里叶级数并不太用到。它已 被下节将介绍的离散傅里叶变换取代，而且两者的计算程序是完全相同的。所以如果写成 MATLAB 子程序，它可以命名为 DFT。 在这 5 行 DFT 子程序前面，加上输入语句： xn=input(xn= ); N=length(xn); 在子程序后面，加上绘图语句： sub
15、plot(2,1,1),stem(xn,.); subplot(2,1,2),stem(abs(Xk),.) 就得到本题的程序 pla604。 运行 pla604，并按提示输入 xn，即可在子图 1 上得到 xn 序列的波形，并在子图 2 上得到此信号的幅 频特性。例如输入序列为 xn=ones(1,64),zeros(1,192)，则得到的时域 波形及其频谱图如图 7-3。 计算每一个 DFT 分量X(i)，需要做N次复数乘法和N-1 次复数加法。而要完成整个 DFT，求出X，则需要做NN次复数乘法和N(N-1)次复数加法。另外，它占的数据存储量 为NN。在频谱计算中，N 的典型值是 102
16、4。所以参与运算的数据达百万，需要的存储量 0100200 0 0.5 1 0100200 0 50 100 图 7-3 序列 x(n)及其频谱 第 2 章 时域中的离散信号和系统 4 4 超过 16M， 乘法运算的次数超过数百万。因此，在低档的微机上，当N较大时，MATLAB 会拒绝执行这个程序，并警告内存不足。实际上，MATLAB 内置的是加快计算并节省内存快 速傅里叶变换(FFT)程序。 从这里也可以看出，矩阵给我们提供了一个组织海量数据进行运算的最好方法，对上 百万数据进行赋值、运算和绘图，只用了几行程序。线性代数的重要性之所以在近 20 年来 可与微积分相妣美，这是一个重要原因。如果
17、学线性代数而不涉及工程计算实践，那就无 法真正体会线性代数的精髓所在。 7.4 显示器色彩制式转换问题 彩色显示器使用红(R)、绿(G)和蓝(B)光的叠成效应生成颜色。 显示器屏幕的内表面由微粒 象素组成, 每个微粒包括三个荧光点: 红、绿、蓝。 电子枪位于屏幕的后方, 向屏幕上每个点 发射电子束。计算机从图形应用程序或扫描仪发出数字信号到电子枪, 这些信号控制电子枪设 置的电压强度. 不同 RGB 的强度组合将产生不同的颜色. 电子枪由偏转电磁场帮助瞄准以确 保快速精确地屏幕刷新。 图 7-4 彩色显示器的工作原理 颜色模型规定一些属性或原色, 将颜色分解成不同属性的数字化组合. 这就是色彩
18、制式的 转换问题。 观察者在屏幕上实际看到的色彩要受色彩制式和屏幕上荧光点数量的影响。. 因此每家计 算机屏幕制造商都必须在(R, G, B)数据和国际通行的 CIE 色彩标准之间进行转换, CIE 标准使用 三原色, 分别称为 X, Y 和 Z。 针对短余辉荧光点的一类典型转换是 0.610.290.150 0.350.590.063 0.040.120.787 R G B = X Y Z A=. 计算机程序把用 CIE 数据(X, Y, Z)表示的色彩信息流发送到屏幕， 求屏幕上的电子枪将这些数据 转换成(R, G, B)数据的方程。现在要根据 CIE 数据(X, Y, Z)计算对应的(R
19、, G, B)数据, 就是等式 A = 中 A 和 已知, 求 。如果 A 是可逆矩阵, 则由 A = 可得 = A1 . 解：解：在 Matlab 命令窗口输入以下命令 A = 0.61,0.29,0.15;0.35,0.59,0.063;0.04,0.12,0.787; % B = inv(A), 程序执行的结果为： B = 2.2586 -1.0395 -0.3473 -1.3495 2.3441 0.0696 0.0910 -0.3046 1.2777 可见将 CIE 数据(X, Y, Z)转换成(R, G, B)数据的方程为 =B 在彩色电视技术中，还有许多地方会用到线性变换.比如民
20、用电视信号的发送并不直接采用 第 2 章 时域中的离散信号和系统 5 5 (R, G, B)数据， 而 是 使用向量(Y, I, Q)来描述每种颜色.其中 Y 称为亮度信号， I 和 Q 为色差信号， 这样的制式可以达到彩色和黑白的兼容.如果屏幕是黑白的, 则只用到了 Y，这比 CIE 数据能提 供更好的单色图象。YIQ 与“标准”RGB 色彩之间的对应如下 Y I Q = 0.2990.5870.114 0.5960.2750.321 0.2120.5280.311 R G B 它的逆变换矩阵留给读者自行计算。 7.5 人员流动问题 某试验性生产线每年一月份进行熟练工与非熟练工的人数统计,
21、然后将1/6熟练工支援其 他生产部门, 其缺额由招收新的非熟练工补齐。新、老非熟练工经过培训及实践至年终考核有 2/5成为熟练工，假设第一年一月份统计的熟练工和非熟练工各占一半，求以后每年一月份统 计的熟练工和非熟练工所占百分比。 设第 n 年一月份统计的熟练工和非熟练工所占百分比分别为 xn和 yn, 记成向量 n n x y = n Z. 因为第一年统计的熟练工和非熟练工各占一半, 所以 1 1 x y = 1 Z= 1/2 1/2 。为了求以后每年的熟 练工和非熟练工所占百分比, 先求 n+1 Z与 n Z的关系式。根据已知条件可得: xn+1 = (1 1 6 )xn + 2 5 (
22、1 6 xn + yn) = 9 10 xn + 2 5 yn, yn+1 = (1 2 5 )( 1 6 xn + yn) = 1 10 xn + 3 5 yn, 即 9/10 2/5 1/103/5 = nn+1n ZZ = AZ n+1 Z= A n Z= A2 n 1 Z= = An 1 Z. 将 A 对角化成为= -1 APP可得到简明易算的结果： n+1 Z= ( ) = 11 n -1n-1 PPPP ZZ。 键入命令P,D = eig(A)，得 0.97010.7071 0.2425 0.70 10 , 0 0.571 = PD，于是便有： 1 1 1 0.97010.7071
23、0.97010.70710.80.3/(2 ) 0.2425 0.70710.2425 100.5 0.50.7071 0.20.3/(02 )0.5 n n n n n n x y + + = + , 当 n 不断增加时， 1n x + ， 1n y + 分别趋向于 0.8 和 0.2。 7.6 二氧化碳分子结构的振动频率 二氧化碳分子可看成中间一个碳原子，左右分别以弹簧（化学键）联接两个氧原子，构成 一个三质量振动系统。按牛顿定律，其动力学方程为： 第 2 章 时域中的离散信号和系统 6 6 12112 23221123 33223 () ()()2 () o c o m xk xxkxk
24、x m xk xxk xxkxkxkx m xk xxkxkx = = + = =+ = = 其中： 222 312 123 222 , d xd xd x xxx dtdtdt = k=14.2e2 为化学键的弹性常数，mo=16*1.6605e-27, mc=12*1.6605e-27 分别为氧原子和碳原子的质量， x1,x2,x3为各原子的轴向位移。 今要求其分子振动的频率。 解：令 111122223333 ,yxyxyxyxyxyx = 设 h=k/mo=5.3448e+028，则 k/mc=h*mo/mc=h*4/3，将方程组写成矩阵形式，便可得到： 11211 1111 2213
25、22 2222 33323 3 2333 0000 100000 20 4 /30-8 /30 4 /3 001000 0000- 000010 o c o m ykxkxyyhh xyxx m ykxkxkxyyhhh xxxy hhyym ykxkx xxxy = + = = + = = = = d dt X = AX 这是标准的一阶线性微分方程组的矩阵形式。它在初始条件 X=X0下的解为： *= -1 00 XX exp(At)XP exp(t) P 其中 P， 分别为 A 的特征向量和特征值。因为题目只要求求振动频率而不是全部过程， 所以也不需要求特征向量，这里就采用分步的算法，先用
26、f=poly(A)求特征方程的系数多项式， 再用 L=roots(f) 求它的特征根。 于是可编制 MATLAB 程序 pla709，核心语句如下： k=14.2e2, mo=16*1.6605e-27, mc=12*1.6605e-27, h=k/mo A=0,-h,0,h,0,0;1,0,0,0,0,0;0,h*4/3,0,-h*8/3,0,h*4/3;0,0,1,0,0,0;0,0,0,h,0,-h;0,0,0,0,1,0 f=poly(A), L=roots(f) 程序运行的结果是： A = 1.0e+029 * 0 -0.5345 0 0.5345 0 0 0 0 0 0 0 0 0
27、 0.7126 0 -1.4253 0 0.7126 0 0 0.0000 0 0 0 0 0 0 0.5345 0 -0.5345 0 0 0 0 0.0000 0 f = 1.0e+071* 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 2.5711 L= 1.0e+014 * -0.0000 4.4269i；-0.0000 2.3119i；0.0000 0.0000i L 是人为缩减显示格式后的形式，f 是特征多项式的系数向量，由于各元素值差别太大而 无法用一个矩阵显示，只好单个显示，依次键入 f(1),f(7)，得知 f(1)=1, f(2)
28、=-0.0119, f(3)=2.4942e+029, f(4) = -6.0409e+027, f(5) = 1.0474e+058, f(6) = -2.8876e+056, f(7) = 2.5711e+071 L 是矩阵 A 的特征值，矩阵是六阶的，共有六个特征值，除去两个零后，剩下两对共轭虚 题 7-8 图 二氧化碳分子结构 第 2 章 时域中的离散信号和系统 7 7 数 1.0e+014 * 4.4269i，2.3119i，分别表示了两种振动角频率【弧度/秒 】。 将 L 除以 2 后可得到赫兹为单位的两种振动频率 70.456 和 36.795 THz 【太赫兹】 ， 用这个 频
29、率除以光速 3*108【米/秒】将得到两种波长 4.258 和 8.533m【微米】 。 7.7 二自由度机械振动 图 7-9 表示了一个由两个质量、两个弹簧和两个 阻尼器构成的二自由度振动系统， 今要求在给定初始 位置和初始速度下的运动。 解：解： 设 x1 和 x2 分别表示两个质量关于它们平衡 位置的偏移，则此二自由度振动系统的一般方程为： 1 11212212122 22221221 ()()0 ()()0 m xcc xc xkkxk x m xc xxkxx += += (7.10.1) 可写成矩阵形式： MX+CX+KX = 0 (7.10.2) 其中 11221221 1222
30、22 0 , 0 mccckkkx mcckkx + = M =CKX (7.10.4) 这是一个四阶的常系数齐次微分方程组，给定其初始位置 X0和初始速度 Xd0： 10101 0 000 20202 0 , d d d xxx xxx = XXX (7.10.4) 引进符号 xd是为了在 MATLAB 编程中给导数以变量名的需要。用解析法求这个方程的解 非常麻烦，只有当 C=0，即无阻尼时，系统可解耦为两种独立的振动模态，通常书上只给出解 耦简化后的解。 有了 MATLAB 工具，根本无需设 C0，也无需解耦，就可以用很简洁的程序求出其数值 解。其基本思路是把原始非常化成典型的四个一阶方程
31、构成的状态方程组： Y = AY (7.10.5) 这个方程在初始条件 Y=Y0下的解为 A t 0 Y = Y e。用 MATLAB 表示为 Y=Y0*expm(A*t)。 其中 expm 表示把(A*t)看成矩阵来求其指数。已经知道标量 x 求指数的方法，求矩阵指数虽然 要麻烦一些，但概念是相仿的。我们不必都弄清细节。MATLAB 提供了 expm, expm1, 等多种 函数供用户调用。所以我们只要把 Y,Y0 和 A 找到就行。 先把方程（2）写成两个一阶矩阵方程： d d ddd X = XXX0I = X-M K -M C XX = X = -M CX -M KX (7.10.6)
32、 于是 0 0 dd0 XX0I Y =, Y =, A = XX-M K-M C (7.10.7) 对于本题的二自由度系统， 11 22 d d xx xx = d XX和， 所 以 Y 和 Y0 都是 41 的单列变量； 图 7-8 两自由度机械振动模型 第 2 章 时域中的离散信号和系统 8 8 由于 A 中的四个元素都是 22 方阵，A 是 44 方阵。对于更多自由度的系统，公式(7)都是正 确的。 下面给出二自由度系统的一个数值例， 设 m1=1; m2=9; k1 = 4; k2=2; c1 和 c2可由用户输入。 求在初始条件 x0 = 1;0和 xd0 = 0;-1下，系统的输
33、出 x1,x2 曲线。 根据上面的模型可以写出程序pla710， 运行此程序， 输入c1=0.2,c2=0.5所得的结果见图7-9。 从中可清楚地看到振动的两种模态。特别是 x1 的运动反映了两种模态的迭合。给出不同的初始 条件，各模态的幅度也会变化。输入 c1=0,c2=0 所得的结果见图 7-10，这时两种振动模态是可 解耦的。可以看出，用计算机解题时，问题和数据的复杂性对解题的过程并无根本影响。通常 我们选择比较简单，且有解析解的数据作为检验程序之用。一旦确定程序正确，它就可用在很 复杂的情况。例如作者就曾把 b 本例题的核心语句用在一个五自由度的系统上，解一个十阶的 线性方程组，同样得
34、出满意的结果。 7.8 FIR 数字滤波器最优化设计12 设给定滤波器的预期幅频特性 D()在=i(i=1,2,K)的 K 个频点上的值，若实际滤 波器的脉冲响应的长度为 N=2L+1，则其幅特性与预期特性相拟合的方程为： 0 ()( )cos() L iiii n Ad nnDD = = (i=1, 2, , K) 这 K 个方程联立。由于() i A是 L+1 个谐波的线性组合，这方程组中的待定参数 d(n) 有 L+1 个。 如果 K=L+1，即方程数与未知数的数目相等，则这个线性方程组是适定的，恰好可以 解了这些系数。 如果 KL+1，即方程数比未知数的数目多，形成了所谓超定方程组。那
35、就不可能找到 精确满足这些方程组的系数 d(n)，只能找到最近似地满足这些方程的最小二乘解。 如果 KL+1，则形成了不定方程组，那会有无穷个解，工程上是没有意义的。所以只 考虑 KL+1 的情况。 e= =D Pd 图 7-10 二自由度振动波形 c1=c2=0 第 2 章 时域中的离散信号和系统 9 9 其中：e 为单列 K 元误差向量，D 为预期幅特性在样本点列上的 K 元单列向量，d 为 L+1 元待定系数单列向量，而 P 则是由 cos（ni）组成的 K(L+1)的系数矩阵。 11110 22221 1 cos() cos( 1 cos() cos( , 1 cos() cos( K
36、KKKL eDLd eDLd eDLd = eDPd ） ） ） 前面已经得到它的最小二乘解的公式： () -1 TT d = P PP D 用 MATLAB 计算此式特别方便，可以表示为 d= =inv（P* *P）* *P* *D，也可以更简便地 用左除运算 d=PD 来完成。 例例 7.8 举一个数字例，要求设计长度为 N=9( 有 5 个待定系数)的 FIR 滤波器，要使 它在 0 之间的八个频点上逼近预期的低通幅特性 D： =0, 0.33, 0.67, 1, 1.5, 2, 2.5, 3.14; D=1,1,1,1,0,0,0,0; 解：解：列出方程组的完整形式： 1111 222
37、2 3333 4444 5555 6666 7777 888 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 coscos2cos3cos4 1 2 3 4 5 6 7 88 (0) (1) (2)* (3) (4) D D d D d D d D d D d D D = P d = D 看来系数矩阵 P 的赋值比较麻烦，其实，回想求离散傅里叶变换时的做法，可以利用 列矩阵 w=w1; w
38、2; ; w8乘行矩阵0:4形成一个 85 矩阵，然后求余弦得到 P。这可以用 MATLAB 语句 P=cos(w*0:4)轻而易举地列出。因此用下列几行 MATLAB 程序 ag1010 就 可以方便地完成最小二乘最优滤波器的设计。 N=9; L=floor(N-1)/2); % 列出待求系数序数 w=0,0.33,0.67,1,1.5,2,2.5,3.14; % 列出频率向量 D=1,1,1,1,0,0,0,0; % 列出预期幅特性向量 P=cos(w*0:L); % 列出P 矩阵 d=inv(P*P)*P*D % 求出待求系数 程序运行的结果为 d = 0.3981 0.6039 0.2
39、137 -0.1205 -0.1506 知道 d(n)以后， 就可以计算滤波器的频 率特性进行校核。得到的幅频特性见图 7-11。因为 L=4，最高的谐波次数为 4，在 0 之间幅特性摆动最多两个周期，达到 峰值的次数应为四次，这从图上也看得很 清楚。 图 7-11 设计出的最优化滤波器的频率特性 第 2 章 时域中的离散信号和系统 10 10 7.9 弹性梁的柔度矩阵 设简支梁如图 7-12 所示，在梁的 三个位置分别施加力 f1， f2和 f3后， 在 该处产生的综合变形为图示的 y1，y2 和 y3，通常称为挠度。根据虎克定律， 在材料未失去弹性的范围内，力与它 引起的变形呈线性关系，可
40、以写出完 全的矩阵形式为： 11112131 22122232 33132333 ydddf ydddf ydddf = y = D*f 不难从这个矩阵乘法关系中得知矩阵 D 的物理意义。假如只施加一个力 f1，其余两个力 f2 和 f3为零，则引起的挠度为：y1=d11*f1，y2=d21*f1，y3=d31*f1，如果加的力为一个单位，则在 1， 2，3 三个位置引起的挠度分别为 d11，d21和 d31。可以用同样的方法来理解其他的 d，利用这个 概念， 可以用实测的方法来得到矩阵 D 中的各个元素。这些挠度元素愈大，表明这个梁愈柔软， 所以矩阵 D 被称作柔度矩阵。 柔度矩阵的逆就是刚
41、度矩阵 K，K= =D 1，在本例中它也是一个 33 矩阵。 11112131 22122232 33132333 fkkky fkkky fkkky = 其物理意义为造成单位挠度所需要的力。这些力愈大，表明这个梁愈刚劲。 在图 7-3 所示的情况下，柔度矩阵的 9 个元素中，没有一个为零。这说明它们之间相互耦 合。在位置 1 施加力不但会引起该处的变形，同时也必然引起位置 2 和 3 处的变形；在其他位 置加力亦然。如果矩阵变为对角矩阵，那也就意味着互相独立。位置1 加的力只引起位置1 处的挠 度，常识告诉我们，在图7-3 所示的简支梁上是不可能出现这种情况的。 现在举一个数字的例子，设柔度矩阵 .005 .002 .001 .002 .004 .003 .