常见插值算法研究

      最近在做心电相关的项目，由于单片机的处理能力有限，在接收心电信号之后需要对数据进行压缩（其实是取一些特征点），然后后期再进行显示。但是在手持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, float x)
{
     float yResult = 0.0;
     float LValue[80];
     int k,m;
     float temp1,temp2;
     for(k=0;k<n;k++) { temp1 = 1.0; temp2 = 1.0; for(m=0;m<n;m++) { if(m==k) { continue; } temp1 *=(x-arrX[m]); temp2 *=(arrX[k]-arrX[m]); } LValue[k]=temp1/temp2; } for(int i=0;i<n;i++) { yResult+=arrY[i]*LValue[i]; } return yResult; } 

 

 

二. 牛顿插值算法

 

     牛顿插值法的公式如下：

        

     

      其中，f[x]的求解如下：

              

 

     牛顿均值插值法，这篇帖子写的不错， http://blog.csdn.net/johnf_nash/article/details/12622963。耐心看完，自己动手算一下。

     牛顿插值算法的关键在于获得系数 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<length;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,int x) { float Result = arrY[0]; float tmp = 1.0; int i,j; for(i=1;i<n;i++) { for(j=0;j<=i;j++) { if(j!=i) { continue; } tmp *= (x-arrX[i-1]); Result += factor[i][j]*tmp; } } return Result; } 

 

三.   切比雪夫插值

切比雪夫插值的提出是为了改善插值在差值区间上的最大误差控制

                       

      关键：切比雪夫插值节点的选取

               对于一个插值区间 [a, b], 如果要在这个插值区间上选取n个点作为插值基点,使得上面的最大误差

         最小，则基点的选法如下:

                                                

               对于 i = 1, 2, ...., n. 下面的不等式在区间 [a, b] 上满足

                                                    

      切比雪夫的实现就是在牛顿插值中选取合适的基点，并且要用到cos函数。这个算法感兴趣的可以自己研究一下。

 

     提醒：每一种算法都会有自己的适用范围，我一开始以为这些插值算法可以解决我的需要，但是出来的结果差的太远，两点连线竟是最好的。仔细想想，真正测出的心电图形并不能用一个曲线函数来描述，所以，仅仅依靠几个特征点得出的心电函数，与真实的测量图形相差甚远，所以，做之前考虑清楚才是最重要的。

 

转载于:https://www.cnblogs.com/wzhangsdu/p/4198322.html