Matalab插值详解和源码

转载：Matalab插值详解和源码 - 知乎 (zhihu.com)

插值法

插值法又称“内插法”，是利用函数f (x)在某区间中已知的若干点的函数值，作出适当的特定函数，在区间的其他点上用这特定函数的值作为函数f (x)的近似值，这种方法称为插值法。如果这特定函数是多项式，就称它为插值多项式。常见分段线性插值法和样条差值，样条插值误差更小。

1 线性插值法

线性插值法是指使用连接两个已知量的直线来确定在这两个已知量之间的一个未知量的值的方法。假设我们已知坐标(x0,y0)与(x1,y1),要得到[x0,x1]区间内某一位置x在直线上的y值。根据图中所示，我们得到两点式直线方程：

假设方程两边的值为α,则：

则：y = y0 + α(y1 − y0)(即已知x就可得a值然后再得y值。

实际上，即使x不在x0到x1之间并且α也不是介于0到1之间，这个公式也是成立的。在这种情况下，这种方法叫作线性外插—参见外插值。

已知y求x的过程与以上过程相同。

2 样条插值法

常用三次样条插值法，具体概念比较复杂，可以自行百度。

3 拉格朗日差值法

略。

MATLAB中的使用-interp1函数：

MATLAB中的插值函数为interp1，其调用格式为： yi=interp1(x,y,xi,'method') 其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量，
'method'表示采用的插值方法，MATLAB提供的插值方法有几种：
'nearest'是最邻近插值， 'linear'线性插值； 'spline'三次样条插值； 'pchip'立方插值．缺省时表示线性插值
注意：所有的插值方法都要求x是单调的，并且xi不能够超过x的范围。

x = 0:2*pi;    
y = sin(x);    
xx = 0:0.5:2*pi;    
  
% interp1对sin函数进行分段线性插值，调用interp1的时候，默认的是分段线性插值    
y1 = interp1(x,y,xx,'linear');    
subplot(2,2,1);  
plot(x,y,'o',xx,y1,'r')    
title('分段线性插值')    
    
% 临近插值    
y2 = interp1(x,y,xx,'nearest');    
subplot(2,2,2);  
plot(x,y,'o',xx,y2,'r');    
title('临近插值')    
    
%球面线性插值    
y3 = interp1(x,y,xx,'spline');    
subplot(2,2,3);  
plot(x,y,'o',xx,y3,'r')    
title('球面插值')    
    
%三次多项式插值法    
y4 = interp1(x,y,xx,'pchip');    
subplot(2,2,4);  
plot(x,y,'o',xx,y4,'r');    
title('三次多项式插值')

例：环境温度

例如：在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为
12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，
推测中午1点（即13点）时的温度．

x = 0:2:24;  
y = [12   9   9   10   18  24   28   27   25   20  18  15  13];  
a = 13;  
y1 = interp1(x,y,a,'spline')  
% 结果为：  27.8725  

% 若要得到一天24小时的温度曲线，则：  
xi = 0:1/3600:24;  
% 插值点可以是向量，则返回的也就是对应的向量  
yi = interp1(x,y,xi, 'spline');  
plot(x,y,'o' ,xi,yi);

一天24小时温度曲线

语法形式	说明
y=interp1(x,Y,xi)	由已知点集(x,Y)插值计算xi上的函数值
y=interp1(Y,xi)	相当于x=1:length(Y)的interp(x,Y,xi)
y=interp1(x,Y,xi,method)	用指定插值方法计算插值点xi上的函数值
y=interp1(x,Y,xi,method,’pp’)	用指定方法插值，但返回结果为分段多项式

Method	方法描述
‘nearest’	最邻近插值:插值点处函数值与插值点最邻近的已知点函数值相等
‘liner’	分段线性插值：插值点处函数值由连接其最邻近的两侧点的线性函数预测。Matlab中interp1的默认方法。
‘spline’	样条插值：默认为三次样条插值。可用spline函数替代即y=spline(x0,y0,x)
‘pchip’	三次Hermite多项式插值，可用pchip函数替代
‘cubic’	同’pchip’，三次Hermite多项式插值

Csape函数

csape函数的用法如下：
pp = csape(x,y,conds,valconds)
其中(x,y)为数据向量，conds表示变界类型， valconds表示边界值。
边界类型(conds)可为：
'complete'，给定边界一阶导数.
'not-a-knot'，非扭结条件,不用给边界值.
'periodic'，周期性边界条件,不用给边界值.
'second'，给定边界二阶导数.
'variational'，自然样条(边界二阶导数为0)
边界类型(valconds)可为：
'complete'，给定边界一阶导数.
'not-a-knot'，非扭结条件,不用给边界值.
'periodic'，周期性边界条件,不用给边界值.
'second'，给定边界二阶导数.
'variational'，自然样条(边界二阶导数为0)

例：机床加工

待加工零件的外形根据工艺要求由一组数据(x,y)给出，用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x,y)坐标。

依据下表给出的(x,y)数据求x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求出x=0处的曲线斜率和13<=x<=15范围内y的最小值。

x 0 3 5 7 9 11 12 13 14 15
y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6

x0=[0  3  5  7  9  11  12  13  14  15];
y0=[0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6];
x=0:0.1:15;
y1=interp1(x0,y0,x);%分段线性插值 得到x=0:0.1:15相应的函数值向量y1
y2=interp1(x0,y0,x,'spline');%三次样条插值 得到x=0:0.1:15相应的函数值向量y2
pp1=csape(x0,y0);%csape函数三次样条插值，默认边界条件为拉格朗日边界条件 
y3=fnval(pp1,x);%得到x=0:0.1:15相应的函数值向量y3
pp2=csape(x0,y0,'second');%设置边界条件为二阶导数，默认为[0,0]
y4=fnval(pp2,x);%得到x=0:0.1:15相应的函数值向量y4
[x',y1',y2',y3',y4'];%得到一个151*4的矩阵可以观察y1,y2,y3,y4随x的变化

subplot(1,3,1)%将当前图窗分为1*3网络并在1处创建坐标区（即当前选定1区绘图）
plot(x0,y0,'+',x,y1);%(x0,y0)所代表的点用+表示
title('分段线性插值');

subplot(1,3,2)%选定2区绘图
plot(x0,y0,'+',x,y2);
title('三次样条插值-interp1函数')

subplot(1,3,3)%选定3区绘图
plot(x0,y0,'+',x,y3);
title('三次样条插值-scape函数')

%求x=0处斜率
dx=diff(x);%dx即∆x
dy=diff(y3);%dy即∆y
dy_dx=dy./dx;%(./ 对应元素相除)
dy_dx0=dy_dx(1);%x=0处的导数就是dy./dx的第一个元素:dy0/dx0

%求13<=x<=15范围内y的最小值
ytemp=y3(131:151);%取出13<=x<=15区间内对应的y3值（因为y3误差更小，数据分析更准确）
ymin=min(ytemp);
index=find(y3==ymin);%y3中最小值对应的序列
xmin=x(index);%并由序列找到对应的x值
[xmin,ymin]%输出

运行程序：