1、第2章MATLAB的數(shù)值計(jì)算,所謂數(shù)值計(jì)算，就是指計(jì)算的表達(dá)式、變量中不得包含未經(jīng)定義的自由變量。 本章主要內(nèi)容： 21 變量及其賦值 22矩陣和數(shù)組的算術(shù)運(yùn)算 23線性方程組的求解 24矩陣函數(shù) 25多項(xiàng)式運(yùn)算 2. 6 練習(xí),2.1 變量及其賦值,在命令窗口可直接輸入變量并賦值，無(wú)須聲明其類型和維數(shù)，MATLAB將自動(dòng)處理。 MATLAB 的變量名稱由英文大小寫(xiě)52個(gè)字母、數(shù)字和下劃線等符號(hào)組成，并且第一個(gè)字符必須是英文字母，變量和常量最長(zhǎng)允許31個(gè)字符。須注意的是，在默認(rèn)狀態(tài)下，MATLAB對(duì)字母大小寫(xiě)敏感 。,MATLAB數(shù)據(jù)格式與顯示形式,MATLAB既可用傳統(tǒng)的十進(jìn)制數(shù)表達(dá)數(shù)值，也

2、可以用科學(xué)計(jì)數(shù)表達(dá)數(shù)值，用e代表10的指數(shù)形式，用i和j來(lái)代表虛數(shù)。 MATLAB內(nèi)部數(shù)據(jù)格式只有一種，是IEEE浮點(diǎn)標(biāo)準(zhǔn)的雙精度二進(jìn)制（64位），相應(yīng)于十進(jìn)制的16位有效數(shù)，范圍為10-308 10+308。 為了人機(jī)交互的友好性，數(shù)據(jù)輸出顯示格式有8種 ?？捎貌藛芜x項(xiàng)或format命令選擇。,MATLAB定義的特別變量及其意義,22矩陣和數(shù)組的算術(shù)運(yùn)算,2.2.1 矩陣的定義,MATLAB中的變量或常量都代表矩陣； 矩陣變量中的元素用“( )”中的數(shù)字來(lái)注明其下標(biāo)。維數(shù)不同，數(shù)字個(gè)數(shù)也不同，數(shù)字之間用“ ，”隔開(kāi) 使用下標(biāo)，用戶可以單獨(dú)給元素賦值，如x(2)=1.7321，a(2,3)=6

3、等,實(shí)例：訪問(wèn)矩陣中的元素,x=(0:pi/4:2*pi); y=sin(x); A=x,y %先構(gòu)造矩陣A 上面構(gòu)造的A由列向量x和在x基礎(chǔ)上生成的列向量y合并組成。由定義，x是一長(zhǎng)度為9的列向量，因此這樣構(gòu)造的A是一個(gè)9 2的矩陣。 b=A(2,1) 此命令訪問(wèn)了A的第二行第一列元素。,創(chuàng)建特殊矩陣的函數(shù),zeros() 矩陣元素都是零 ones() 矩陣元素都是1 注意： () 中的數(shù)字用于定義矩陣維數(shù)， ()內(nèi)不是數(shù)字的情況參看Help。 rand() 矩陣元素是0到1之間均勻分布的隨機(jī)數(shù) randn() 矩陣元素是正態(tài)分布的隨機(jī)數(shù),2.2.2 矩陣運(yùn)算,MATLAB提供矩陣算術(shù)運(yùn)算有

4、： 加 (+)、減 (-)、乘 (*)、除 (/)、冪次方 ()、轉(zhuǎn)置() 注意，對(duì)于矩陣來(lái)說(shuō)，乘法有左乘右乘之分，A *B時(shí)要求A的列數(shù)等于B的行數(shù)，即滿足矩陣相乘的條件。 MATLAB中，為求解線性方程組，引入了矩陣除法。 MATLAB中的矩陣除法有左除()和右除(/)兩種。 以下是矩陣運(yùn)算的合法命令形式： A+B 、A-B、A*B、A/B、AB、Ap（此處p只能是標(biāo)量，A必須為方陣）、A,矩陣的左除和右除,Case 1：方程為 A *x = B 若A非奇異，用左除解得 x = A B 如A為方陣，則A 形式上相當(dāng)于A的逆陣。但須注意，MATLAB中，求解線性方程組并不采用求逆的方法。 C

5、ase 2：方程為 x *A = B 若A非奇異，用右除解得 x = B / A 如A為方陣， / A形式上相當(dāng)于A的逆矩陣。,2.2.3 數(shù)組運(yùn)算,和矩陣運(yùn)算符比較，除了加減符號(hào)外，其余的數(shù)組運(yùn)算符號(hào)均須加上一個(gè) “ . ”符號(hào)，即： 加 (+)、減 (-)、乘 (.*)、左除 (./)、右除(.)、冪次方 (.) 數(shù)組運(yùn)算是矩陣中元素對(duì)元素的操作，因而進(jìn)行數(shù)組運(yùn)算時(shí)，加、減、乘、除均要求A、B有相同的維數(shù)。對(duì)加減運(yùn)算而言，A和B中有一個(gè)是標(biāo)量也可以。 數(shù)組相除也有左除和右除之分，其中A./B表示A的元素為被除數(shù)，而A.B 則表示A的元素是除數(shù)。 數(shù)組運(yùn)算A .p中，如p是標(biāo)量，則此運(yùn)算是指

6、對(duì)A中每個(gè)元素按冪次p進(jìn)行乘方運(yùn)算。但若數(shù)組運(yùn)算A.B中，A為標(biāo)量，B為矩陣，則A.B的含義為A B(i)。,范數(shù),范數(shù)：信號(hào)處理問(wèn)題中，信號(hào)常表示為線性空間中的一個(gè)向量，而向量的長(zhǎng)度是一個(gè)非常重要的量。利用上述的數(shù)組運(yùn)算A.p，向量X的長(zhǎng)度的平方很容易表示為 sum(X.2) 范數(shù)（norm）概念可容易地用來(lái)度量向量長(zhǎng)度。 MATLAB 中設(shè)置了專用的norm函數(shù) norm(X，p)= sum(X.p)(1/p) 注意，這里用到了上面的sum()函數(shù)。 常用的p值為1，2，。p=2時(shí)返回X的二階范，即通常意義下向量的歐氏長(zhǎng)度。,2.3 線性方程組的求解,實(shí)際情況下，通常遇到的線性方程組的形式

7、為 A * X = B 其中，系數(shù)矩陣A不必為方陣 設(shè)其為mn維矩陣 。 (1)mn,系統(tǒng)為超定，X=AB，給出系統(tǒng)最小二乘解； (2)m=n,如矩陣A非奇異，X=AB，給出系統(tǒng)的精確解； (3)mn,系統(tǒng)為欠定，X=AB，給出系統(tǒng)的一組通解，不具惟一性。,例：均衡濾波器,圖中是一個(gè)5抽頭的橫向?yàn)V波器，用作均衡濾波器， C1表示抽頭因子構(gòu)成的列向量,其輸入為h(n) 。設(shè)h(n)長(zhǎng)為5點(diǎn) h=-0.1,0.2,0.9,-0.3,0.15 理想情況下，均衡濾波器輸出y(n)的理想值應(yīng)為 y = 0 0 0 0 1 0 0 0 0 長(zhǎng)度為9點(diǎn)。我們的任務(wù)是確定抽頭因子的數(shù)值。,例：均衡濾波器的最小

8、均方誤差(MMSE)解與迫零(Zero forcing)解,MMSE解： 1. 求出h的卷積矩陣A = convmtx(h ,5) ； 2. 令y1=y，左除求解超定方程組 A * (C1) = (y1) 即得MMSE解，即最小二乘解C1。 Zero forcing解： 1. 去掉A中的上下各兩行，構(gòu)成方陣B = A(3:7;1:5) 2. 令y2=0 0 1 0 0，左除求解 B * (C2) = (y2) 即得Zero forcing解C2。,運(yùn)用范數(shù)概念比較兩個(gè)解,對(duì)于得到的C1，C2，均衡濾波器輸出分別為 y11 = conv(h, C1 ); y22 = conv(h,C2 ) 長(zhǎng)度

9、相同，均為9點(diǎn)。 y11，y22與理想輸出y之間的誤差分別為y-y11和y-y22，這兩個(gè)誤差向量的二階范數(shù)就是方均根誤差。 也可用均方誤差即二階范的平方來(lái)衡量。 Minimum MSE解優(yōu)于Zero forcing解（Why？）,2.4 矩陣函數(shù),本節(jié)給出了一些重要的矩陣函數(shù)，在信號(hào)處理中，這些函數(shù)常會(huì)用到。 對(duì)初學(xué)者而言，本節(jié)內(nèi)容如感困難，可略過(guò)不學(xué)，不影響后續(xù)內(nèi)容的學(xué)習(xí)。 矩陣范數(shù)、秩與條件數(shù)。 行列式與逆矩陣。 三種矩陣分解。 特征值與特征向量。 奇值分解（SVD）。,2.5 多項(xiàng)式運(yùn)算,多項(xiàng)式在工程中有著廣泛的應(yīng)用，如實(shí)驗(yàn)數(shù)據(jù)的多項(xiàng)式擬合表達(dá)等。,2.5.1 多項(xiàng)式的定義,1直接定義

10、 令a(x) 代表一個(gè)如下的多項(xiàng)式 在MATLAB中，這個(gè)多項(xiàng)式可以用它的系數(shù)向量A=a(1)，a(2)，a(n)，a(n+1)來(lái)表示。 2 用命令poly創(chuàng)建 如果A是矩陣，則poly(A)將創(chuàng)建A矩陣的特征多項(xiàng)式；如果A是向量，則poly(A)將創(chuàng)建以A中各元素為根的多項(xiàng)式。,2.5.2 多項(xiàng)式的四則運(yùn)算,1.多項(xiàng)式相乘,以兩個(gè)多項(xiàng)式為例： conv函數(shù)本用于兩個(gè)有限長(zhǎng)序列的卷積，但卷積和的運(yùn)算也符合多項(xiàng)式相乘規(guī)則。因此在命令窗口鍵入 a = 2，4，6，8，b = 3，6，9， c = conv(a,b) 后，可得 c = 6 24 60 96 102 72 而,2.多項(xiàng)式相加,MATL

11、AB規(guī)定，只有長(zhǎng)度相同即維數(shù)相同的向量才能相加。因此，必須把短的向量前面補(bǔ)以若干個(gè)零元素，也即對(duì)多項(xiàng)式高次項(xiàng)系數(shù)以零表示，才能用MATLAB的矩陣加法運(yùn)算符。 可以編寫(xiě)一個(gè)子程序polyadd.m，以使計(jì)算機(jī)自動(dòng)完成兩個(gè)多項(xiàng)式符合相加的條件。關(guān)于子程序的編寫(xiě)，下一章將會(huì)學(xué)習(xí)到。 function y=polyadd(x1,x2) n1=length(x1);n2=length(x2); if n1n2 x2=zeros(1,n1-n2),x2; elseif n1n2 x1=zeros(1,n2-n1),x1; end, y=x1+x2;,3.多項(xiàng)式相除,相除是相乘的逆運(yùn)算，不能除盡時(shí)會(huì)有余子

12、式產(chǎn)生。 故在MATLAB中，多項(xiàng)式相除后會(huì)返回一項(xiàng)余子式。能夠除盡時(shí)，余子式全為零。 函數(shù)形式為： q,r=deconv(c,a) 其中q是商，r是余子式。 注意deconv的原意是反卷積，工程中也稱之為逆濾波。但由于deconv允許出現(xiàn)余子式，這個(gè)函數(shù)不能保證得到滿意的逆濾波解。,2.5.3 多項(xiàng)式求導(dǎo)、求根和求值,1.多項(xiàng)式求導(dǎo)數(shù)（polyder） 函數(shù)形式： e=polyder(c) 對(duì)上面的c(x),使用這一函數(shù)得到 e = 30 96 180 192 102 表示,2.多項(xiàng)式求根,函數(shù)形式： ra = roots (a) 其中a表示多項(xiàng)式a(x)。 注意，前面已經(jīng)學(xué)習(xí)過(guò)根據(jù)多項(xiàng)式的

13、根創(chuàng)建多項(xiàng)式系數(shù)的函數(shù)poly，實(shí)際上，這是roots的逆運(yùn)算。即 a = poly (ra) 這里ra是一個(gè)由根構(gòu)成的向量。,3.多項(xiàng)式求值（polyval）,將多項(xiàng)式a中的自變量x賦予值xv時(shí)，該多項(xiàng)式的值可用F=polyval(a,xv)求得。 其中xv可以是復(fù)數(shù)，也可以是矩陣或數(shù)組，對(duì)于后者，polyval對(duì)輸入變?cè)倪\(yùn)算結(jié)果也是矩陣或數(shù)組。 這一性質(zhì)可以用于求線性系統(tǒng)的頻率特性。,例2-1： LTI系統(tǒng)頻率響應(yīng),求此系統(tǒng)的頻率響應(yīng)并畫(huà)出幅頻與相頻特性。 解： 設(shè)a為該系統(tǒng)分母系數(shù)向量，b為系統(tǒng)分子系數(shù)向量，即： a=1 2 3 4， b=1 2 3 該系統(tǒng)的頻率特性為 程序中使用po

14、lyval函數(shù)求出了系統(tǒng)幅頻特性和相頻特性在頻率數(shù)組w各個(gè)取值處的數(shù)值。 MATLAB程序見(jiàn)教材p3233。頻率特性圖見(jiàn)圖2-2。,例2-1 設(shè)有LTI系統(tǒng)的系統(tǒng)函數(shù)為：,2.5.4 多項(xiàng)式擬合,函數(shù)：p=polyfit(x,y,n) 用于在最小均方誤差準(zhǔn)則下對(duì)數(shù)據(jù)進(jìn)行n階多項(xiàng)式曲線擬合。 其中x，y是已知的N個(gè)數(shù)據(jù)點(diǎn)坐標(biāo)向量，其長(zhǎng)度均為N，現(xiàn)欲用一個(gè)n階多項(xiàng)式來(lái)表示x，y之間的函數(shù)關(guān)系，待定的是該n階多項(xiàng)式的系數(shù)。 多項(xiàng)式階數(shù)n越大，擬合程度越好，設(shè)p為用于擬合的多項(xiàng)式系數(shù)個(gè)數(shù)，p = n + 1 。 例2-2：測(cè)量放大器3dB通頻帶的測(cè)定。參見(jiàn)pp3335。,2. 6 練習(xí),命令窗口中，鍵入：demo，出現(xiàn)MATLAB的demo演示窗口，然后選擇： 1 MATLAB下的M