均匀三次B样条曲线插值Matlab程序
&%============均匀三次B样条曲线插值===========
%定义变量:
%X:原始资料,d:控制顶点
%n:数据条数,k:B样条的次数
X=load('data.txt');
n=length(X); %得数据维数;
%A:方程系数-----------------------------
A=zeros(n+2);
A(1,1)=1;A(1,2)=-2;A(1,3)=1;
A(n+2,n)=1;A(n+2,n+1)=-2;A(n+2,n+2)=1;
for i=2:(n+1)
& & A(i,i-1)=1;
& & A(i,i)=4;
& & A(i,i+1)=1;
%e:方程右边.得到的控制点首尾与插值数据相同.且与控制首末连线相切
e(n+2,:)=[0,0];
【声明】:黑吧安全网()登载此文出于传递更多信息之目的，并不代表本站赞同其观点和对其真实性负责，仅适于网络安全技术爱好者学习研究使用，学习中请遵循国家相关法律法规。如有问题请联系我们，联系邮箱，我们会在最短的时间内进行处理。
上一篇：【】【】样条插值是一种工业设计中常用的、得到平滑曲线的一种插值方法，三次样条又是其中用的较为广泛的一种。本篇介绍力求用容易理解的方式，介绍一下三次样条插值的原理，并附C语言的实现代码。
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个未知系数。
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)推出
(i = 0, 1, &, n-2)推出
由此可得：
(i = 0, 1, &, n-2)推出
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
#define S_FUNCTION_LEVEL 2
#include "simstruc.h"
#include "malloc.h"
//方便使用变量定义数组大小
static void mdlInitializeSizes(SimStruct *S)
/*参数只有一个，是n乘2的定点数组[xi, yi]:
* [ x1，y1;
ssSetNumSFcnParams(S, 1);
if (ssGetNumSFcnParams(S) != ssGetSFcnParamsCount(S)) return;
ssSetNumContStates(S, 0);
ssSetNumDiscStates(S, 0);
if (!ssSetNumInputPorts(S, 1)) return;
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)
#define MDL_START
#if defined(MDL_START)
static void mdlStart(SimStruct *S)
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);
*y = ssGetOutputPortSignal(S,0); //输出y
int_T step = 0;
//输入x在定点数中的位置
for (i = 0; i & mapSize[0]; i++)
if (x[0] &= map[i] && x[0] & map[i + 1])
cubic_getval(&yval, mapSize, map, x[0], step);
//自然边界的三次样条曲线函数
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));
/* M矩阵的系数
*[B0, C0, ...
*[A1, B1, C1, ...
A2, B2, C2, ...
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矩阵
//计算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]);
void TDMA(real_T* X, const int_T n, real_T* A, real_T* B, real_T* C, real_T* D)
//上三角矩阵
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] /
D[i] = (D[i] - A[i] * D[i-1]) /
//直接求出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)
#define MDL_DERIVATIVES
#if defined(MDL_DERIVATIVES)
static void mdlDerivatives(SimStruct *S)
static void mdlTerminate(SimStruct *S)
MATLAB_MEX_FILE
#include "simulink.c"
#include "cg_sfun.h"
以y=sin(x)为例, &x步长为1，x取值范围是[0,10]。对它使用三次样条插值，插值前后对比如下：
阅读(...) 评论() 上传我的文档
 下载
 收藏
如有需要或疑问，请联系站内信或我的QQ：
 下载此文档
正在努力加载中...
MATLAB实现三次样条插值函数
下载积分：400
内容提示：论文
MATLAB实现三次样条插值函数
文档格式：DOC|
浏览次数：213|
上传日期： 14:21:04|
文档星级：
该用户还上传了这些文档
MATLAB实现三次样条插值函数
官方公共微信求问，MATLAB来做三次样条插值，如何得到插值的函数表达式_百度知道
求问，MATLAB来做三次样条插值，如何得到插值的函数表达式
2.98 0.38]:1.0].81 0; y=[0.64 0.92 0x=[0.2:0
提问者采纳
9200-0，三次样条插值每段是三次多项式.7292 -0; x=[0.2..8 0.1042 -0; pp.3 0.].38].64 0.6875 -0;&&gt,y)&y=[0.98 0; pp=spline(x.2. 0.;&gt.9800-0;&gt.0458 0.coefsans =-0&gt.8100-0.6400 返回的是三次样条插值函数每段的系数.92 0
你好，请问如何解释：spline(x,y,0.3)ans =
0.9559但是 -0.^3 -0.^2 -0.+ 0.5而非0.9559
你好，得出来的分段函数是这样的：-0.1042*(t-0.2)^ 3-0.5625*(t-0.2)^2 -0.1833*(t-0.2) + 0.9800-0.1042*(t-0.4)^ 3 -0.6250*(t-0.4)^2
-0.4208*(t-0.4) + 0.9200-0.7292*(t-0.6)^ 3 -0.6875*(t-0.6)^2 -0.6833*(t-0.6)+ 0.8100-0.7292*(t-0.8)^ 3 -1.1250*(t-0.8)^2
-1.0458*(t-0.8)+ 0.6400建议可以看一下数值分析方面的书籍。
提问者评价
万分感谢！！！请问有推荐的MATLAB在数值分析上应用的书吗？
其他类似问题
figure(3)plot(x;plot(x.2,y)clcclearx=[0.92 0.0],x).98 0.38]:0,'Y=ppval(P,Y.81 0.2,y)hold onP=spline(x.64 0:1;y=[0;k&#39
为您推荐：
三次样条插值的相关知识
其他1条回答
1.81 0:0.2.98 0.0].64 0; y=[0=[0.2.38].92 0
等待您来回答
下载知道APP
随时随地咨询
出门在外也不愁matlab中三次样条插值能解决什么问题?
特锐德QJ80
插值的时候,根据你选择插值的对象不同,有时需要选择不同的方法,期望达到最合适的插值效果.你在interp方程里面的method可以看到.spline fitting的最大好处就是可以用较低的阶数（比较简单的多项式）来对比较复杂的数据进行拟合,这样的好处就是可以避免使用高阶多项式的时候引入一些不必要的误差（你可以理解为噪音,比如傅氏级数中的高阶项）.所以,如果你期望插值的对象有可能比较复杂的数量关系,并非简单的线性或者二次曲面,spline fitting可以带来相对较好的插值结果,同时计算得也比较快.
为您推荐：
其他类似问题
扫描下载二维码