Java学习者论坛

 找回密码
 立即注册

QQ登录

只需一步,快速开始

手机号码,快捷登录

恭喜Java学习者论坛(www.javaxxz.com)已经为数万Java学习者服务超过7年了!积累会员资料超过10000G+
成为本站VIP会员,下载本站10000G+会员资源,购买链接:http://item.taobao.com/item.htm?id=44171550842
成为荣耀会员,分享5TB资料及站长学习指导,购买链接:https://item.taobao.com/item.htm?id=44435180049
技术售后:点击这里给我发消息 资料售后:点击这里给我发消息 ①群:Java学习者群②javaxxz.com ②群:Java学习者群③javaxxz.com 求职招聘群:Java求职与招聘 精英群:Java学习者精英群
JavaEE 49期就业班视频教程Java从菜鸟到大神的学习路线之基础篇Java从菜鸟到大神的学习路线之实战篇Java从菜鸟到大神的学习路线之高级篇

价值两万达内2017年最新Java整套视频

Java开发视频教程下载

大数据开发视频教程

前端开发视频教程下载

安卓开发视频教程下载

Java亿级流量电商系统视频教程

互联网架构师视频教程

年薪50万Spark2.0从入门到精通

年薪50万!人工智能学习路线教程

年薪50万!大数据从入门到精通学习路线年薪50万!机器学习入门到精通视频教程
查看: 100|回复: 0

[默认分类] 三次样条插值(Cubic Spline Interpolation)及代码实现(C语言)

[复制链接]
  • TA的每日心情
    开心
    2018-8-25 14:10
  • 签到天数: 222 天

    [LV.7]常住居民III

    发表于 2018-7-13 16:36:19 | 显示全部楼层 |阅读模式
    样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法,三次样条又是其中用的较为广泛的一种。本篇介绍力求用容易理解的方式,介绍一下三次样条插值的原理,并附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"  //方便使用变量定义数组大小
    5. static void mdlInitializeSizes(SimStruct *S)
    6. {
    7.     /*参数只有一个,是n乘2的定点数组[xi, yi]:
    8.      * [ x1,y1;
    9.      *   x2, y2;
    10.      *   ..., ...;
    11.      *   xn, yn;
    12.     */
    13.     ssSetNumSFcnParams(S, 1);
    14.     if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
    15.     ssSetNumContStates(S, 0);
    16.     ssSetNumDiscStates(S, 0);
    17.     if (!ssSetNumInputPorts(S, 1)) return;  //输入是x
    18.     ssSetInputPortWidth(S, 0, 1);
    19.     ssSetInputPortRequiredContiguous(S, 0, true);
    20.     ssSetInputPortDirectFeedThrough(S, 0, 1);
    21.     if (!ssSetNumOutputPorts(S, 1)) return;  //输出是S(x)
    22.     ssSetOutputPortWidth(S, 0, 1);
    23.     ssSetNumSampleTimes(S, 1);
    24.     ssSetNumRWork(S, 0);
    25.     ssSetNumIWork(S, 0);
    26.     ssSetNumPWork(S, 0);
    27.     ssSetNumModes(S, 0);
    28.     ssSetNumNonsampledZCs(S, 0);
    29.     ssSetSimStateCompliance(S, USE_DEFAULT_SIM_STATE);
    30.     ssSetOptions(S, 0);
    31. }
    32. static void mdlInitializeSampleTimes(SimStruct *S)
    33. {
    34.     ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
    35.     ssSetOffsetTime(S, 0, 0.0);
    36. }
    37. #define MDL_INITIALIZE_CONDITIONS
    38. #if defined(MDL_INITIALIZE_CONDITIONS)
    39.   static void mdlInitializeConditions(SimStruct *S)
    40.   {
    41.   }
    42. #endif
    43. #define MDL_START
    44. #if defined(MDL_START)
    45.   static void mdlStart(SimStruct *S)
    46.   {
    47.   }
    48. #endif /*  MDL_START */
    49. static void mdlOutputs(SimStruct *S, int_T tid)
    50. {
    51.     const real_T *map = mxGetPr(ssGetSFcnParam(S,0));  //获取定点数据
    52.     const int_T *mapSize = mxGetDimensions(ssGetSFcnParam(S,0));  //定点数组维数
    53.     const real_T *x = (const real_T*) ssGetInputPortSignal(S,0);  //输入x
    54.     real_T       *y = ssGetOutputPortSignal(S,0); //输出y
    55.     int_T step = 0;  //输入x在定点数中的位置
    56.     int_T i;
    57.     real_T yval;
    58.     for (i = 0; i < mapSize[0]; i++)
    59.     {
    60.         if (x[0] >= map[i] && x[0] < map[i + 1])
    61.         {
    62.             step = i;
    63.             break;
    64.         }
    65.     }
    66.    
    67.     cubic_getval(&yval, mapSize, map, x[0], step);
    68.     y[0] = yval;
    69.    
    70. }
    71. //自然边界的三次样条曲线函数
    72. void cubic_getval(real_T* y, const int_T* size, const real_T* map, const real_T x, const int_T step)
    73. {
    74.     int_T n = size[0];
    75.    
    76.     //曲线系数
    77.     real_T* ai = (real_T*)malloc(sizeof(real_T) * (n-1));
    78.     real_T* bi = (real_T*)malloc(sizeof(real_T) * (n-1));
    79.     real_T* ci = (real_T*)malloc(sizeof(real_T) * (n-1));
    80.     real_T* di = (real_T*)malloc(sizeof(real_T) * (n-1));
    81.    
    82.     real_T* h = (real_T*)malloc(sizeof(real_T) * (n-1));  //x的??
    83.    
    84.     /* M矩阵的系数
    85.      *[B0, C0, ...
    86.      *[A1, B1, C1, ...
    87.      *[0,  A2, B2, C2, ...
    88.      *[0, ...             An-1, Bn-1]
    89.      */
    90.     real_T* A = (real_T*)malloc(sizeof(real_T) * (n-2));
    91.     real_T* B = (real_T*)malloc(sizeof(real_T) * (n-2));
    92.     real_T* C = (real_T*)malloc(sizeof(real_T) * (n-2));
    93.     real_T* D = (real_T*)malloc(sizeof(real_T) * (n-2)); //等号右边的常数矩阵
    94.     real_T* E = (real_T*)malloc(sizeof(real_T) * (n-2)); //M矩阵
    95.    
    96.     real_T* M = (real_T*)malloc(sizeof(real_T) * (n));  //包含端点的M矩阵
    97.    
    98.     int_T i;
    99.    
    100.     //计算x的步长
    101.     for ( i = 0; i < n -1; i++)
    102.     {
    103.         h[i] = map[i + 1] - map[i];
    104.     }
    105.    
    106.     //指定系数
    107.     for( i = 0; i< n - 3; i++)
    108.     {
    109.         A[i] = h[i]; //忽略A[0]
    110.         B[i] = 2 * (h[i] + h[i+1]);
    111.         C[i] = h[i+1]; //忽略C(n-1)
    112.     }
    113.    
    114.     //指定常数D
    115.     for (i = 0; i<n - 3; i++)
    116.     {
    117.         D[i] = 6 * ((map[n + i + 2] - map[n + i + 1]) / h[i + 1] - (map[n + i + 1] - map[n + i]) / h[i]);
    118.     }
    119.    
    120.    
    121.     //求解三对角矩阵,结果赋值给E
    122.     TDMA(E, n-3, A, B, C, D);
    123.    
    124.     M[0] = 0; //自然边界的首端M为0
    125.     M[n-1] = 0;  //自然边界的末端M为0
    126.     for(i=1; i<n-1; i++)
    127.     {
    128.         M[i] = E[i-1]; //其它的M值
    129.     }
    130.    
    131.     //?算三次?条曲?的系数
    132.     for( i = 0; i < n-1; i++)
    133.     {
    134.         ai[i] = map[n + i];
    135.         bi[i] = (map[n + i + 1] - map[n + i]) / h[i] - (2 * h[i] * M[i] + h[i] * M[i + 1]) / 6;
    136.         ci[i] = M[i] / 2;
    137.         di[i] = (M[i + 1] - M[i]) / (6 * h[i]);
    138.     }
    139.    
    140.     *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]);
    141.    
    142.     free(h);
    143.     free(A);
    144.     free(B);
    145.     free(C);
    146.     free(D);
    147.     free(E);
    148.     free(M);
    149.     free(ai);
    150.     free(bi);
    151.     free(ci);
    152.     free(di);
    153. }
    154. void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D)
    155. {
    156.     int_T i;
    157.     real_T tmp;
    158.     //上三角矩阵
    159.     C[0] = C[0] / B[0];
    160.     D[0] = D[0] / B[0];
    161.     for(i = 1; i<n; i++)
    162.     {
    163.         tmp = (B[i] - A[i] * C[i-1]);
    164.         C[i] = C[i] / tmp;
    165.         D[i] = (D[i] - A[i] * D[i-1]) / tmp;
    166.     }
    167.     //直接求出X的最后一个值
    168.     X[n-1] = D[n-1];
    169.     //逆向迭代, 求出X
    170.     for(i = n-2; i>=0; i--)
    171.     {
    172.         X[i] = D[i] - C[i] * X[i+1];
    173.     }
    174. }
    175. #define MDL_UPDATE
    176. #if defined(MDL_UPDATE)
    177.   static void mdlUpdate(SimStruct *S, int_T tid)
    178.   {
    179.   }
    180. #endif
    181. #define MDL_DERIVATIVES
    182. #if defined(MDL_DERIVATIVES)
    183.   static void mdlDerivatives(SimStruct *S)
    184.   {
    185.   }
    186. #endif
    187. static void mdlTerminate(SimStruct *S)
    188. {
    189. }
    190. #ifdef  MATLAB_MEX_FILE  
    191. #include "simulink.c"
    192. #else
    193. #include "cg_sfun.h"
    194. #endif
    复制代码



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

    回复

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    .

    QQ|手机版|Java学习者论坛

    GMT+8, 2018-11-19 12:56 , Processed in 0.389505 second(s), 26 queries .

    Powered by Discuz! X3.4

    © 2001-2017 Comsenz Inc.

    快速回复 返回顶部 返回列表