1.4 数据插值
我们接触的数据有时是存在缺失的，那么根据前后的数据，如何较精确地知道缺失的数据是什么？这就需要用到数据插值的方法。数据插值包括一维数据插值和二维数据插值，不同的插值方法得到的图形数据精度是不一样的。
1.4.1 一维数据插值
MATLAB提供的函数interp1()，用于进行一维插值，所谓一维插值就是根据二维平面散点图进行平滑插值，达到拟合和预测的目的。
函数interp1()可以有3种调用形式。
（1）yi=interp1(x, y, xi)：根据数据[x, y]，用分段线性插值算法得到xi各点上的函数近似值yi，yi为尽可能逼近的最小误差对应的因变量值。当y是向量时，则对y向量插值，得到结果yi是与xi同样大小的向量；当y是矩阵时，则对y的逐列向量插值，得到结果yi是一矩阵，它的列数与y的列数相同，行数与xi的大小相同。
（2）yi=interp1(y, yi)：格式调用方法与yi=interp1(x, y, xi)相同，只是插值节点不同，此时用节点序号x=1:n，n=size(y)。
（3）yi=interp1(x, y, xi, method)：函数调用与yi=interp1(x, y, xi)和yi=interp1(y, yi)相同，只是需要指定输入参数method的具体算法。在MATLAB中，method有6种形式如下：
>> help interp1
Vq = interp1(X,V,Xq,METHOD) specifies alternate methods.
X：x
V: y
Xq: xi
METHOD: Method,插值方法
'nearest' - 最近邻插值
'linear' - 线性插值
'spline' - 样条曲线插值
'pchip' - 分段立方插值
'cubic' - 立方插值
'v5cubic' - 如果X不是等距节点，不进行外插，采用样条插值
其中，最常用的是第5种插值方式，当然用户可以选择不同的插值方法去实现插值操作。
1）nearest：最近邻插值，即在插值节点的近邻区域内取函数值，该插值函数为一个阶梯函数，即：
具体的MATLAB最近邻插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10]; % x
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7]; % y
figure('color',[1,1,1]) % 先建一个图形窗口
hold on % 图像保持画图句柄
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10; % x
yi = interp1(x,y,xi,'nearest'); % 最近邻插值
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y');
legend('原始数据','nearest插值')
toc % 计时结束
运行程序输出结果如图1-25所示。
图1-25 nearest插值
2）linear：分段线性插值，当不指定插值方法时，默认采用线性插值方法。
分段线性插值是在每个小区间
上采用简单的线性插值。在区间
上的子插值多项式为：
分段线性插值方法较为常用，在实际计算中，处理速度较快，但是对海量数据本身而言，以及非线性问题，处理误差较大，线性插值方法获得的曲线不是平滑的，因此，应根据实际情况要求选择。
表1-6所示为待插值数据。
表1-6 源数据
具体的MATLAB线性插值程序如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7]; % y值
figure('color',[1,1,1]) % 先建图形窗口
hold on
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10;
yi = interp1(x,y,xi); % 默认情况下为线性插值
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y'); % xy轴标记
legend('原始数据','linear插值')
toc % 计时结束
运行程序输出图形如图1-26所示。
图1-26 linear插值
3）spline：三次样条插值，即在每个分段（子区间）内构造一个三次多项式，使其插值函数满足差值条件外，还要求在各个节点处具有光滑的条件（导数存在）。
三次样条函数s（x）在每个子区间
上可由4个系数唯一确定。因此，s（x）在［a，b］上有4n个待定系数。
由于
，则有：
为了确定s（x），通常还需要补充边界条件，常用的边界条件分为三类。
第一种边界条件：给定两边界节点处的一阶导数
，并要求s（x）满足
。
第二种边界条件：给定两边界节点处的二阶导数
，并要求s（x）满足
。
特别地，若
，则所得的样条称为自然样条。
第三种边界条件：被插函数f（x）是以
为周期的周期函数，要求s（x）满足
，
，
。
具体的三次样条插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7]; % y轴
figure('color',[1,1,1])
hold on
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10;
yi = interp1(x,y,xi,'spline'); % 样条插值
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y');
legend('原始数据','spline插值')
toc % 计时结束
运行程序输出图形如图1-27所示。
图1-27 spline插值
4）pchip：分段立方插值。分段立方插值的思想就是将插值的数据分段，针对每一段采用立方插值的方法进行拟合逼近，具体的MATLAB代码如下：
x = [1:10];
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7];
figure('color',[1,1,1])
hold on
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10;
yi = interp1(x,y,xi,'pchip'); % 分段立方插值
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y');
legend('原始数据','pchip插值')
运行程序输出图形如图1-28所示。
图1-28 pchip插值
5）cubic：立方插值。cubic插值法插值精度较高，插值曲线较平滑。
具体的MALTAB cubic插值程序如下：
tic % 运算计时
x = [1:10];
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7];
figure('color',[1,1,1])
hold on
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10;
yi = interp1(x,y,xi,'cubic'); % 双立方插值
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y');
legend('原始数据','cubic插值')
toc % 计时结束
运行程序得到如图1-29所示图形。
图1-29 cubic插值
6）v5cubic：MATLABv5版本cubic插值方法。v5cubic是一种旧MATLAB版本下的cubic插值方法，类似于三次样条插值方法，该插值方法要求插值点为均匀间隔的插值点。这种v5cubic方法逐渐被cubic取代。
具体的MALTAB v5cubic插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [1, 1, 7, 4, 0, 6, 3, 0, 8, 7];
figure('color',[1,1,1])
hold on
plot(x,y,'o--','linewidth',2) % 原数图形
xi = 1:0.1:10;
yi = interp1(x,y,xi,'v5cubic');
hold on % 图形保持画图
grid on % 网格化
plot(xi,yi,'r*--','linewidth',2) % 原数图形
xlabel('x'); ylabel('y');
legend('原始数据','v5cubic插值')
toc % 计时结束
运行程序得到如图1-30所示图形。
图1-30 v5cubic插值
1.4.2 二维数据插值
二维数据插值是构造一个二元插值函数
去近似
，即曲面插值。与interp1类似，指令griddata可以根据数据表[x, y, z]，用各种不同的算法计算[xi, yi]各点上的函数近似值zi。
（1）zi = griddata(x, y, z, xi, yi)：本指令格式根据数据表[x, y, z]，用双线性差值算法计算坐标平面x-y上[xi, yi]各点的二元函数近似值zi，这里x可以是一个行向量，它与矩阵z的各列向量相对应；y可以为一个列向量，它与z的各行向量相对应。对于[xi, yi]与zi间的对应关系，则和[x, y]与z的关系相同。
（2）zi =griddata (x, y, z, xi, yi, method)：本指令格式的调用方法与griddata(x, y, z, xi, yi)格式相同，只是规定在格式中指定具体的算法。
在MATLAB中，method有如下4种形式。
Vq = griddata(...,METHOD) specifies alternate methods． The default
is linear interpolation． Available methods are:
METHOD: Method,插值方法
'nearest' - 最近邻插值
'linear' - 线性插值
'cubic' - 立方插值
表1-7所示为待插值数据。
表1-7 源数据
绘制表1-7所示数据的三维曲面，程序如下：
% Designed by Yu Shengwei From SWJTU University
% 2014年12月29日
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [0, 0, 1, 1, 3, 4, 6, 7, 7, 8];
z = [8, 2, 8, 0, 3, 7, 8, 2, 4, 1];
[X,Y] = meshgrid(x,y); % 栅格
Z = z'*ones(1,length(x));
mesh(X,Y,Z)
toc % 计时结束
运行程序输出图形如图1-31所示。
图1-31 原始数据图
1）采用最近邻插值nearest。
具体的MALTAB nearest插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [0, 0, 1, 1, 3, 4, 6, 7, 7, 8]; % y值
z = [8, 2, 8, 0, 3, 7, 8, 2, 4, 1]; % z值
Z = z'*ones(1,length(x));
figure('color',[1,1,1]) % 先建一个图形窗口
xi = min(x):(max(x)-min(x))/50:max(x); % 插值数据
yi = min(y):(max(y)-min(y))/50:max(y); % 插值数据
yi = yi'; % 转置
Zi = griddata(x,y,Z,xi,yi,'nearest'); % 插值
mesh(xi,yi,Zi) % 画图
grid on % 网格化
xlabel('x'); ylabel('y'); zlabel('z') % xyz轴标记
title('nearest插值')
toc % 计时结束
运行程序得到如图1-32所示图形。
图1-32 nearest插值
2）采用线性插值linear。
具体的MALTAB linear插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [0, 0, 1, 1, 3, 4, 6, 7, 7, 8]; % y值
z = [8, 2, 8, 0, 3, 7, 8, 2, 4, 1]; % z值
Z = z'*ones(1,length(x));
figure('color',[1,1,1]) % 先建一个图形窗口
xi = min(x):(max(x)-min(x))/50:max(x); % 插值数据
yi = min(y):(max(y)-min(y))/50:max(y); % 插值数据
yi = yi'; % 转置
Zi = griddata(x,y,Z,xi,yi,'linear'); % 插值
mesh(xi,yi,Zi) % 画图
grid on % 网格化
xlabel('x'); ylabel('y'); zlabel('z')
title('linear插值')
toc % 计时结束
运行程序得到如图1-33所示图形。
图1-33 linear插值
3）采用cubic插值。
具体的MALTAB cubic插值程序如下：
clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
tic % 运算计时
x = [1:10];
y = [0, 0, 1, 1, 3, 4, 6, 7, 7, 8];
z = [8, 2, 8, 0, 3, 7, 8, 2, 4, 1];
Z = z'*ones(1,length(x));
figure('color',[1,1,1])
xi = min(x):(max(x)-min(x))/50:max(x); % 插值数据
yi = min(y):(max(y)-min(y))/50:max(y); % 插值数据
yi = yi'; % 转置
Zi = griddata(x,y,Z,xi,yi,'cubic'); % 插值
mesh(xi,yi,Zi) % 画图
grid on % 网格化
xlabel('x'); ylabel('y'); zlabel('z')
title(' cubic插值')
toc % 计时结束
运行程序得到如图1-34所示图形。
图1-34 cubic插值