二維曲線擬合
相關基礎理論知識
我們根據平面上的一些離散點繪製出一條近似曲線。如果曲線通過所有點，稱為插值；如果曲線不一定通過點，而是以某種方式逼近這些點，稱為擬合。
構造擬合曲線，通常有以下幾種方法：
（1）最小二乘法；
（2）分段擬合法；
（3）…
最小二乘法
根據解的存在情況，線性方程可以分為：
- 有唯一解的恰定方程組；
- 解不存在的超定方程組；
(Ax=b，A為n×m矩陣，如果A列滿秩，且n>m,方程組沒有精確解,常用於資料擬合)； - 有無窮多解的欠定方程組；
(Ax=b，A為n×m矩陣，如果A行滿秩，且n<m)；
在MATLAB中求解超定方程，有以下幾種方法：
- 左除是建立在奇異值分解基礎之上得到最小二乘法的解，因此最可靠；
x=A\b;
- 最小二乘法求解;
x=lsqnonneg(A,b);
- 廣義逆，解不一定滿足Ax=b，x只是最小二乘意義上的解;
x=pinv(A);
廣義逆（偽逆矩陣）
廣義逆法（偽逆矩陣）是建立在對原超定方程直接進行 householder變換的基礎上，其演算法可靠性稍遜與（SVD）奇異值求解，但速度較快；以B=pinv(A)為例，函式返回矩陣A的偽逆矩陣。如果矩陣A是可逆（非奇異）的，那麼pinv(A)與inv(A)的結果是一樣的，而且pinv比inv效率低。但如果矩陣A是非方陣或奇異矩陣，則inv(A)不存在，但pinv(A)仍然存在，並表現出一些與逆矩陣類似的性質。
【定義】
令A是任意m x n矩陣，若G滿足下述條件（Moore-penrose條件），稱矩陣G是A的廣義逆矩陣：
（1）GAG = G；
（2）AGA = A；
（3）AG為hermitian矩陣，即（AG）^H=AG；
（4）GA為hermitian矩陣，即（GA）^H=GA；
【測試】
在Matlab中，用以下幾種方式求逆：
（1）直接求解：InvA = inv(A’*A)*A’; %求導，令導數為0，結果如下: InvA=(ATA)-1AT
（2）SVD分解：[U,Λ,V]=svd(A)
（3）QR分解：[Q,R]=qr(A)
（4）LU分解：[L,U]=lu(A)
a= 結果：
矩陣分解
由上面例子可知，我們可以通過矩陣分解來求逆矩陣，接下來介紹幾種矩陣分解的方式：
-
特徵值分解（Eigen Value Decomposition，EVD）
特徵值分解，是指將矩陣分解為由特徵值和特徵向量表示的矩陣之積的方法。由（λ I - A）v=0（λ 為N 階方陣A的特徵值，v為特徵值 λ 對應的特徵向量）求解特徵值和特徵向量：
（1）特徵值：det | λ I - A | = 0（特徵多項式），可以得到：
這樣，我們便可以得到矩陣A的l種不同的特徵值，每種特徵值有ni個重複，共有N個特徵值。
（2）特徵向量：對於單根特徵值，（λi I - A）vi =0，即可求得該特徵值對應的特徵向量；對於有重根的特徵值，可用以下公式依次迭代求解：
這樣，特徵分解將矩陣分解為如下形式：A=QΛQ-1；若無重根，Λ=diag（λ1，λ2，…，λn）,若有重根，Λ為Jordan標準型。這裡的Q是以特徵向量為列向量的N x N矩陣。 -
Householder變換法（QR Factorization,QR分解）
這是目前求一般矩陣全部特徵值的最有效並廣泛應用的方法，精度高於偽逆法。一般矩陣先經過正交相似變化成為Heisenberg矩陣，然後再應用QR方法求特徵值和特徵向量。
A=QR （其中：A為m x n的矩陣，n>=m，Q為m×m的正交矩陣，R為m×n的上三角矩陣），我們稱其為A的QR分解。 -
奇異值分解（singular value decomposition,SVD）
對於矩陣A ϵ Fm x n，A=UΛVT，其中，U ϵ Fm x m酉矩陣，被稱為左奇異向量；Λ ϵ Fm x n的半正定對角矩陣，Λ 對角線上的元素 λi為原矩陣A的奇異值；V T ϵ F n×n 是V 的共軛轉置，被稱為右奇異向量，這樣的分解稱為奇異值分解。
Matlab樣條工具箱(Curve Fitting Toolbox)與曲線擬合
樣條函式可分為4類：cs* 三次樣條、係數為t^n 的pp* 分段多項式樣條、係數為基函式B_n^i(t)的sp* B樣條、 rp* 有理B樣條。樣條操作包括：函式操作（求值，算術運算，求導求積分等）、節點操作（節點重數的調節，設定，修改等）。
- 匯入資料，開啟工具箱選擇合適的模型擬合
程式碼：
x0=data(1,:);
y0=data(2,:);
cftool;
結果：
2. 樣條工具箱函式
（1）. 三次樣條函式
csapi 插值生成三次樣條函式
csape 生成給定約束條件下的三次樣條函式
csaps 平滑生成三次樣條函式
cscvn 生成一條內插引數的三次樣條曲線
getcurve 動態生成三次樣條曲線
（2）. 分段多項式樣條函式
ppmak 生成分段多項式樣條函式
ppual 計算在給定點處的分段多項式樣條函式值
（3）. B樣條函式
spmak 生成B樣條函式
spcrv 生成均勻劃分的B樣條函式
spapi 插值生成B樣條函式
spap2 用最小二乘法擬合生成B樣條函式
spaps 對生成的B樣條曲線進行光滑處理
spcol 生成B樣條函式的配置矩陣
（4）. 有理樣條函式
rpmak 生成有理樣條函式
rsmak 生成有理樣條函式
（5）. 樣條操作函式
fnval 計算在給定點處的樣條函式值
fmbrk 返回樣條函式的某一部分（如斷點或係數等）
fncmb 對樣條函式進行算術運算
fn2fm 把一種形式的樣條函式轉化成另一種形式的樣條函式
fnder 求樣條函式的微分(即求導數)
fndir 求樣條函式的方向導數
fnint 求樣條函式的積分
fnjmp 在間斷點處求函式值
fnplt 畫樣條曲線圖
fnrfn 在樣條曲線中插入斷點。
fntlr 生成tarylor係數或taylor多項式
（6）. 樣條曲線端點和節點處理函式
augknt 在已知節點陣列中新增一個或多個節點
aveknt 求出節點陣列元素的平均值
brk2knt 增加節點陣列中節點的重次
knt2brk 從節點陣列中求得節點及其重次
knt2mlt 從節點陣列中求得節點及其重次
sorted 求出節點陣列的元素在另一節點陣列中屬於第幾個分量
aptknt 求出用於生成樣條曲線的節點陣列
newknt 對分段多項式樣條函式進行重分佈
optknt 求出用於內插的最優節點陣列
chbpnt 求出用於生成樣條曲線的合適節點陣列
樣條線擬合
spline函式（piecewise polynomial function）
特點：
設樣條函式為S；匯入資料點為：knots（結）（k組資料，則有k個結，k-1段，k-1個多項式函式）；相鄰兩個knots之間分段定義多項式函式Si（i=0,1,……,k-1），則S=S1+ S2+……+Sk-1；
>>x0=data(1,:);
>>y0=data(2,:);
>>x1=1:0.1:7;%設定插值的自變數
>>n=spline(x0,y0);%返回函式的結構體資訊
>>y1=spline(x0,y0,x1);%返回插值計算的函式值
>>plot(x0,y0,'ro');
>>plot(x1,y1,'r-');
>>num=4;
>>x=0:0.1:7;
%最大值所在方程函式
>>y=f2.coefs(num,1).*(x-num).^3+f2.coefs(num,2).*(x-num).^2+f2.coefs(num,3).*(x-num)+f2.coefs(num,4); >>x_solve=x(find(y==max(y)));
結果：
n =
form: 'pp' %樣條函式型別
breaks: [1 2 3 4 5 6 7] %分段區間節點矩陣
coefs: [6x4 double] %各分段區間上的插值多項式係數矩陣
pieces: 6 %分段數
order: 4 %項式階數
dim: 1 %矩陣維度
x_solve =
4.3000
polyfit(x,y,N)函式
特點：對資料x，y擬合N階多項式
>>x0=data(1,:);
>>y0=data(2,:);
>>x1=1:0.1:n;
>>f1=polyfit(x0,y0,2); %polynomial fitting
>>y2=polyval(f1,x1); %Calculate dependent variable
>>plot(x1,y2,'g-');
結果：
二次B樣條曲線
特點：已知三個平面離散點就可以定義二次拋物線段，多於3個離散點，則採用二次B樣條曲線分段擬合。
設有n+1個離散點，記為Pi（i=0,1，…,n）
矩陣形式為：
三次B樣條曲線
矩陣形式為：
二次Beier曲線
矩陣形式為：
