MATLAB 曲線(xiàn)擬合與插值

1、第11章 曲線(xiàn)擬合與插值在大量的應(yīng)用領(lǐng)域中，人們經(jīng)常面臨用一個(gè)解析函數(shù)描述數(shù)據(jù)(通常是測(cè)量值 的任務(wù)。對(duì)這個(gè)問(wèn)題有兩種方法。在插值法里，數(shù)據(jù)假定是正確的，要求以某種方法描述數(shù)據(jù)點(diǎn)之間所發(fā)生的情況。這種方法在下一節(jié)討論。這里討論的方法是曲線(xiàn)擬合或回歸。人們?cè)O(shè)法找出某條光滑曲線(xiàn)，它最佳地?cái)M合數(shù)據(jù)，但不必要經(jīng)過(guò)任何數(shù)據(jù)點(diǎn)。圖11.1說(shuō)明了這兩種方法。標(biāo)有'o' 的是數(shù)據(jù)點(diǎn)；連接數(shù)據(jù)點(diǎn)的實(shí)線(xiàn)描繪了線(xiàn)性?xún)?nèi)插，虛線(xiàn)是數(shù)據(jù)的最佳擬合。11.1 曲線(xiàn)擬合曲線(xiàn)擬合涉及回答兩個(gè)基本問(wèn)題：最佳擬合意味著什么？應(yīng)該用什么樣的曲線(xiàn)？可用許多不同的方法定義最佳擬合，并存在無(wú)窮數(shù)目的曲線(xiàn)。所以，從這里開(kāi)始，

2、我們走向何方？正如它證實(shí)的那樣，當(dāng)最佳擬合被解釋為在數(shù)據(jù)點(diǎn)的最小誤差平方和，且所用的曲線(xiàn)限定為多項(xiàng)式時(shí)，那么曲線(xiàn)擬合是相當(dāng)簡(jiǎn)捷的。數(shù)學(xué)上，稱(chēng)為多項(xiàng)式的最小二乘曲線(xiàn)擬合。如果這種描述使你混淆，再研究圖11.1。虛線(xiàn)和標(biāo)志的數(shù)據(jù)點(diǎn)之間的垂直距離是在該點(diǎn)的誤差。對(duì)各數(shù)據(jù)點(diǎn)距離求平方，并把平方距離全加起來(lái)，就是誤差平方和。這條虛線(xiàn)是使誤差平方和盡可能小的曲線(xiàn)，即是最佳擬合。最小二乘這個(gè)術(shù)語(yǔ)僅僅是使誤差平方和最小的省略說(shuō)法。 圖11.1 2階曲線(xiàn)擬合在MATLAB 中，函數(shù)polyfit 求解最小二乘曲線(xiàn)擬合問(wèn)題。為了闡述這個(gè)函數(shù)的用法，讓我們以上面圖11.1中的數(shù)據(jù)開(kāi)始。» x=0 .1 .

3、2 .3 .4 .5 .6 .7 .8 .9 1;» y=-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;為了用polyfit ，我們必須給函數(shù)賦予上面的數(shù)據(jù)和我們希望最佳擬合數(shù)據(jù)的多項(xiàng)式的階次或度。如果我們選擇n=1作為階次，得到最簡(jiǎn)單的線(xiàn)性近似。通常稱(chēng)為線(xiàn)性回歸。相反，如果我們選擇n=2作為階次，得到一個(gè)2階多項(xiàng)式。現(xiàn)在，我們選擇一個(gè)2階多項(xiàng)式。» n=2; % polynomial order» p=polyfit(x, y, np =polyfit 的輸出是一個(gè)多項(xiàng)式系數(shù)的行向量。其解是y

4、= 9.8108x 2 20.1293x 0.0317。為了將曲線(xiàn)擬合解與數(shù)據(jù)點(diǎn)比較，讓我們把二者都繪成圖。» xi=linspace(0, 1, 100; % x-axis data for plotting» z=polyval(p, xi;為了計(jì)算在xi 數(shù)據(jù)點(diǎn)的多項(xiàng)式值，調(diào)用MATLAB 的函數(shù)polyval 。» plot(x, y, ' o ' , x, y, xi, z, ' : ' 畫(huà)出了原始數(shù)據(jù)x 和y ，用'o' 標(biāo)出該數(shù)據(jù)點(diǎn)，在數(shù)據(jù)點(diǎn)之間，再用直線(xiàn)重畫(huà)原始數(shù)據(jù)，并用點(diǎn)' : '線(xiàn)

5、，畫(huà)出多項(xiàng)式數(shù)據(jù)xi 和z 。» xlabel(' x ', ylabel(' y=f(x ', title(' Second Order Curve Fitting '將圖作標(biāo)志。這些步驟的結(jié)果表示于前面的圖11.1中。多項(xiàng)式階次的選擇是有點(diǎn)任意的。兩點(diǎn)決定一直線(xiàn)或一階多項(xiàng)式。三點(diǎn)決定一個(gè)平方或2階多項(xiàng)式。按此進(jìn)行，n+1數(shù)據(jù)點(diǎn)唯一地確定n 階多項(xiàng)式。于是，在上面的情況下，有11個(gè)數(shù)據(jù)點(diǎn)，我們可選一個(gè)高達(dá)10階的多項(xiàng)式。然而，高階多項(xiàng)式給出很差的數(shù)值特性，人們不應(yīng)選擇比所需的階次高的多項(xiàng)式。此外，隨著多項(xiàng)式階次的提高，近似變得不夠光滑

6、，因?yàn)檩^高階次多項(xiàng)式在變零前，可多次求導(dǎo)。例如，選一個(gè)10階多項(xiàng)式» pp=polyfit(x, y, 10 ;» format short e % change display format» pp.' % display polynomial coefficients as a columnans =-4.6436e+0052.2965e+006-4.8773e+0065.8233e+006-4.2948e+0062.0211e+006-6.0322e+0051.0896e+005-1.0626e+0044.3599e+002要注意在現(xiàn)在情況下，多項(xiàng)式

7、系數(shù)的規(guī)模與前面的2階擬合的比較。還要注意在最小(-4.4700e-001和最大(5.8233e+006系數(shù)之間有7個(gè)數(shù)量級(jí)的幅度差。將這個(gè)解作圖，并把此圖與原始數(shù)據(jù)及2階曲線(xiàn)擬合相比較，結(jié)果如何呢？» zz=polyval(pp, xi; % evaluate 10th order polynomial» plot(x, y, ' o ' , xi, z, ' : ' , xi, zz % plot data» xlabel(' x ', ylabel(' y=f(x ', title('

8、2nd and 10th Order curve Fitting '在下面的圖11.2中，原始數(shù)據(jù)標(biāo)以'o' ，2階曲線(xiàn)擬合是虛線(xiàn)，10階擬合是實(shí)線(xiàn)。注意，在10階擬合中，在左邊和右邊的極值處，數(shù)據(jù)點(diǎn)之間出現(xiàn)大的紋波。當(dāng)企圖進(jìn)行高階曲線(xiàn)擬合時(shí)，這種紋波現(xiàn)象經(jīng)常發(fā)生。根據(jù)圖11.2，顯然， 越多就越好 的觀念在這里不適用。 圖11.2 2階和10階曲線(xiàn)擬合11.2 一維插值正如在前一節(jié)對(duì)曲線(xiàn)擬合所描述的那樣，插值定義為對(duì)數(shù)據(jù)點(diǎn)之間函數(shù)的估值方法，這些數(shù)據(jù)點(diǎn)是由某些集合給定。當(dāng)人們不能很快地求出所需中間點(diǎn)的函數(shù)值時(shí)，插值是一個(gè)有價(jià)值的工具。例如，當(dāng)數(shù)據(jù)點(diǎn)是某些實(shí)驗(yàn)測(cè)量的結(jié)果

9、或是過(guò)長(zhǎng)的計(jì)算過(guò)程時(shí)，就有這種情況?；蛟S最簡(jiǎn)單插值的例子是MATLAB 的作圖。按缺省，MATLAB 用直線(xiàn)連接所用的數(shù)據(jù)點(diǎn)以作圖。這個(gè)線(xiàn)性插值猜測(cè)中間值落在數(shù)據(jù)點(diǎn)之間的直線(xiàn)上。當(dāng)然，當(dāng)數(shù)據(jù)點(diǎn)個(gè)數(shù)的增加和它們之間距離的減小時(shí)，線(xiàn)性插值就更精確。例如，» x1=linspace(0, 2*pi, 60;» x2=linspace(0, 2*pi, 6;» plot(x1, sin(x1, x2, sin(x2, ' - '» xlabel(' x ', ylabel(' sin(x ', title('

10、; Linear Interpolation ' 圖11.3 線(xiàn)性插值圖11.3是sine 函數(shù)的兩個(gè)圖，一個(gè)在數(shù)據(jù)點(diǎn)之間用60個(gè)點(diǎn)，它比另一個(gè)只用6個(gè)點(diǎn)更光滑和更精確。如曲線(xiàn)擬合一樣，插值要作決策。根據(jù)所作的假設(shè)，有多種插值。而且，可以在一維以上空間中進(jìn)行插值。即如果有反映兩個(gè)變量函數(shù)的插值，z=f(x, y，那么就可在x 之間和在y 之間，找出z 的中間值進(jìn)行插值。MATLAB 在一維函數(shù)interp1和在二維函數(shù)interp2中，提供了許多的插值選擇。其中的每個(gè)函數(shù)將在下面闡述。為了說(shuō)明一維插值，考慮下列問(wèn)題，12小時(shí)內(nèi)，一小時(shí)測(cè)量一次室外溫度。數(shù)據(jù)存儲(chǔ)在兩個(gè)MATLAB 變量中

11、。» hours=1:12; % index for hour data was recorded» temps=5 8 9 15 25 29 31 30 22 25 27 24; % recorded temperatures» plot(hours, temps, hours, temps,' + ' % view temperatures» title(' Temperature '» xlabel(' Hour ', ylabel(' Degrees Celsius ' 圖

12、11.4 在線(xiàn)性插值下室外溫度曲線(xiàn)正如圖11.4看到的，MATLAB 畫(huà)出了數(shù)據(jù)點(diǎn)線(xiàn)性插值的直線(xiàn)。為了計(jì)算在任意給定時(shí)間的溫度，人們可試著對(duì)可視的圖作解釋。另外一種方法，可用函數(shù)interp1。t =t =22» t=interp1(hours, temps, 3.2 6.5 7.1 11.7 % find temp at many points!t =interp1的缺省用法是由interp1(x, y, xo來(lái)描述，這里x 是獨(dú)立變量(橫坐標(biāo) ，y 是應(yīng)變量(縱坐標(biāo) ，xo 是進(jìn)行插值的一個(gè)數(shù)值數(shù)組。另外，該缺省的使用假定為線(xiàn)性插值。若不采用直線(xiàn)連接數(shù)據(jù)點(diǎn)，我們可采用某些更光滑的

13、曲線(xiàn)來(lái)擬合數(shù)據(jù)點(diǎn)。最常用的方法是用一個(gè)3階多項(xiàng)式，即3次多項(xiàng)式，來(lái)對(duì)相繼數(shù)據(jù)點(diǎn)之間的各段建模，每個(gè)3次多項(xiàng)式的頭兩個(gè)導(dǎo)數(shù)與該數(shù)據(jù)點(diǎn)相一致。這種類(lèi)型的插值被稱(chēng)為3次樣條或簡(jiǎn)稱(chēng)為樣條。函數(shù)interp1也能執(zhí)行3次樣條插值。t =t =» t=interp1(hours, temps, 3.2 6.5 7.1 11.7, ' spline 't =注意，樣條插值得到的結(jié)果，與上面所示的線(xiàn)性插值的結(jié)果不同。因?yàn)椴逯凳且粋€(gè)估計(jì)或猜測(cè)的過(guò)程，其意義在于，應(yīng)用不同的估計(jì)規(guī)則導(dǎo)致不同的結(jié)果。一個(gè)最常用的樣條插值是對(duì)數(shù)據(jù)平滑。也就是，給定一組數(shù)據(jù)，使用樣條插值在更細(xì)的間隔求值。例如，

14、» h=1:0.1:12; % estimate temperature every 1/10 hour» t=interp1(hours, temps, h, ' spline ' ;» plot(hours, temps, ' - ' , hours, temps, ' + ' , h, t % plot comparative results» title(' Springfield Temperature '» xlabel(' Hour ', ylabel

15、(' Degrees Celsius '在圖11.5中，虛線(xiàn)是線(xiàn)性插值，實(shí)線(xiàn)是平滑的樣條插值，標(biāo)有' + '的是原始數(shù)據(jù)。如要求在時(shí)間軸上有更細(xì)的分辨率，并使用樣條插值，我們有一個(gè)更平滑、但不一定更精確地對(duì)溫度的估計(jì)。尤其應(yīng)注意，在數(shù)據(jù)點(diǎn)，樣條解的斜率不突然改變。作為這個(gè)平滑插值的回報(bào)，3次樣條插值要求更大量的計(jì)算，因?yàn)楸仨氄业?次多項(xiàng)式以描述給定數(shù)據(jù)之間的特征。關(guān)于樣條的更詳細(xì)信息可見(jiàn)下一章。 圖11.5 在不同插值下室外溫度曲線(xiàn)在討論二維插值之前，了解interp1所強(qiáng)制的二個(gè)強(qiáng)約束是很重要的。首先，人們不能要求有獨(dú)立變量范圍以外的結(jié)果，例如，interp1

16、(hours, temps, 13.5 導(dǎo)致一個(gè)錯(cuò)誤，因?yàn)閔ours 在1到12之間變化。其次，獨(dú)立變量必須是單調(diào)的。即獨(dú)立變量在值上必須總是增加的或總是減小的。在我們的例子里，hours 是單調(diào)的。然而，如果我們已經(jīng)定義獨(dú)立變量為一天的實(shí)際時(shí)間，» time_of_day=7:12 1:6 % start at 7AM,end at 6PMtime_of_day =7 8 9 10 11 12 1 2 3 4 5 6則獨(dú)立變量將不是單調(diào)的，因?yàn)閠ime_of_day增加到12，然后跌到1，再然后增加。如果用time_of_day代替interp1中的hours ，將會(huì)返回一個(gè)錯(cuò)誤。

17、同樣的理由，人們不能對(duì)temps 插值來(lái)找出產(chǎn)生某溫度的時(shí)間(小時(shí) ，因?yàn)閠emps 不是單調(diào)的。11.3 二維插值二維插值是基于與一維插值同樣的基本思想。然而，正如名字所隱含的，二維插值是對(duì)兩變量的函數(shù)z=f(x, y 進(jìn)行插值。為了說(shuō)明這個(gè)附加的維數(shù)，考慮一個(gè)問(wèn)題。設(shè)人們對(duì)平板上的溫度分布估計(jì)感興趣，給定的溫度值取自平板表面均勻分布的格柵。采集了下列的數(shù)據(jù)：» width=1:5; % index for width of plate (i.e.,the x-dimension» depth=1:3; % index for depth of plate (i,e,th

18、e y-dimension» temps=82 81 80 82 84; 79 63 61 65 81; 84 84 82 85 86 % temperature datatemps =82 81 80 82 8479 63 61 65 8184 84 82 85 86如同在標(biāo)引點(diǎn)上測(cè)量一樣，矩陣temps 表示整個(gè)平板的溫度分布。temps 的列與下標(biāo)depth 或y-維相聯(lián)系，行與下標(biāo)width 或x-維相聯(lián)系(見(jiàn)圖11.6 。為了估計(jì)在中間點(diǎn)的溫度，我們必須對(duì)它們進(jìn)行辨識(shí)。» wi=1:0.2:5; % estimate across width of plate&#

19、187; d=2; % at a depth of 2» zlinear=interp2(width, depth, temps, wi, d ; % linear interpolation» zcubic=interp2(width, depth, temps, wi,d, ' cubic ' ; % cubic interpolation» plot(wi, zlinear, ' - ' , wi, zcubic % plot results» xlabel(' Width of Plate ', y

20、label(' Degrees Celsius '» title( ' Temperature at Depth = ' num2str(d 另一種方法，我們可以在兩個(gè)方向插值。先在三維坐標(biāo)畫(huà)出原始數(shù)據(jù)，看一下該數(shù)據(jù)的粗糙程度(見(jiàn)圖11.7 。» mesh(width, depth, temps % use mesh plot» xlabel(' Width of Plate ', ylabel(' Depth of Plate '» zlabel(' Degrees Celsius

21、 ', axis(' ij ', grid 圖11.6 在深度d=2處的平板溫度 圖11.7 平板溫度然后在兩個(gè)方向上插值，以平滑數(shù)據(jù)。» di=1:0.2:3; % choose higher resolution for depth» wi=1:0.2:5; % choose higher resolution for width» zcubic=interp2(width, depth, temps, wi, di, ' cubic ' ; % cubic» mesh(wi, di, zcubic»

22、 xlabel(' Width of Plate ', ylabel(' Depth of Plate '» zlabel(' Degrees Celsius ', axis(' ij ', grid上面的例子清楚地證明了，二維插值更為復(fù)雜，只是因?yàn)橛懈嗟牧恳３指?。interp2的基本形式是interp2(x, y, z, xi, yi, method。這里x 和y 是兩個(gè)獨(dú)立變量，z 是一個(gè)應(yīng)變量矩陣。x 和y 對(duì)z 的關(guān)系是z(i, : = f(x, y(i 和 z(:, j = f(x(j, y.也就是，當(dāng)x

23、 變化時(shí)，z 的第i 行與y 的第i 個(gè)元素y(i相關(guān)，當(dāng)y 變化時(shí)，z 的第j 列與x 的第j 個(gè)元素x(j相關(guān)，。xi 是沿x-軸插值的一個(gè)數(shù)值數(shù)組；yi 是沿y-軸插值的一個(gè)數(shù)值數(shù)組。 圖11.8 二維插值后的平板溫度可選的參數(shù)method 可以是 'linear' ，'cubic' 或'nearest' 。在這種情況下，cubic 不意味著3次樣條，而是使用3次多項(xiàng)式的另一種算法。linear 方法是線(xiàn)性插值，僅用作連接圖上數(shù)據(jù)點(diǎn)。nearest 方法只選擇最接近各估計(jì)點(diǎn)的粗略數(shù)據(jù)點(diǎn)。在所有的情況下，假定獨(dú)立變量x 和y 是線(xiàn)性間隔和單調(diào)

24、的。關(guān)于這些方法的更多的信息，可請(qǐng)求在線(xiàn)幫助，例如，» helpinterp2，或參閱MATLAB 參考手冊(cè)。11.4 M 文件舉例雖然對(duì)于許多應(yīng)用，函數(shù)interp1和interp2是很有用的，但它們限制為對(duì)單調(diào)向量進(jìn)行插值。在某些情況，這個(gè)限制太嚴(yán)格。例如，考慮下面的插值：» x=linspace(0, 5;» y=1-exp(-x.*sin(2*pi*x;» plot(x, y 圖11.9 函數(shù)1-exp(-x.*sin(2*pi*x的曲線(xiàn)函數(shù)interp1可用來(lái)在任何值或x 的值上估計(jì)y 值。yi =然而，interp1不能找出對(duì)應(yīng)于某些y 值的

25、x 值。例如，如在圖11.9上所示，考慮尋找y=1.1處的x 值： 圖11.10 給y 值在函數(shù)曲線(xiàn)上求x 的值» plot(x, y, 0, 5, 1.1 1.1 從圖11.10上，我們看到有四個(gè)交點(diǎn)。使用interp1，我們得到：? Error using => table1First column of the table must be monotonic.這個(gè)函數(shù)interp1失敗，由于y 不是單調(diào)的。在本章精通MATLAB 工具箱所說(shuō)明的M 文件例子，消除了單調(diào)性的要求。» table=x; y.' ; % create column orient

26、ed table from dataxi =這里使用了線(xiàn)性插值，函數(shù)mminterp 估計(jì)了y=1.1處的四個(gè)點(diǎn)。由于函數(shù)mminterp 的一般性質(zhì)，要插值的數(shù)據(jù)是由面向列矩陣給出，在上面的例子中稱(chēng)作為表(table 。第二個(gè)輸入?yún)⒘渴潜凰阉骶仃噒able 的列，第三個(gè)參量是要找的值。這個(gè)精通MATLAB 工具箱函數(shù)的主體由下面給出：function y=mminterp(tab, col, val% MMINTERP 1-D Table Search by Linear Interpolation.% Y=MMINTERP(TAB,COL,VAL linearly interpolates

27、 the table% TAB searching for the scalar value VAL in the column COL.% All crossings are found and TAB(:,COL need not be monotonic.% Each crossing is returned as a separate row in Y and Y has as% many columns as TAB.Naturally,the column COL of Y contains% the value VAL. If VAL is not found in the ta

28、ble,Y=.% Copyright (c 1996 by Prentice-Hall,Inc.rt, ct=size(tab;if length(val > 1, error(' VAL must be a scalar. ', endif col>ct|col < 1, error(' Chosen column outside table width. ', endif rt < 2, error(' Table too small or not oriented in columns. ', endabove=ta

29、b(: , col > val; % True where > VALbelow=tab(: , col < val; % True where < VALequal=tab(: , col = = val; % True where = VALif all(above = = 0 | all(below = = 0, % handle simplest casey=tab(find(equal, : ; returnendpslope=find(below(1:rt-1&above(2:rt; % indices where slope is +nslope=

30、find(below(2:rt&above(1:rt-1; % indices where slope is -ib=sort(pslope; nslope+1; % put indices below in orderia=sort(nslope; pslope+1; % put indices above in orderie=find(equal; % indices where equal to valtmp,ix=sort( ib, ie ; % find where equals fit in resultieq=ix > length(ib; % True wher

31、e equals values fitry=length(tmp; % # of rows in result yy=zeros(ry, ct; % poke data into a zero matrixalpha=(val-tab(ib,col./(tab(ia,col-tab(ib,col;alpha=alpha(: , ones(1, ct; % duplicate for all columnsy(ieq, : =alpha.*tab(ia, : +(1-alpha.*tab(ib, : ; % interpolated valuesy(ieq, : =tab(ie, : ; % e

32、qual valuesy( : , col=val*ones(ry, 1; % remove roundoff error正如所見(jiàn)的，mminterp 利用了find 和sort 函數(shù)、邏輯數(shù)組和數(shù)組操作技術(shù)。沒(méi)有For 循環(huán)和While 循環(huán)。不論用其中哪一種技術(shù)來(lái)實(shí)現(xiàn)將使運(yùn)行變慢，尤其對(duì)大的表。注意mminterp 與含有大于或等于2的任意數(shù)列的表一起工作，如同函數(shù)interp1一樣。而且，在這種情況下，插值變量可以是任意的列。例如，» z=sin(pi*x; % add more data to table» table=x; y; z.' ;» t

33、=mminterp(table, 2, 1.1 % same interpolation as earliert =» t=mminterp(table, 3, -.5 % second third column nowt =這些最后的結(jié)果估計(jì)了x 和y 在z= -.5處的值。盡管逐條地對(duì)函數(shù)mminterp 解釋如何工作是很有幫助的，但這樣做要求有更多的篇幅和時(shí)間。解釋mminterp 如何工作最容易的方法是創(chuàng)建一個(gè)小表格，然后，在重要的語(yǔ)句末尾刪除分號(hào)以后，調(diào)用函數(shù)。這樣，中間值將幫助用戶(hù)理解函數(shù)是如何找到與所需值相符的數(shù)據(jù)值以及如何執(zhí)行插值。前面已闡述了interp1的用法。當(dāng)

34、用于線(xiàn)性插值時(shí)，只要所要求的插值點(diǎn)的個(gè)數(shù)少，interp1工作很好。在要求許多插值點(diǎn)情況下，由于所用的算法，interp1工作較慢。為了克服這個(gè)問(wèn)題，精通MATLAB 工具箱包括了函數(shù)mmtable ，它的幫助文本是：»help tableMMTABLE 1-D Monotonic Table Search by Linear Interpolation.YI=MMTABLE(TAB,COL,VALS linearly interpolates the table TABsearching for values VALS in the column COL.TAB(:,COL must be monotonic, but need NOT be equally spaced.YI has as many rows as VALS and as many columns TAB NaNs are returned where VALS are outside the range of TAB(:,