第5章-MATLAB數(shù)據(jù)分析與多項(xiàng)式計(jì)算1

MATLAB程序設(shè)計(jì)教程（第二版）劉衛(wèi)國(guó)主編

中國(guó)水利水電出版社第5章

MATLAB數(shù)據(jù)分析與多項(xiàng)式計(jì)算

MATLAB數(shù)據(jù)統(tǒng)計(jì)處理

MATLAB數(shù)據(jù)插值與曲線擬合

MATLAB離散傅立葉變換

MATLAB多項(xiàng)式計(jì)算第5章MATLAB數(shù)據(jù)分析與多項(xiàng)式計(jì)算5.1數(shù)據(jù)統(tǒng)計(jì)處理5.2數(shù)據(jù)插值5.3曲線擬合5.4離散傅立葉變換5.5多項(xiàng)式計(jì)算5.1數(shù)據(jù)統(tǒng)計(jì)處理5.1.1最大值和最小值5.1.2求和與求積5.1.3平均值和中值5.1.4累加和與累乘積5.1.5標(biāo)準(zhǔn)方差與相關(guān)系數(shù)5.1.6排序5.1數(shù)據(jù)統(tǒng)計(jì)處理5.1.1最大值和最小值MATLAB提供的求數(shù)據(jù)序列的最大值和最小值的函數(shù)分別為max和min，兩個(gè)函數(shù)的調(diào)用格式和操作過程類似。1．求向量的最大值和最小值求一個(gè)向量X的最大值的函數(shù)有兩種調(diào)用格式，分別是：y=max(X)：返回向量X的最大值存入y，如果X中包含復(fù)數(shù)元素，則按模取最大值。(2)[y,I]=max(X)：返回向量X的最大值存入y，最大值的序號(hào)存入I，如果X中包含復(fù)數(shù)元素，則按模取最大值。求向量X的最小值的函數(shù)是min(X)，用法和max(X)完全相同。例5-1求向量x的最大值。x=[-43,72,9,16,23,47];y=max(x)%求向量x中的最大值[y,l]=max(x)%求向量x中的最大值及其該元素的位置以上是對(duì)行向量進(jìn)行操作，對(duì)列向量的操作是一樣的。如：y=max(x’)%求向量x中的最大值[y,l]=max(x’)2．求矩陣的最大值和最小值求矩陣A的最大值的函數(shù)有3種調(diào)用格式，分別是：(1)max(A)：返回一個(gè)行向量，向量的第i個(gè)元素是矩陣A的第i列上的最大值。(2)[Y,U]=max(A)：返回行向量Y和U，Y向量記錄A的每列的最大值，U向量記錄每列最大值的行號(hào)。(3)max(A,[],dim)：dim取1或2。dim取1時(shí)，該函數(shù)和max(A)完全相同；dim取2時(shí)，該函數(shù)返回一個(gè)列向量，其第i個(gè)元素是A矩陣的第i行上的最大值。求最小值的函數(shù)是min，其用法和max完全相同。例5-2分別求3×4矩陣x中各列和各行元素中的最大值，并求整個(gè)矩陣的最大值和最小值。X=[1842;9625;3671];Y=max(X)

%求矩陣X中各列元素的最大值[Y,l]=max(X)%求矩陣X中各列元素的最大值及其下標(biāo)[Y,l]=max(X,[],1)%同上[Y,l]=max(X,[],2)%求矩陣X中各行元素的最大值及其下標(biāo)max(max(X))%整個(gè)矩陣的最大值min(min(X))%整個(gè)矩陣的最小值3．兩個(gè)向量或矩陣對(duì)應(yīng)元素的比較函數(shù)max和min還能對(duì)兩個(gè)同型（同維）的向量或矩陣進(jìn)行比較，調(diào)用格式為：(1)U=max(A,B)：A,B是兩個(gè)同型的向量或矩陣，結(jié)果U是與A,B同型的向量或矩陣，U的每個(gè)元素等于A,B對(duì)應(yīng)元素的較大者。(2)U=max(A,n)：n是一個(gè)標(biāo)量，結(jié)果U是與A同型的向量或矩陣，U的每個(gè)元素等于A對(duì)應(yīng)元素和n中的較大者。min函數(shù)的用法和max完全相同。比較出a，b兩個(gè)同型距陣對(duì)應(yīng)元素的大小，較大的元素組成u距陣f是一個(gè)標(biāo)量，下面是距陣u與f的比較結(jié)果5.1.2求和與求積數(shù)據(jù)序列求和與求積的函數(shù)是sum和prod，其使用方法類似。設(shè)X是一個(gè)向量，A是一個(gè)矩陣，函數(shù)的調(diào)用格式為：sum(X)：求向量X各元素的和。prod(X)：求回向量X各元素的乘積。sum(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的元素和。prod(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的元素乘積。求距陣a的列元素和與積sum(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于sum(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的各元素之和。prod(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于prod(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的各元素乘積。各行的乘積作為s的元素求距陣a的所有元素之積5.1.3平均值和中值求數(shù)據(jù)序列平均值的函數(shù)是mean，求數(shù)據(jù)序列中值的函數(shù)是median。兩個(gè)函數(shù)的調(diào)用格式為：mean(X)：返回向量X的算術(shù)平均值。median(X)：返回向量X的中值.mean(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的算術(shù)平均值。median(A)：返回一個(gè)行向量，其第i個(gè)元素是A的第i列的中值。mean(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于mean(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的算術(shù)平均值。median(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于median(A)；當(dāng)dim為2時(shí)，返回一個(gè)列向量，其第i個(gè)元素是A的第i行的中值。例5-5分別求向量x與y的平均值和中值。求向量x的算術(shù)平均值Median命令可生成x的中值5.1.4累加和與累乘積在MATLAB中，使用cumsum和cumprod函數(shù)能方便地求得向量和矩陣元素的累加和與累乘積向量：累加和：向量X=(x1,x2,…,xn)的累加和生成一個(gè)同維向量Y=(y1,y2,…,yn),其中yi=x1+x2+…+xi累乘積：向量X=(x1,x2,…,xn)的累加和生成一個(gè)同維向量Y=(y1,y2,…,yn),其中yi=x1*x2*…*xicumsum(X)：返回向量X累加和向量。cumprod(X)：返回向量X累乘積向量。cumsum(A)：返回一個(gè)矩陣，其第i列是A的第i列的累加和向量。cumprod(A)：返回一個(gè)矩陣，其第i列是A的第i列的累乘積向量。cumsum(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于cumsum(A)；當(dāng)dim為2時(shí)，返回一個(gè)矩陣，其第i行是A的第i行的累加和向量。cumprod(A,dim)：當(dāng)dim為1時(shí)，該函數(shù)等同于cumprod(A)；當(dāng)dim為2時(shí)，返回一個(gè)向量，其第i行是A的第i行的累乘積向量。求向量x的累加和與累乘積5.1.5標(biāo)準(zhǔn)方差與相關(guān)系數(shù)1．求標(biāo)準(zhǔn)方差向量x=(x1,x2,…,xn),標(biāo)準(zhǔn)方差有兩個(gè)計(jì)算公式：或在MATLAB中，提供了計(jì)算數(shù)據(jù)序列的標(biāo)準(zhǔn)方差的函數(shù)std。對(duì)于向量X，std(X)返回一個(gè)標(biāo)準(zhǔn)方差。對(duì)于矩陣A，std(A)返回一個(gè)行向量，它的各個(gè)元素便是矩陣A各列或各行的標(biāo)準(zhǔn)方差。std函數(shù)的一般調(diào)用格式為：

Y=std(A,flag,dim)其中dim取1或2。當(dāng)dim=1時(shí)，求各列元素的標(biāo)準(zhǔn)方差；當(dāng)dim=2時(shí)，則求各行元素的標(biāo)準(zhǔn)方差。flag取0或1，當(dāng)flag=0時(shí)，按σ1所列公式計(jì)算標(biāo)準(zhǔn)方差，當(dāng)flag=1時(shí)，按σ2所列公式計(jì)算標(biāo)準(zhǔn)方差。缺省flag=0，dim=1。例5-7對(duì)二維矩陣x，從不同維方向求出其標(biāo)準(zhǔn)方差。x=[456;148];y1=std(x,0,1)y2=std(x,1,1)y3=std(x,0,2)y2=std(x,1,2)2．相關(guān)系數(shù)MATLAB提供了corrcoef函數(shù)，可以求出數(shù)據(jù)的相關(guān)系數(shù)矩陣。corrcoef函數(shù)的調(diào)用格式為：corrcoef(X)：返回從矩陣X形成的一個(gè)相關(guān)系數(shù)矩陣。此相關(guān)系數(shù)矩陣是方陣，大小與矩陣X的列數(shù)一樣。它把矩陣X的每列作為一個(gè)變量，然后求它們的相關(guān)系數(shù)。corrcoef(X,Y)：在這里，X,Y是向量，它們與corrcoef([X,Y])的作用一樣。例5-8生成滿足正態(tài)分布的10000×5隨機(jī)矩陣，然后求各列元素的均值和標(biāo)準(zhǔn)方差，再求這5列隨機(jī)數(shù)據(jù)的相關(guān)系數(shù)矩陣。命令如下：X=randn(10000,5);M=mean(X)D=std(X)R=corrcoef(X)5.1.6排序MATLAB中對(duì)向量X是排序函數(shù)是sort(X)，函數(shù)返回一個(gè)對(duì)X中的元素按升序排列的新向量。sort函數(shù)也可以對(duì)矩陣A的各列或各行重新排序，其調(diào)用格式為：

[Y,I]=sort(A,dim,mode)其中dim指明對(duì)A的列還是行進(jìn)行排序,默認(rèn)則按列排；若dim=1，則按列排；若dim=2，則按行排。Y是排序后的矩陣，而I記錄Y中的元素在A中位置。mode是指明按升序還是降序排列，默認(rèn)為升序，

'ascend'是升序，'descend'是降序。對(duì)距陣x的各列按升序排列對(duì)距陣x的各行按升序排列5.2數(shù)據(jù)插值在生產(chǎn)和科學(xué)實(shí)驗(yàn)中，我們不知道函數(shù)f(x)的具體表達(dá)式，只能獲得它的若干點(diǎn)的函數(shù)值或微商值，即只給出的f(x)一張數(shù)據(jù)表。如果根據(jù)這張數(shù)據(jù)表，構(gòu)造一個(gè)簡(jiǎn)單函數(shù)g(x)，使之滿足數(shù)據(jù)表中的數(shù)據(jù)。這樣的函數(shù)g(x)就是的逼近函數(shù)f(x)。這種逼近問題就稱為插值問題。

xy機(jī)翼下輪廓線例已知飛機(jī)下輪廓線上數(shù)據(jù)如下，求x每改變0.1時(shí)的y值。設(shè)未知函數(shù)為f(x),求函數(shù)g(x)去逼近函數(shù)f(x)，要求：使得g(x)與f(x)在上面對(duì)應(yīng)位置x上的值相等，即若x=11,則g(x)=f(x)=2.0一維插值函數(shù)：y1=interp1(x，y，x1，'method')插值方法被插值點(diǎn)插值節(jié)點(diǎn)(已知點(diǎn))‘nearest’

：最鄰近插值‘linear’

：線性插值；‘spline’

：三次樣條插值；‘cubic’

：立方插值。缺省時(shí)：分段線性插值。注意：所有的插值方法都要求x是單調(diào)的，并且x1不能夠超過x的范圍。5.2.1一維數(shù)據(jù)插值‘X,Y是兩個(gè)等長(zhǎng)的已知向量，分別描述采樣點(diǎn)和樣本值，X1是一個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)，Y1是一個(gè)與X1等長(zhǎng)的插值結(jié)果。例：在1-12的11小時(shí)內(nèi)，每隔1小時(shí)測(cè)量一次溫度，測(cè)得的溫度依次為：5，8，9，15，25，29，31，30，22，25，27，24。試估計(jì)每隔1/10小時(shí)的溫度值。hours=1:12;temps=[589152529313022252724];h=1:0.1:12;t=interp1(hours,temps,h,'spline');%(直接輸出數(shù)據(jù)將是很多的)plot(hours,temps,'+',h,t,hours,temps,'r:')%作圖xlabel('Hour'),ylabel('DegreesCelsius')注意：X1的取值范圍不能超出X的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。例5-10用不同的插值方法計(jì)算在π/2點(diǎn)的值。

x=0:0.2:pi;y=sin(x);interp1(x,y,pi/2,’linear’)interp1(x,y,pi/2,’nearest’)interp1(x,y,pi/2,’spline’)interp1(x,y,pi/2,’cubic’)MATLAB中有一個(gè)專門的3次樣條插值函數(shù)Y1=spline(X,Y,X1)，其功能及使用方法與函數(shù)Y1=interp1(X,Y,X1,‘spline’)完全相同。例5-11某觀測(cè)站測(cè)得某日6:00時(shí)至18:00時(shí)之間每隔2小時(shí)的室內(nèi)外溫度(℃)，用3次樣條插值分別求得該日室內(nèi)外6:30至17:30時(shí)之間每隔2小時(shí)各點(diǎn)的近似溫度(℃)。設(shè)時(shí)間變量h為一行向量，溫度變量t為一個(gè)兩列矩陣，其中第一列存放室內(nèi)溫度，第二列儲(chǔ)存室外溫度。命令如下：h=6:2:18;t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';XI=6.5:2:17.5YI=interp1(h,t,XI,‘spline’)%用3次樣條插值計(jì)算5.2.2二維數(shù)據(jù)插值在MATLAB中，提供了解決二維插值問題z=f(x,y)的函數(shù)interp2，其調(diào)用格式為：

Z1=interp2(X,Y,Z,X1,Y1,'method')其中X,Y是兩個(gè)向量，分別描述兩個(gè)參數(shù)的采樣點(diǎn)，Z是與參數(shù)采樣點(diǎn)對(duì)應(yīng)的函數(shù)值，X1,Y1是兩個(gè)向量或標(biāo)量，描述欲插值的點(diǎn)。Z1是根據(jù)相應(yīng)的插值方法得到的插值結(jié)果。method的取值與一維插值函數(shù)相同。X,Y,Z也可以是矩陣形式。同樣，X1,Y1的取值范圍不能超出X,Y的給定范圍，否則，會(huì)給出“NaN”錯(cuò)誤。例5-13某實(shí)驗(yàn)對(duì)一根長(zhǎng)10米的鋼軌進(jìn)行熱源的溫度傳播測(cè)試。用x表示測(cè)量點(diǎn)0:2.5:10(米)，用h表示測(cè)量時(shí)間0:30:60(秒)，用T表示測(cè)試所得各點(diǎn)的溫度(℃)。試用線性插值求出在一分鐘內(nèi)每隔20秒、鋼軌每隔1米處的溫度TI。命令如下：x=0:2.5:10;h=[0:30:60]';T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];meshz(x,h,T)xi=[0:10];hi=[0:20:60]';TI=interp2(x,h,T,xi,hi)figure(2),meshz(xi,hi,TI)gridonholdon例：測(cè)得平板表面3*5網(wǎng)格點(diǎn)處的溫度分別為：828180828479636165818484828586試作出平板表面的溫度分布曲面z=f(x,y)的圖形。輸入以下命令：x=1:5;y=1:3;temps=[8281808284;7963616581;8484828586];mesh(x,y,temps)1.先在三維坐標(biāo)畫出原始數(shù)據(jù)，畫出粗糙的溫度分布曲圖.2．以平滑數(shù)據(jù),在x、y方向上每隔0.2個(gè)單位的地方進(jìn)行插值.再輸入以下命令:xi=1:0.2:5;yi=1:0.2:3;zi=interp2(x,y,temps,xi',yi,'cubic');mesh(xi,yi,zi)畫出插值后的溫度分布曲面圖.

插值函數(shù)griddata格式為:

cz

=griddata（x，y，z，cx，cy，‘method’）用MATLAB作散點(diǎn)數(shù)據(jù)的插值計(jì)算

要求cx取行向量，cy取為列向量。被插值點(diǎn)插值方法插值節(jié)點(diǎn)被插值點(diǎn)的函數(shù)值‘nearest’最鄰近插值‘linear’雙線性插值‘cubic’雙三次插值'v4'-Matlab提供的插值方法缺省時(shí),雙線性插值例在某海域測(cè)得一些點(diǎn)(x,y)處的水深z由下表給出，船的吃水深度為5英尺，在矩形區(qū)域（75，200）*（-50，150）里的哪些地方船要避免進(jìn)入。clearx=[129140103.588185.5195105.5157.5107.57781162162117.5];y=[7.5141.52314722.5137.585.5-6.5-81356.5-66.584-33.5];z=[-4-8-6-8-6-8-8-9-9-8-8-9-4-9];cx=75:0.5:200;cy=-70:0.5:150;cz=griddata(x,y,z,cx,cy','cubic');4.作出水深小于5的海域范圍,即z=-5的等高線...3作海底曲面圖meshz(cx,cy,cz),rotate3d%3d旋轉(zhuǎn)xlabel('X'),ylabel('Y'),zlabel('Z')%pausefigure(2),contour3(cx,cy,cz,[-5-5],'k');gridon

%只畫z=-5的等高線holdonplot3(x,y,z,'+')xlabel('X'),ylabel('Y')擬合問題引例1溫度t(0C)20.532.751.073.095.7電阻R()7658268739421032已知熱敏電阻數(shù)據(jù)：求600C時(shí)的電阻R。

設(shè)

R=at+ba,b為待定系數(shù)5.3曲線擬合擬合問題引例2

t(h)0.250.511.523468c(g/ml)19.2118.1515.3614.1012.899.327.455.243.01已知一室模型快速靜脈注射下的血藥濃度數(shù)據(jù)(t=0注射300mg)求血藥濃度隨時(shí)間的變化規(guī)律c(t).作半對(duì)數(shù)坐標(biāo)系(semilogy)下的圖形曲線擬合問題的提法已知一組（二維）數(shù)據(jù)，即平面上n個(gè)點(diǎn)（xi,yi)i=1,…n,

尋求一個(gè)函數(shù)（曲線）y=f(x),

使f(x)

在某種準(zhǔn)則下與所有數(shù)據(jù)點(diǎn)最為接近，即曲線擬合得最好。

+++++++++xyy=f(x)(xi,yi)

i

i為點(diǎn)（xi,yi)與曲線y=f(x)的距離擬合與插值的關(guān)系

函數(shù)插值與曲線擬合都是要根據(jù)一組數(shù)據(jù)構(gòu)造一個(gè)函數(shù)作為近似，由于近似的要求不同，二者的數(shù)學(xué)方法上是完全不同的。實(shí)例：下面數(shù)據(jù)是某次實(shí)驗(yàn)所得，希望得到X和f之間的關(guān)系？問題：給定一批數(shù)據(jù)點(diǎn)，需確定滿足特定要求的曲線或曲面解決方案：若不要求曲線（面）通過所有數(shù)據(jù)點(diǎn)，而是要求它反映對(duì)象整體的變化趨勢(shì)，這就是數(shù)據(jù)擬合，又稱曲線擬合或曲面擬合。若要求所求曲線（面）通過所給所有數(shù)據(jù)點(diǎn)，就是插值問題；最臨近插值、線性插值、樣條插值與曲線擬合結(jié)果：用MATLAB解擬合問題1、線性最小二乘擬合2、非線性最小二乘擬合用MATLAB作線性最小二乘擬合1.作多項(xiàng)式f(x)=a1xm+…+amx+am+1擬合,可利用已有程序:a=polyfit(x,y,m)2.對(duì)超定方程組可得最小二乘意義下的解。，用3.多項(xiàng)式在x處的值y可用以下命令計(jì)算：

y=polyval（a，x）輸出擬合多項(xiàng)式系數(shù)a=[a1,…am,

am+1](數(shù)組)）輸入同長(zhǎng)度的數(shù)組X，Y擬合多項(xiàng)式次數(shù)即要求出二次多項(xiàng)式:中的使得:例對(duì)下面一組數(shù)據(jù)作二次多項(xiàng)式擬合1）輸入以下命令：

x=0:0.1:1;y=[-0.4471.9783.286.167.087.347.669.569.489.3011.2];A=polyfit(x,y,2)z=polyval(A,x);plot(x,y,'k+',x,z,'r')%作出數(shù)據(jù)點(diǎn)和擬合曲線的圖形2）計(jì)算結(jié)果：Ａ=-9.810820.1293-0.0317用多項(xiàng)式擬合的命令1.lsqcurvefit已知數(shù)據(jù)點(diǎn)：

xdata=（xdata1，xdata2，…，xdatan），

ydata=（ydata1，ydata2，…，ydatan）

用MATLAB作非線性最小二乘擬合

Matlab的提供了兩個(gè)求非線性最小二乘擬合的函數(shù)：lsqcurvefit和lsqnonlin。兩個(gè)命令都要先建立M-文件fun.m，在其中定義函數(shù)f(x)，但兩者定義f(x)的方式是不同的,可參考例題.

lsqcurvefit用以求含參量x（向量）的向量值函數(shù)F(x,xdata)=（F（x，xdata1），…，F(xiàn)（x，xdatan））T中的參變量x(向量),使得

輸入格式為:(1)x=lsqcurvefit(‘fun’,x0,xdata,ydata);(2)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);(3)x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);(4)[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(5)[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);(6)[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);fun是一個(gè)事先建立的定義函數(shù)F(x,xdata)

的M-文件,自變量為x和xdata說明：x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);迭代初值已知數(shù)據(jù)點(diǎn)選項(xiàng)見無約束優(yōu)化

lsqnonlin用以求含參量x（向量）的向量值函數(shù)

f(x)=(f1(x),f2(x),…,fn(x))T

中的參量x，使得

最小。其中fi（x）=f（x，xdatai，ydatai）

=F(x,xdatai)-ydatai

2.lsqnonlin已知數(shù)據(jù)點(diǎn)：xdata=（xdata1，xdata2，…，xdatan）

ydata=（ydata1，ydata2，…，ydatan）輸入格式為：

1）x=lsqnonlin（‘fun’，x0）；

2）x=lsqnonlin

（‘fun’，x0，options）；

3）x=lsqnonlin

（‘fun’，x0，options，‘grad’）；

4）[x，options]=lsqnonlin

（‘fun’，x0，…）；

5）[x,options,funval]=lsqnonlin(‘fun’，x0，…)；說明：x=lsqnonlin

（‘fun’，x0，options）；fun是一個(gè)事先建立的定義函數(shù)f(x)的M-文件，自變量為x迭代初值選項(xiàng)見無約束優(yōu)化

例2用下面一組數(shù)據(jù)擬合中的參數(shù)a，b，k該問題即解最優(yōu)化問題：

1）編寫M-文件curvefun1.mfunctionf=curvefun1(x,tdata)f=x(1)+x(2)*exp(-0.02*x(3)*tdata)%其中x(1)=a;x(2)=b；x(3)=k;2）輸入命令tdata=100:100:1000cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];x0=[0.2,0.05,0.05];x=lsqcurvefit('curvefun1',x0,tdata,cdata)f=curvefun1(x,tdata)

F(x，tdata)=

，x=(a，b，k)解法1.用命令lsqcurvefit3）運(yùn)算結(jié)果為：f=0.00430.00510.00560.00590.00610.00620.00620.00630.00630.0063x=0.0063-0.00340.25424）結(jié)論：a=0.0063,b=-0.0034,k=0.2542

解法2

用命令lsqnonlin

f(x)=F(x,tdata,ctada)=x=（a，b，k）1）編寫M-文件curvefun2.mfunctionf=curvefun2(x)

tdata=100:100:1000;

cdata=1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];f=x(1)+x(2)*exp(-0.02*x(3)*tdata)-cdata2）輸入命令:

x0=[0.2,0.05,0.05];x=lsqnonlin('curvefun2',x0)f=curvefun2(x)函數(shù)curvefun2的自變量是x，cdata和tdata是已知參數(shù)，故應(yīng)將cdata

tdata的值寫在curvefun2.m中3）運(yùn)算結(jié)果為

f=1.0e-003*(0.2322-0.1243-0.2495-0.2413-0.1668-0.07240.02410.11590.20300.2792x=0.0063-0.00340.2542可以看出,兩個(gè)命令的計(jì)算結(jié)果是相同的.4）結(jié)論：即擬合得a=0.0063b=-0.0034k=0.25425.4離散傅立葉變換5.4.1離散傅立葉變換算法簡(jiǎn)要在某時(shí)間片等距地抽取N個(gè)抽樣時(shí)間tm處得樣本值f(tm)，記做f(m),m=0,1,2...N-1,稱向量F(k)(k=0,12...N-1)為f(m)的一個(gè)離散傅里葉變換，因在MATLAB中不允許有零下標(biāo)，m下標(biāo)均移動(dòng)1，公式更改為變換對(duì)6.4.2離散傅立葉變換的實(shí)現(xiàn)一維離散傅立葉變換函數(shù)，其調(diào)用格式與功能為：(1)fft(X)：返回向量X的離散傅立葉變換。設(shè)X的長(zhǎng)度(即元素個(gè)數(shù))為N，若N為2的冪次，則為以2為基數(shù)的快速傅立葉變換，否則為運(yùn)算速度很慢的非2冪次的算法。對(duì)于矩陣X，fft(X)應(yīng)用于矩陣的每一列。(2)fft(X,N)：計(jì)算N點(diǎn)離散傅立葉變換。它限定向量的長(zhǎng)度為N，若X的長(zhǎng)度小于N，則不足部分補(bǔ)上零；若大于N，則刪去超出N的那些元素。對(duì)于矩陣X，它同樣應(yīng)用于矩陣的每一列，只是限定了向量的長(zhǎng)度為N。(3)fft(X,[],dim)或fft(X,N,dim)：這是對(duì)于矩陣而言的函數(shù)調(diào)用格式，前者的功能與FFT(X)基本相同，而后者則與FFT(X,N)基本相同。只是當(dāng)參數(shù)dim=1時(shí)，該函數(shù)作用于X的每一列；當(dāng)dim=2時(shí)，則作用于X的每一行。

值得一提的是，當(dāng)已知給出的樣本數(shù)N0不是2的冪次時(shí)，可以取一個(gè)N使它大于N0且是2的冪次，然后利用函數(shù)格式fft(X,N)或fft(X,N,dim)便可進(jìn)行快速傅立葉變換。這樣，計(jì)算速度將大大加快。相應(yīng)地，一維離散傅立葉逆變換函數(shù)是ifft。ifft(F)返回F的一維離散傅立葉逆變換；ifft(F,N)為N點(diǎn)逆變換；ifft(F,[],dim)或ifft(F,N,dim)則由N或dim確定逆變換的點(diǎn)數(shù)或操作方向。例6-15給定數(shù)學(xué)函數(shù)x(t)=12sin(2π×10t+π/4)+5cos(2π×40t)取N=128，試對(duì)t從0~1秒采樣，用fft作快速傅立葉變換，繪制相應(yīng)的振幅-頻率圖。在0~1秒時(shí)間范圍內(nèi)采樣128點(diǎn)，從而可以確定采樣周期和采樣頻率。由于離散傅立葉變換時(shí)的下標(biāo)應(yīng)是從0到N-1，故在實(shí)際應(yīng)用時(shí)下標(biāo)應(yīng)該前移1。又考慮到對(duì)離散傅立葉變換來說，其振幅|F(k)|是關(guān)于N/2對(duì)稱的，故只須使k從0到N/2即可。程序如下：N=128;%采樣點(diǎn)數(shù)T=1;%采樣時(shí)間終點(diǎn)t=linspace(0,T,N);%給出N個(gè)采樣時(shí)間ti(I=1:N)x=12*sin(2*pi*10*t+pi/4)+5*cos(2*pi*40*t);%求各采樣點(diǎn)樣本值xdt=t(2)-t(1);%采樣周期f=1/dt;%采樣頻率(Hz)X=fft(x);%計(jì)算x的快速傅立葉變換XF=X(1:N/2+1);%F(k)=X(k)(k=1:N/2+1)f=f*(0:N/2)/N;%使頻率軸f從零開始plot(f,abs(F),'-*')%繪制振幅-頻率圖xlabel('Frequency');ylabel('|F(k)|')6.5多項(xiàng)式計(jì)算P(x)=a0x(n)+a1x(n-1)+.....a(n-1)x+anA=[a0a1aan]6.5.1多項(xiàng)式的四則運(yùn)算1．多項(xiàng)式的加減運(yùn)算(見課本p146)加法：a=[006-1];b=[006-1]c=a+b2．多項(xiàng)式乘法運(yùn)算函數(shù)conv(P1,P2)用于求多項(xiàng)式P1和P2的乘積。這里，P1、P2是兩個(gè)多項(xiàng)式系數(shù)向量。例6-16求多項(xiàng)式x4+8x3-10與多項(xiàng)式2x2-x+3的乘積。其運(yùn)算結(jié)果為2x6+15x5-5x4+24x3-20x2+10x-303．多項(xiàng)式除法函數(shù)[Q,r]=deconv(P1,P2)用于對(duì)多項(xiàng)式P1和P2作除法運(yùn)算。其中Q返回多項(xiàng)式P1除以P2的商式，r返回P1除以P2的余式。這里