C语言和MATLAB一维线性插值程序实现对比

在实际的数据处理或者产品开发时，经常会使用到数据的插值，其中一维线性插值是使用比较多的，数学原理比较简单，编程实现也方便，今天主要是以实际的例子介绍一下C语言的一维线性插值实现与MATLAB的一维线性插值函数interp1的对比。

一次函数的5种形式

斜截式：y=kx+b，k是斜率,y是一次函数与y轴的交点 
点斜式：y-y0=k(x-x0) ，k是斜率，(x0,y0）是已知的点   
两点式：(y-y1)(x-x2)=(y-y2)(x-x1),(k = (y2-y1)/(x2-x1))
(x1,y1),(x2,y2)的已知的点
截距式：x/a+y/b=1（其中a,b分别为该直线在x轴和y轴上的截距）
一般式：Ax+By+C=0.

1.C语言的一维线性插值程序

已知两组线性关系的数据为

求x = 4时，对应的y值是多少

主程序mian.c

#include <stdio.h>
#include <math.h>
#include "method.h"
double X[4] = {2,5,8,20};
double Y[4] = {2*RATIO, 3*RATIO, 6*RATIO, 7*RATIO};
int main()
{
    double x = 4;
    double y = interp1(x, X, Y, sizeof(X)/sizeof(X[0]));
    printf("(x = %lf,y = %lf) \n",x,y);/*Factor:0.01*/
    return 0;
}

method.c

#include "method.h"
double interp1(double x, const double *Xtable, const double *Ytable, u32 Len)
{
    u32 i = 0;
    /*数组长度*/
  u32 j = Len - 1;
    double t = 0;
    /*超出x的下限，此时反回YTable[0]*/
    if (x <= Xtable[0])
    {
      return Ytable[0];
  }
    /*超出x的上限，此时反回YTable[end]*/    
    if (x >= Xtable[Len - 1])
    {
      return Ytable[Len - 1];
  }
    
  /*寻找数据表中最接近的数据的两个端点 x(i) x(x) x(j)*/    
    while (i < j - 1)
    {
        t = (i + j) / 2;
        if (x < Xtable[(u32)t])
        {
          j = (u32)t;  
    } 
        else
        {
          i = (u32)t;  
    }     
    }
    t = (x - Xtable[i]) / (Xtable[j] - Xtable[i]);
    return Ytable[i] + t * (Ytable[j] - Ytable[i]);
    /* (x0,y0) = (xtable[i],ytable[i])    (x1,y1) = (xtable[i],ytable[i]) (x2,y2) = (xtable[j],ytable[j])
    (x-xtable[i])*(Ytable[j] - Ytable[i]) / (Xtable[j] - Xtable[i])= K*(x-x0)+y0  
  y-y0 = k*(x - x0)*/
    
}

method.h

#ifndef _METHOD_H
#define _METHOD_H


typedef  unsigned char u8;
typedef  unsigned short u16;
typedef  unsigned int u32;
typedef  unsigned long int u64;


#define RATIO  100


extern double interp1(double x, const double *Xtable, const double *Ytable, u32 Len);


#endif 

运行结果

2.MATLAB的interp1函数

一维插值函数为interp1

调用格式：
y = interp1(X,Y,X1,method)
该式可以根据X，Y的值来计算函数在X1处的值。
其中X，Y是两个等长的已知向量，分别表示采样点和采样值。
X1是一个向量或标量，表示要插值的点。
method参数表示用于插值的方法，常用的取值由以下几种方法：
(1) linear: 线形插值，默认方法。
将与插值点靠近的两个数据点用直线连接，然后在直线上选取对应插值点
的数据。
(2) nearest: 最近点插值。
选择最近样本点的值作为插值数据。
(3) pchip: 分段3次埃尔米特插值。
采用分段三次多项式，除满足插值条件，还需满足在若干节点处相邻段插值
函数的一阶导数相等，使得曲线光滑的同时，还具有保形性。
(4) spline: 3次样条插值。
每一个分段的内构造一个三次多项式，使其插值函数除满足插值条件外，
还要求在各节点处具有连续的一阶和二阶导数。

以上四种方法的区别：

线形插值和最近点插值方法比较简单。其中线形插值方法的计算量与样本点n 无关。n越大，误差越小。

3次埃尔米特插值和3次样条插值都能保证曲线的光滑性。相比较而言，3次埃尔米特插值具有保形性，而3次样条插值要求其二阶导数也连续，所以插值函数的性态更好。

实例1

程序

clc;
clear all;
close all;
x0 = [2 5 8 20];
y0 = [200 300 600 700];
xi = 4;
y = interp1(x0,y0,xi,'linear');
fprintf("x = %f, y = %f\n",xi,y)

运行结果

x = 4.000000, y = 266.666667

实例2

程序

clc;
clear all;
close all;
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];
x1 = 0:0.1:15;
y1 = interp1(x,y,x1,'spline');
y2 = interp1(x,y,x1,'linear');
y3 = interp1(x,y,x1,'nearest');
y4 = interp1(x,y,x1,'pchip');
figure;
plot(x,y,'r*');
hold on;
plot(x1,y1,'b');
hold on;
plot(x1,y2,'g')
hold on;
plot(x1,y3,'black')
hold on;
plot(x1,y4,'r')
xlabel('x');
ylabel('y');
grid on;
legend('原始数据','spline 3次样条插值','linear 线形插值',...
    'nearest 最近点插值','pchip 分段3次埃尔米特插值','location','southeast')

运行结果

3.参考内容

[1] 今日头条作者酒足饭饱抡大锤的文章《用c语言实现matlab的一维插值函数interp1》，文章链接为：https://www.toutiao.com/article/7197987953841734159/?wid=1729415202883

本文内容来源于网络，仅供参考学习，如内容、图片有任何版权问题，请联系处理，24小时内删除。

作 者 | 郭志龙

编 辑 | 郭志龙
校 对 | 郭志龙