最近在做心电相关的项目，由于单片机的处理能力有限，在接收心电信号之后需要对数据进行压缩(其实是取一些特征点)，然后后期再进行显示。但是在手持ARM上进行显示的时候，通过这些残缺的数据绘出心电图形是很困难的，这就要进行插值处理，所以进行了一些插值算法相关的研究，常用的插值算法是拉格郎日插值和牛顿插值算法，切比雪夫插值算法。
一. 拉格朗日插值
拉格郎日的原理，下面这篇帖子写的比较好，感兴趣可以看一下，对于应用型达人来讲，不必对原理过于纠结。
http://www.360doc.com/content/11/1006/13/4539198_153780761.shtml
复杂的公式：
其中每个
为拉格朗日基本多项式(或称插值基函数)，其表达式为：
拉格朗日基本多项式
的特点是在
上取值为1，在其它的点
上取值为0。
拉格郎日插值的本质就是通过一些离散的点预测一个函数，然后得到未知点所对应的函数值。很容易想到的就是n阶多项式。
f(x) = a0*x^n+a1*x^(n-1)+……an
那我们求得系数就可以了。实际上，我们正是通过拉格朗日公式来解决这个问题，原理见上面贴子。
下面是使用该插值算法时所遇到的问题：
1. 真实的心电信号可能都是大于0的值。如果通过所有的离散点进行拟合，则从得到的拟合函数中所获得的插值容易产生负值，一开始不理解，还以为是自己算法问题，后来想想，为了使得到的函数穿过所有这些离散点，得到负值是完全有可能的。这就是常见的拉格朗日问题，当次数较高时，会出现数值不稳定的现象(龙格现象)。解决办法就是分段低阶插值。
2. 我改用了3点插值，3点插值也会出现负值，虽然频率少了一些。因为3点其实构造的就算拋物线，当点的分布不均匀(2,3点离的很近)时，很容易会出现负值。如下：
3. 所以，最后只能采用两点差值来实现趋势的描绘，这是纯折线的，而且就是线性插值，其实完全没有必要用拉格朗日了。
……然后发现拉格朗日在实际中的应用很少，因为存在数值不稳定性的特性。
测量数据会有误差，构造的函数还会有误差，计算机计算浮点数也有误差，所以，最终可能会造成很大的误差。所以，虽然我们有了这个理论基础，但是在实际应用的时候，仍然需要根据特殊的情况进行编写，以适应我们的应用。比如可以让距离间隔差不多的几个点来进行3点插值，这样抛物线会平滑，就不会产生负值了，当然，这只是针对我的这种情况而言。再比如，可以进行3点滑动插值，每次插值完成后，移除第一个点，既每次插值都会保留上次所使用的插值点。
4. 由于拉格朗日插值算法在插值点增加或者减少时，需要重新计算基本多项式，这样效率会很低。为了解决这个问题，引入了牛顿插值法.
拉格朗日代码实现：
//arrX[N],arrY[N]分别存放的是已知插值节点(Xi,Yi)中的Xi,Yi,参数n为已知插值节点的个数,而参数x为待求解的插值节点的横坐标值
float Lagrange(float arrX[],float arrY[],int n, floatx)
{
float yResult = 0.0;
float LValue[80];
intk,m;
floattemp1,temp2;
for(k=0;k
{
temp1 = 1.0;
temp2 = 1.0;
for(m=0;m
{
if(m==k)
{
continue;
}
temp1 *=(x-arrX[m]);
temp2 *=(arrX[k]-arrX[m]);
}
LValue[k]=temp1/temp2;
}
for(int i=0;i
{
yResult+=arrY[i]*LValue[i];
}
returnyResult;
}
二. 牛顿插值算法
牛顿插值法的公式如下：
其中，f[x]的求解如下：
牛顿插值算法的关键在于获得系数f(Xk)。我们就以作者所取的3个点(0, 1), (2, 2), (3, 4)来探讨。根据公式，我们需要取得f(X1), f(X1,X2), f(X1,X2,X3).首先我们分别计算一下值，这是数学计算。
f(X1) = 1；
f(X1,X2) = (f(X2)- f(X1))/(X2- X1) = 1/2
f(X2,X3) = 2. 计算同上，虽然这不是有效系数，但是在进一步求系数的时候会用到。
f(X1,X2,X3) = (f(X2, X3)- f(X2, X3))/(X3- X1) = 1/2
所以，f(X) = 1+1/2*(X-0)+1/2(X-0)(X-2)
这就是理论学习中牛顿插值算法所谓的差商原理。
这一步数学计算的目的可以加深大脑对计算过程的理解，推荐演算一遍。
下面是程序实现，仅仅有理论是不够的，还得写出代码。
首先，我们需要两个数组，X[3],Y[3],分别存放3个横坐标和对应的纵坐标值。
X[3] = {0,2,3}, Y[3] = {1,2,4}.
关键的问题还是求系数，作者给我们提供的思路是使用二维数组来存放，我推演一遍。二维数组a[i][j]的大小根据自己的需要来定。
先给出结果：
当j=0的时候，a[i][0] = Y[i].
当i>=1,j>=1的时候，a[i][j] = (a[i][j-1]-a[i-1][j-1])/(X[i]-X[i-j])
自己归纳一下就可以得到上面的结论。
a[0][0] = Y[0] = 1;
a[1][0] = Y[1] = 2; a[1][1] = (a[1][0]-a[0][0])/(X[1]-X[0]) = (2-1)/(2-0) = 1/2
a[2][0] = Y[2] = 4; a[2][1] = (a[2][0]-a[1][0])/(X[2]-X[1]) = (4-2)/(3-2) = 2
a[2][2] = (a[2][1]-a[1][1])/(X[2]-X[0]) = (2-1/2)/3 = 1/2
所以，f(X) = 1+1/2*(X-0)+1/2(X-0)(X-2)
有效的系数全部存在了a[i][i](i=j的时候)中。
有了上面的二维数组，就可以很容易的实现牛顿插值算法了。
牛顿代码实现：
插值系数计算：
//求得系数，length为已知的插值点个数
for(i=0;i
{
for(j=0;j<=i;j++)
{
if(j==0)
{
a[i][j] =Y[i];
continue;
}
a[i][j] = (a[i][j-1]-a[i-1][j-1])/(X[i]-X[i-j]);
}
}
//牛顿插值算法
float NewTon(float arrX[], float arrY[],float factor[][40],int n,intx)
{
float Result = arrY[0];
float tmp = 1.0;
inti,j;
for(i=1;i
{
for(j=0;j<=i;j++)
{
if(j!=i)
{
continue;
}
tmp *= (x-arrX[i-1]);
Result += factor[i][j]*tmp;
}
}
returnResult;
}
三. 切比雪夫插值
切比雪夫插值的提出是为了改善插值在差值区间上的最大误差控制
关键：切比雪夫插值节点的选取
对于一个插值区间 [a, b], 如果要在这个插值区间上选取n个点作为插值基点,使得上面的最大误差
最小，则基点的选法如下:
对于 i = 1, 2, ...., n. 下面的不等式在区间 [a, b] 上满足
切比雪夫的实现就是在牛顿插值中选取合适的基点，并且要用到cos函数。这个算法感兴趣的可以自己研究一下。
提醒：每一种算法都会有自己的适用范围，我一开始以为这些插值算法可以解决我的需要，但是出来的结果差的太远，两点连线竟是最好的。仔细想想，真正测出的心电图形并不能用一个曲线函数来描述，所以，仅仅依靠几个特征点得出的心电函数，与真实的测量图形相差甚远，所以，做之前考虑清楚才是最重要的。