三次样条插值（CubicSplineInterpolation）及代码实现（C语
⾔）
样条插值是⼀种⼯业设计中常⽤的、得到平滑曲线的⼀种插值⽅法，三次样条⼜是其中⽤的较为⼴泛的⼀种。本篇介绍⼒求⽤容易理解的⽅式，介绍⼀下三次样条插值的原理，并附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，代码如下：
#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);
ssSetOfftTime(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的??
小学古诗80首
/
* 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* 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)
{
