我正在使用其中一种提议的算法,但结果非常糟糕.

我实施了wiki algorithm

在Java中(代码如下). x(0)是points.get(0),y(0)是值[points.get(0)],α是alfa,μ是mi.其余部分与wiki伪代码相同.

    public void createSpline(double[] values, ArrayList<Integer> points){
    a = new double[points.size()+1];

    for (int i=0; i <points.size();i++)
    {
        a[i] = values[points.get(i)];

    }

    b = new double[points.size()];
    d = new double[points.size()];
    h = new double[points.size()];

    for (int i=0; i<points.size()-1; i++){
        h[i] = points.get(i+1) - points.get(i);
    }

    alfa = new double[points.size()];

    for (int i=1; i <points.size()-1; i++){
        alfa[i] = (double)3 / h[i] * (a[i+1] - a[i]) - (double)3 / h[i-1] *(a[i+1] - a[i]);
    }

    c = new double[points.size()+1];
    l = new double[points.size()+1];
    mi = new double[points.size()+1];
    z = new double[points.size()+1];

    l[0] = 1; mi[0] = z[0] = 0;

    for (int i =1; i<points.size()-1;i++){
        l[i] = 2 * (points.get(i+1) - points.get(i-1)) - h[i-1]*mi[i-1];
        mi[i] = h[i]/l[i];
        z[i] = (alfa[i] - h[i-1]*z[i-1])/l[i];
    }

    l[points.size()] = 1;
    z[points.size()] = c[points.size()] = 0;

    for (int j=points.size()-1; j >0; j--)
    {
        c[j] = z[j] - mi[j]*c[j-1];
        b[j] = (a[j+1]-a[j]) - (h[j] * (c[j+1] + 2*c[j])/(double)3) ;
        d[j] = (c[j+1]-c[j])/((double)3*h[j]);
    }


    for (int i=0; i<points.size()-1;i++){
        for (int j = points.get(i); j<points.get(i+1);j++){
            //                fk[j] = values[points.get(i)];
            functionResult[j] = a[i] + b[i] * (j - points.get(i)) 
                                + c[i] * Math.pow((j - points.get(i)),2)
                                + d[i] * Math.pow((j - points.get(i)),3);
        }
    }

}

我得到的结果如下：

但它应该类似于：

我也试图按照另一种方式实现算法：http://www.geos.ed.ac.uk/~yliu23/docs/lect_spline.pdf

首先,他们展示了如何做线性样条,这很容易.我创建了计算A和B系数的函数.然后,他们通过添加二阶导数来扩展线性样条. C和D系数也很容易计算.

但是当我试图计算二阶导数时,问题就出现了.我不明白他们是如何计算的.

所以我只实现了线性插值.结果是：

有谁知道如何修复第一个算法或解释我如何计算第二个算法中的二阶导数？

解决方法:

最近,Unser,Thévenaz等人在一系列论文中描述了立方b样条,其中包括

M. Unser, A. Aldroubi, M. Eden, “Fast B-Spline Transforms for Continuous Image Representation and Interpolation”, IEEE Trans. Pattern Anal. Machine Intell., vol. 13, n. 3, pp. 277-285, Mar. 1991.

M. Unser, “Splines, a Perfect Fit for Signal and Image Processing”, IEEE Signal Proc. Mag., pp. 22- 38, Nov. 1999.

和

P. Thévenaz, T. Blu, M. Unser, “Interpolation Revisited,” IEEE Trans. on Medical Imaging, vol. 19, no. 7, pp. 739-758, July 2000.

以下是一些指导原则.

什么是样条线？

样条曲线是平滑连接在一起的分段多项式.对于度数n的样条,每个段是n次多项式.连接这些部件使得样条连续到其在结的n-1度的导数,即多项式部件的连接点.

如何构造样条线？

零阶样条曲线如下

所有其他样条曲线都可以构造为

卷积采取n-1次.

立方样条

最流行的样条是三次样条,其表达式为

样条插值问题

给定在离散整数点k处采样的函数f(x),样条插值问题是确定以下列方式表示的近似值s(x)到f(x)

其中ck是插值系数,s(k)= f(k).

样条预过滤

不幸的是,从n = 3开始,

这样ck不是插值系数.它们可以通过求解通过强制s(k)= f(k)得到的线性方程组来确定,即

这样的方程可以以卷积形式重铸并且在变换的z空间中求解为

哪里

因此,

以这种方式进行总是优于通过例如LU分解提供线性方程组的解.

可以通过注意到来确定上述等式的解

哪里

第一部分代表因果过滤器,而第二部分代表反义过滤器.它们都在下图中说明.

因果过滤器

Anticausal过滤器

在上图中,

滤波器的输出可以用以下递归方程表示

可以通过首先确定c-和c的“初始条件”来解决上述等式.在假定周期性镜像输入序列fk时

那么可以证明c的初始条件可以表示为

而c的初始条件可表示为