Go语言视频零基础入门到精通

admin

样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法，三次样条又是其中用的较为广泛的一种。本篇介绍力求用容易理解的方式，介绍一下三次样条插值的原理，并附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

1. #define S_FUNCTION_NAME  cubic
2. #define S_FUNCTION_LEVEL 2
3. #include "simstruc.h"
4. #include "malloc.h"  //方便使用变量定义数组大小
6. static void mdlInitializeSizes(SimStruct *S)
7. {
8.     /*参数只有一个，是n乘2的定点数组[xi, yi]:
9.      * [ x1，y1;
10.      *   x2, y2;
11.      *   ..., ...;
12.      *   xn, yn;
13.     */
14.     ssSetNumSFcnParams(S, 1);
15.     if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
17.     ssSetNumContStates(S, 0);
18.     ssSetNumDiscStates(S, 0);
20.     if (!ssSetNumInputPorts(S, 1)) return;  //输入是x
21.     ssSetInputPortWidth(S, 0, 1);
22.     ssSetInputPortRequiredContiguous(S, 0, true);
23.     ssSetInputPortDirectFeedThrough(S, 0, 1);
25.     if (!ssSetNumOutputPorts(S, 1)) return;  //输出是S(x)
26.     ssSetOutputPortWidth(S, 0, 1);
28.     ssSetNumSampleTimes(S, 1);
29.     ssSetNumRWork(S, 0);
30.     ssSetNumIWork(S, 0);
31.     ssSetNumPWork(S, 0);
32.     ssSetNumModes(S, 0);
33.     ssSetNumNonsampledZCs(S, 0);
35.     ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);
37.     ssSetOptions(S, 0);
38. }
40. static void mdlInitializeSampleTimes(SimStruct *S)
41. {
42.     ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
43.     ssSetOffsetTime(S, 0, 0.0);
44. }
48. #define MDL_INITIALIZE_CONDITIONS
49. #if defined(MDL_INITIALIZE_CONDITIONS)
50.   static void mdlInitializeConditions(SimStruct *S)
51.   {
52.   }
53. #endif
57. #define MDL_START
58. #if defined(MDL_START)
59.   static void mdlStart(SimStruct *S)
60.   {
61.   }
62. #endif /*  MDL_START */
65. static void mdlOutputs(SimStruct *S, int_T tid)
66. {
67.     const real_T *map = mxGetPr(ssGetSFcnParam(S,0));  //获取定点数据
68.     const int_T *mapSize = mxGetDimensions(ssGetSFcnParam(S,0));  //定点数组维数
69.     const real_T *x = (const real_T*) ssGetInputPortSignal(S,0);  //输入x
70.     real_T       *y = ssGetOutputPortSignal(S,0); //输出y
71.     int_T step = 0;  //输入x在定点数中的位置
72.     int_T i;
73.     real_T yval;
75.     for (i = 0; i < mapSize[0]; i++)
76.     {
77.         if (x[0] >= map[i] && x[0] < map[i + 1])
78.         {
79.             step = i;
80.             break;
81.         }
82.     }
83.    
84.     cubic_getval(&yval, mapSize, map, x[0], step);
85.     y[0] = yval;
86.    
87. }
89. //自然边界的三次样条曲线函数
90. void cubic_getval(real_T* y, const int_T* size, const real_T* map, const real_T x, const int_T step)
91. {
92.     int_T n = size[0];
93.    
94.     //曲线系数
95.     real_T* ai = (real_T*)malloc(sizeof(real_T) * (n-1));
96.     real_T* bi = (real_T*)malloc(sizeof(real_T) * (n-1));
97.     real_T* ci = (real_T*)malloc(sizeof(real_T) * (n-1));
98.     real_T* di = (real_T*)malloc(sizeof(real_T) * (n-1));
99.    
100.     real_T* h = (real_T*)malloc(sizeof(real_T) * (n-1));  //x的??
101.    
102.     /* M矩阵的系数
103.      *[B0, C0, ...
104.      *[A1, B1, C1, ...
105.      *[0,  A2, B2, C2, ...
106.      *[0, ...             An-1, Bn-1]
107.      */
108.     real_T* A = (real_T*)malloc(sizeof(real_T) * (n-2));
109.     real_T* B = (real_T*)malloc(sizeof(real_T) * (n-2));
110.     real_T* C = (real_T*)malloc(sizeof(real_T) * (n-2));
111.     real_T* D = (real_T*)malloc(sizeof(real_T) * (n-2)); //等号右边的常数矩阵
112.     real_T* E = (real_T*)malloc(sizeof(real_T) * (n-2)); //M矩阵
113.    
114.     real_T* M = (real_T*)malloc(sizeof(real_T) * (n));  //包含端点的M矩阵
115.    
116.     int_T i;
117.    
118.     //计算x的步长
119.     for ( i = 0; i < n -1; i++)
120.     {
121.         h[i] = map[i + 1] - map[i];
122.     }
123.    
124.     //指定系数
125.     for( i = 0; i< n - 3; i++)
126.     {
127.         A[i] = h[i]; //忽略A[0]
128.         B[i] = 2 * (h[i] + h[i+1]);
129.         C[i] = h[i+1]; //忽略C(n-1)
130.     }
132.    
133.     //指定常数D
134.     for (i = 0; i<n - 3; i++)
135.     {
136.         D[i] = 6 * ((map[n + i + 2] - map[n + i + 1]) / h[i + 1] - (map[n + i + 1] - map[n + i]) / h[i]);
137.     }
138.    
139.    
140.     //求解三对角矩阵，结果赋值给E
141.     TDMA(E, n-3, A, B, C, D);
142.    
143.     M[0] = 0; //自然边界的首端M为0
144.     M[n-1] = 0;  //自然边界的末端M为0
145.     for(i=1; i<n-1; i++)
146.     {
147.         M[i] = E[i-1]; //其它的M值
148.     }
149.    
150.     //?算三次?条曲?的系数
151.     for( i = 0; i < n-1; i++)
152.     {
153.         ai[i] = map[n + i];
154.         bi[i] = (map[n + i + 1] - map[n + i]) / h[i] - (2 * h[i] * M[i] + h[i] * M[i + 1]) / 6;
155.         ci[i] = M[i] / 2;
156.         di[i] = (M[i + 1] - M[i]) / (6 * h[i]);
157.     }
158.    
159.     *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]);
160.    
161.     free(h);
162.     free(A);
163.     free(B);
164.     free(C);
165.     free(D);
166.     free(E);
167.     free(M);
168.     free(ai);
169.     free(bi);
170.     free(ci);
171.     free(di);
173. }
175. void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D)
176. {
177.     int_T i;
178.     real_T tmp;
180.     //上三角矩阵
181.     C[0] = C[0] / B[0];
182.     D[0] = D[0] / B[0];
184.     for(i = 1; i<n; i++)
185.     {
186.         tmp = (B[i] - A[i] * C[i-1]);
187.         C[i] = C[i] / tmp;
188.         D[i] = (D[i] - A[i] * D[i-1]) / tmp;
189.     }
191.     //直接求出X的最后一个值
192.     X[n-1] = D[n-1];
194.     //逆向迭代， 求出X
195.     for(i = n-2; i>=0; i--)
196.     {
197.         X[i] = D[i] - C[i] * X[i+1];
198.     }
199. }
202. #define MDL_UPDATE
203. #if defined(MDL_UPDATE)
204.   static void mdlUpdate(SimStruct *S, int_T tid)
205.   {
206.   }
207. #endif
209. #define MDL_DERIVATIVES
210. #if defined(MDL_DERIVATIVES)
211.   static void mdlDerivatives(SimStruct *S)
212.   {
213.   }
214. #endif
216. static void mdlTerminate(SimStruct *S)
217. {
218. }
220. #ifdef  MATLAB_MEX_FILE  
221. #include "simulink.c"
222. #else
223. #include "cg_sfun.h"
224. #endif

复制代码

3. 例子
以y=sin(x)为例,  x步长为1，x取值范围是[0,10]。对它使用三次样条插值，插值前后对比如下：