樣條插值是一種工業設計中常用的、得到平滑曲線的一種插值方法，三次樣條又是其中用的較為廣泛的一種。本篇介紹力求用容易理解的方式，介紹一下三次樣條插值的原理，並附C語言的實現代碼。
1. 三次樣條曲線原理
假設有以下節點
1.1 定義
樣條曲線
是一個分段定義的公式。給定n+1個數據點，共有n個區間，三次樣條方程滿足以下條件：
a. 在每個分段區間
（i = 0, 1, …, n-1，x遞增）， 
都是一個三次多項式。
b. 滿足
（i = 0, 1, …, n ）
c. 
，導數
，二階導數
在[a, b]區間都是連續的，即
曲線是光滑的。
所以n個三次多項式分段可以寫作：
，i = 0, 1, …, n-1
其中ai, bi, ci, di代表4n個未知系數。
1.2 求解
已知:
a. n+1個數據點[xi, yi], i = 0, 1, …, n
b. 每一分段都是三次多項式函數曲線
c. 節點達到二階連續
d. 左右兩端點處特性（自然邊界，固定邊界，非節點邊界）
根據定點，求出每段樣條曲線方程中的系數，即可得到每段曲線的具體表達式。
插值和連續性：
, 其中 i = 0, 1, …, n-1
微分連續性：
, 其中 i = 0, 1, …, n-2
樣條曲線的微分式：
將步長
帶入樣條曲線的條件：
a. 由
(i = 0, 1, …, n-1)推出
b. 由
(i = 0, 1, …, n-1)推出
c. 由 
(i = 0, 1, …, n-2)推出
由此可得：
d. 由 
(i = 0, 1, …, n-2)推出
設
，則
a. 
可寫為：
，推出
b. 將ci, di帶入 
可得：
c. 將bi, ci, di帶入
(i = 0, 1, …, n-2)可得：
端點條件
由i的取值范圍可知，共有n-1個公式， 但卻有n+1個未知量m 。要想求解該方程組，還需另外兩個式子。所以需要對兩端點x0和xn的微分加些限制。 選擇不是唯一的，3種比較常用的限制如下。
a. 自由邊界(Natural)
首尾兩端沒有受到任何讓它們彎曲的力，即
。具體表示為
和 
則要求解的方程組可寫為：
b. 固定邊界(Clamped)
首尾兩端點的微分值是被指定的，這里分別定為A和B。則可以推出
將上述兩個公式帶入方程組，新的方程組左側為
c. 非節點邊界(Not-A-Knot)
指定樣條曲線的三次微分匹配，即
根據
和
，則上述條件變為
新的方程組系數矩陣可寫為：
右下圖可以看出不同的端點邊界對樣條曲線的影響：
1.3 算法總結
假定有n+1個數據節點
a. 計算步長
(i = 0, 1, …, n-1)
b. 將數據節點和指定的首位端點條件帶入矩陣方程
c. 解矩陣方程，求得二次微分值
。該矩陣為三對角矩陣，具體求法參見我的上篇文章：三對角矩陣的求解。
d. 計算樣條曲線的系數：
其中i = 0, 1, …, n-1
e. 在每個子區間
中，創建方程
2. C語言實現
用C語言寫了一個三次樣條插值(自然邊界)的S-Function，代碼如下：
View Code
#define S_FUNCTION_NAME cubic #define S_FUNCTION_LEVEL 2 #include "simstruc.h" #include "malloc.h" //方便使用變量定義數組大小 static void mdlInitializeSizes(SimStruct *S) { /*參數只有一個，是n乘2的定點數組[xi, yi]: * [ x1，y1; * x2, y2; * ..., ...; * xn, yn; */ ssSetNumSFcnParams(S, 1); if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return; ssSetNumContStates(S, 0); ssSetNumDiscStates(S, 0); if (!ssSetNumInputPorts(S, 1)) return; //輸入是x ssSetInputPortWidth(S, 0, 1); ssSetInputPortRequiredContiguous(S, 0, true); ssSetInputPortDirectFeedThrough(S, 0, 1); if (!ssSetNumOutputPorts(S, 1)) return; //輸出是S(x) ssSetOutputPortWidth(S, 0, 1); ssSetNumSampleTimes(S, 1); ssSetNumRWork(S, 0); ssSetNumIWork(S, 0); ssSetNumPWork(S, 0); ssSetNumModes(S, 0); ssSetNumNonsampledZCs(S, 0); ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE); ssSetOptions(S, 0); } static void mdlInitializeSampleTimes(SimStruct *S) { ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME); ssSetOffsetTime(S, 0, 0.0); } #define MDL_INITIALIZE_CONDITIONS #if defined(MDL_INITIALIZE_CONDITIONS) static void mdlInitializeConditions(SimStruct *S) { } #endif #define MDL_START #if defined(MDL_START) static void mdlStart(SimStruct *S) { } #endif /* MDL_START */ static void mdlOutputs(SimStruct *S, int_T tid) { const real_T *map = mxGetPr(ssGetSFcnParam(S,0)); //獲取定點數據 const int_T *mapSize = mxGetDimensions(ssGetSFcnParam(S,0)); //定點數組維數 const real_T *x = (const real_T*) ssGetInputPortSignal(S,0); //輸入x real_T *y = ssGetOutputPortSignal(S,0); //輸出y int_T step = 0; //輸入x在定點數中的位置 int_T i; real_T yval; for (i = 0; i < mapSize[0]; i++) { if (x[0] >= map[i] && x[0] < map[i + 1]) { step = i; break; } } cubic_getval(&yval, mapSize, map, x[0], step); y[0] = yval; } //自然邊界的三次樣條曲線函數 void cubic_getval(real_T* y, const int_T* size, const real_T* map, const real_T x, const int_T step) { int_T n = size[0]; //曲線系數 real_T* ai = (real_T*)malloc(sizeof(real_T) * (n-1)); real_T* bi = (real_T*)malloc(sizeof(real_T) * (n-1)); real_T* ci = (real_T*)malloc(sizeof(real_T) * (n-1)); real_T* di = (real_T*)malloc(sizeof(real_T) * (n-1)); real_T* h = (real_T*)malloc(sizeof(real_T) * (n-1)); //x的?? /* M矩陣的系數 *[B0, C0, ... *[A1, B1, C1, ... *[0, A2, B2, C2, ... *[0, ... An-1, Bn-1] */ real_T* A = (real_T*)malloc(sizeof(real_T) * (n-2)); real_T* B = (real_T*)malloc(sizeof(real_T) * (n-2)); real_T* C = (real_T*)malloc(sizeof(real_T) * (n-2)); real_T* D = (real_T*)malloc(sizeof(real_T) * (n-2)); //等號右邊的常數矩陣 real_T* E = (real_T*)malloc(sizeof(real_T) * (n-2)); //M矩陣 real_T* M = (real_T*)malloc(sizeof(real_T) * (n)); //包含端點的M矩陣 int_T i; //計算x的步長 for ( i = 0; i < n -1; i++) { h[i] = map[i + 1] - map[i]; } //指定系數 for( i = 0; i< n - 3; i++) { A[i] = h[i]; //忽略A[0] B[i] = 2 * (h[i] + h[i+1]); C[i] = h[i+1]; //忽略C(n-1) } //指定常數D for (i = 0; i<n - 3; i++) { D[i] = 6 * ((map[n + i + 2] - map[n + i + 1]) / h[i + 1] - (map[n + i + 1] - map[n + i]) / h[i]); } //求解三對角矩陣，結果賦值給E TDMA(E, n-3, A, B, C, D); M[0] = 0; //自然邊界的首端M為0 M[n-1] = 0; //自然邊界的末端M為0 for(i=1; i<n-1; i++) { M[i] = E[i-1]; //其它的M值 } //?算三次?條曲?的系數 for( i = 0; i < n-1; i++) { ai[i] = map[n + i]; bi[i] = (map[n + i + 1] - map[n + i]) / h[i] - (2 * h[i] * M[i] + h[i] * M[i + 1]) / 6; ci[i] = M[i] / 2; di[i] = (M[i + 1] - M[i]) / (6 * h[i]); } *y = ai[step] + bi[step]*(x - map[step]) + ci[step] * (x - map[step]) * (x - map[step]) + di[step] * (x - map[step]) * (x - map[step]) * (x - map[step]); free(h); free(A); free(B); free(C); free(D); free(E); free(M); free(ai); free(bi); free(ci); free(di); } void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D) { int_T i; real_T tmp; //上三角矩陣 C[0] = C[0] / B[0]; D[0] = D[0] / B[0]; for(i = 1; i<n; i++) { tmp = (B[i] - A[i] * C[i-1]); C[i] = C[i] / tmp; D[i] = (D[i] - A[i] * D[i-1]) / tmp; } //直接求出X的最后一個值 X[n-1] = D[n-1]; //逆向迭代， 求出X for(i = n-2; i>=0; i--) { X[i] = D[i] - C[i] * X[i+1]; } } #define MDL_UPDATE #if defined(MDL_UPDATE) static void mdlUpdate(SimStruct *S, int_T tid) { } #endif #define MDL_DERIVATIVES #if defined(MDL_DERIVATIVES) static void mdlDerivatives(SimStruct *S) { } #endif static void mdlTerminate(SimStruct *S) { } #ifdef MATLAB_MEX_FILE #include "simulink.c" #else #include "cg_sfun.h" #endif
3. 例子
以y=sin(x)為例, x步長為1，x取值范圍是[0,10]。對它使用三次樣條插值，插值前后對比如下：
