做视频编解码算法的话一般都会接触到一个目前比较通用的客观评价指标，即 Bjøntegaard Delta（BD），如果关注码率变化的话就是 BD-Rate，如果关注 PSNR 变化的话就是BD-PSNR，其实质就是求在一个指标相同时另外一个指标的平均变化。一般我们关注的是在同等客观质量下的码率变化，所以用得比较多的是 BD-Rate。不过问题是怎么去求平均，一种简单的方法就是假设码率（R）和 PSNR（D）具有一定的函数关系，即 R = f ( D ) R=f(D) R=f(D) 或者 D = g ( R ) D=g( R) D=g(R)，一般来说认为是单调的，那么两种编解码算法就对应了两条曲线，要求码率或者 PSNR 的平均变化就可以通过求区间积分的差值来表示了。不过问题又来了，编解码算法那么复杂我们怎么去得到两者的函数关系，至于它们两者有没有函数关系都说不准，所以我们只能再进行一些简化，通常就是测几个 QP 的数据，得到几个离散的点，然后进行插值作为该曲线的估计。所以 BD-Rate 归根结底还是一个插值问题。插值的方法有很多，不过我们还要考虑它们的复杂度以及合理性，其中多项式插值倒是一个不错的选择。多项式插值我们听得最多的应该就是 Lagrange 插值、分段三次插值以及样条（Spline）插值了，所以借着学习 BD-Rate 的机会顺便把这些算法的原理搞懂，以后用到的地方可能也不少。

1 多项式插值算法

给定一系列已知点 ( x k , y k ) (x_k, y_k) (xk,yk)，其中 k = 1 , 2 , . . . , n k=1, 2, ..., n k=1,2,...,n，且 x k x_k xk 各不相同， y k y_k yk 不全为零，不失一般性假设 x k x_k xk 从小到大排列，插值算法要做的是对其他非 x k x_k xk 位置的值进行合理的估计，同时要保持 x k x_k xk 位置的值与 y k y_k yk 一致。由于非 x k x_k xk 位置的值未知，实际上我们可以定义无数种插值方法，然而其中绝大多数得到的结果都是不合常理的，我们也不会接受这些结果。由于多项式函数的形式比较简单，而且易于求导和积分，所以在实际应用中被大量使用。

1.1 Lagrange 插值法

可以证明，给定一系列已知点 ( x k , y k ) (x_k, y_k) (xk,yk)，其中 k = 1 , 2 , . . . , n k=1, 2, ..., n k=1,2,...,n，且 x k x_k xk 各不相同，我们总可得到且仅能得到一个唯一的阶数等于或小于 n − 1 n-1 n−1 的多项式 P ( x ) P(x) P(x)，使得

P ( x k ) = y k , k = 1 , 2 , . . . , n . P({x_k}) = {y_k},{\text{ }}k = 1,2,...,n. P(xk)=yk, k=1,2,...,n.

其可表示为 Lagrange 形式，即

P ( x ) = ∑ i = 1 n ( ∏ j ≠ i x − x j x i − x j ) ⋅ y i . (1-1) P(x) = \sum\limits_{i = 1}^n {\left( {\prod\limits_{j \ne i} {\frac{{x - {x_j}}}{{{x_i} - {x_j}}}} } \right) \cdot {y_i}} . \tag{1-1} P(x)=i=1∑n j=i∏xi−xjx−xj ⋅yi.(1-1)

可以看出， P ( x ) P(x) P(x) 求和项内每一项都为最高阶数为 n − 1 n-1 n−1 的多项式，由于各 x k x_k xk 互不相等，因此不会存在断点，函数在整个实数域中是连续且可导的。另外，也容易证明，当 x = x k x=x_k x=xk 时，求和项中只有 k k k-th 这一项不为 0，且其连乘项为 1，因此有 P ( x k ) = y k P(x_k)=y_k P(xk)=yk。由于阶数等于或小于 n − 1 n-1 n−1 的多项式是唯一的，式 (1-1) 即为我们所求的插值函数，其他未知点的值也可以根据该函数进行计算估计。

进一步，我们注意到式 (1-1) 所定义的 Lagrange 形式使用十分不方便，因此我们更希望将其简化为常用的多项式形式，即

P ( x ) = ∑ i = 0 n − 1 a i x i = a 0 + a 1 x + ⋯ a n − 1 x n − 1 . (1-2) P(x) = \sum\limits_{i = 0}^{n - 1} {{a_i}{x^i}} = {a_0} + {a_1}x + \cdots {a_{n - 1}}{x^{n - 1}}.\tag{1-2} P(x)=i=0∑n−1aixi=a0+a1x+⋯an−1xn−1.(1-2)

我们可以看到，式 (1-2) 中有 n n n 个未知数，而我们刚好有 n n n 个不同的点刚好组成 n n n 个方程，写成矩阵形式有

( x 1 n − 1 ⋯ x 1 1 x 2 n − 1 ⋯ x 2 1 ⋮ ⋱ ⋮ ⋮ x n n − 1 ⋯ x n 1 ) ( a n − 1 ⋮ a 1 a 0 ) = X a = y = ( y 1 y 2 ⋮ y n ) . (1-3) \left( {\begin{array}{c} {x_1^{n - 1}}& \cdots &{{x_1}}&1\\ {x_2^{n - 1}}& \cdots &{{x_2}}&1\\ \vdots & \ddots & \vdots & \vdots \\ {x_n^{n - 1}}& \cdots &{{x_n}}&1 \end{array}} \right)\left( {\begin{array}{c} {{a_{n - 1}}}\\ \vdots \\ {{a_1}}\\ {{a_0}} \end{array}} \right) = {\bf{Xa}} = {\bf{y}} = \left( {\begin{array}{c} {{y_1}}\\ {{y_2}}\\ \vdots \\ {{y_n}} \end{array}} \right).\tag{1-3} x1n−1x2n−1⋮xnn−1⋯⋯⋱⋯x1x2⋮xn11⋮1 an−1⋮a1a0 =Xa=y= y1y2⋮yn .(1-3)

其中 X X X 称为范德蒙 (Vandermonde) 矩阵，即其元素

v k , j = x k n − j . {v_{k,j}} = x_k^{n - j}. vk,j=xkn−j.

可以证明，范德蒙矩阵的秩等于不同的 x k x_k xk 个数，也就是 R ( X ) = n R(X)=n R(X)=n。推广到非方阵，即假设不要求多项式的阶数最高为 n − 1 n-1 n−1，也就是可以考虑更低阶或高阶的多项式，例如最高阶为 m − 1 m-1 m−1，总共有 m m m 个系数，那么矩阵 X X X 的秩 R ( X ) = m i n { n , m } R(X)=min\{n, m\} R(X)=min{n,m}。则：

1） n > m n>m n>m，那么 m = R ( X ) ≤ R ( X , y ) ≤ m + 1 m=R(X) \le R(X, y) \le m+1 m=R(X)≤R(X,y)≤m+1。注意第一个 ≤ \le ≤ 号。当所有样点都来自于同一个最高阶数为 m − 1 m-1 m−1 的多项式时，无论 n n n 有多大，总有 X a = y Xa=y Xa=y 成立，即有 R ( X ) = R ( X , y ) = m R(X)=R(X, y)=m R(X)=R(X,y)=m，所以这时方程组也有唯一解。但绝大多数情况下 X a = y Xa=y Xa=y 是不成立的，这时 R ( X ) < R ( X , y ) = m + 1 R(X)<R(X, y)=m+1 R(X)<R(X,y)=m+1，方程组无解，只能通过最小二乘法进行优化，使得 ∣ ∣ y − X a ∣ ∣ ||y-Xa|| ∣∣y−Xa∣∣ 最小，这时候称为拟合问题。

2） n = m n=m n=m，那么 n = R ( X ) = R ( X , y ) = m n=R(X)=R(X, y)=m n=R(X)=R(X,y)=m，方程有唯一解，这也证明了最高阶为 n − 1 n-1 n−1 的多项式的是唯一的，因此通过式 (1-2) 的待定系数法和式 (1-1) 的 Lagrange 法求得的多项式是等价的。由于多项式系数个数与样点个数一样，这时也称为全阶 (Full-degree) 多项式。注意有时得到的解中 a n − 1 = 0 a_{n-1}=0 an−1=0，即多项式的最高阶小于 n − 1 n-1 n−1，问题就退化到（1）中所有样点都来自于一个低阶多项式的情况。因此直接假设一个最高阶为 n − 1 n-1 n−1 的多项式比人为设置一个低阶多项式更容易有解。

3） n < m n<m n<m，那么 n = R ( X ) = R ( X , y ) < m n=R(X)=R(X, y)<m n=R(X)=R(X,y)<m，方程有无限多解。也就是说，我们可以找到无限多的最高阶数大于等于 n n n 的多项式 P ( x ) P(x) P(x)，使得 P ( x k ) = y k P(x_k)=y_k P(xk)=yk。但一般我们不会这么做，因为这容易造成过拟合的情况，特别是 n n n 比较大的情况。

经过以上的分析可知，我们总可以得到一个经过所有已知点的最高阶数等于或小于 n − 1 n-1 n−1 的多项式，这似乎已经完美解决了我们需要求解的插值问题。然而，我们应该注意到，当 n n n 比较大时，我们求得的多项式阶数也相应比较高，但在实际应用中我们很少会用到非常高阶的多项式，因为这在自然界中十分少见，很容易造成过拟合。另外由于计算机的存储问题，高阶函数非常容易造成数据溢出以及误差扩散问题。还有，由于阶数问题，式 (1-3) 所定义的范德蒙矩阵左右两侧数据量级是极其不平衡的，准确来说其条件数会非常大，求解该方程组的问题属于病态 (ill-conditioning) 问题，虽然我们总可以得到解，但 y y y 的微小变化，就可能导致求得的多项式系数变化剧烈，而由于模数转换以及计算机存储的精度问题， y y y 的微小变化是很常见的。

图1 Lagrange 法插值结果

举个例子，假设有以下数据：

x = [ 1 , 2 , 3 , 4 , 5 , 6 ] ; y = [ 16 , 18 , 21 , 17 , 15 , 12 ] . x = [1, 2, 3, 4, 5, 6]; \\ y = [16, 18, 21, 17, 15, 12]. x=[1,2,3,4,5,6];y=[16,18,21,17,15,12].

利用 Lagrange 法可求得从高阶到低阶的多项式系数为

a = [ 0.2417 , 4.3333 , − 28.9583 , 87.6667 , − 115.8 , 69.0 ] . a = [0.2417, 4.3333, -28.9583, 87.6667, -115.8, 69.0]. a=[0.2417,4.3333,−28.9583,87.6667,−115.8,69.0].

如果对 y y y 添加一个细小的均匀分布的随机噪声，如

y ’ = [ 16.34 , 18.08 , 20.74 , 17.35 , 14.57 , 11.59 ] , y’ = [16.34, 18.08, 20.74, 17.35, 14.57, 11.59], y’=[16.34,18.08,20.74,17.35,14.57,11.59],

这时可得到的多项式系数为

a ’ = [ − 0.1766 , 3.3192 , − 21.8785 , 67.1505 , − 89.3781 , 57.404 ] . a’ = [-0.1766, 3.3192, -21.8785, 67.1505, -89.3781, 57.404]. a’=[−0.1766,3.3192,−21.8785,67.1505,−89.3781,57.404].

可以看到，尽管 y y y 只发生了微小变化，得到的多项式系数却是截然不同的，两者对应的函数图像如图 1 所示。另外，从图 1 可以看出，得到的多项式函数明显有些过拟合了，或者说在某些区间变化过于剧烈，例如在区间 [ 1 , 2 ] [1, 2] [1,2] 和 [ 5 , 6 ] [5, 6] [5,6]，这种过度的变化称为超调量或者过冲 (overshoot)，有可能破坏原本数据的变化趋势，例如区间内的单调性等。因此，对于已知样点数较多的情况下，我们一般不会使用 Lagrange 法进行插值以避免过高的多项式阶数。

1.2 Hermite 插值法

在 1.1 节中我们已经证明最高阶数大于等于 n n n 的插值多项式有无限多个，这是由于待定系数个数大于样点数也就是约束方程个数的原因，因此，通过增加约束方程，我们可以确定具有更高阶数的多项式的具体表达式。因为样点数已经确定，我们只能从各个样点的导数入手，这就是 Hermite 插值法的基本思路。

Hermite 插值法要求找到一个多项式 P ( x ) P(x) P(x)，使得该多项式在各个样点处的函数值以及至少达到 p p p 阶的导数值与给定值相等，也就是

P ( j ) ( x k ) = y k ( j ) , k = 1 , 2 , . . . , n , j = 0 , 1 , . . . , p . {P^{(j)}}({x_k}){\rm{ }} = y_k^{(j)},{\rm{ }}k = 1,2,...,n,{\rm{ }}j = 0,1,...,p. P(j)(xk)=yk(j),k=1,2,...,n,j=0,1,...,p.

由于多项式的各阶导数依然是多项式，那么我们总共可以得到 ( p + 1 ) n (p+1)n (p+1)n 个线性约束方程。为了方程组有且有唯一解，我们需要定义一个最高阶数为 ( p + 1 ) n − 1 (p+1)n-1 (p+1)n−1 的多项式。可以看到，Lagrange 插值法实际是 p = 0 p=0 p=0 的 Hermite 插值法。

在前面我们已经讨论过，过高的阶数会导致超调量的产生，而 Hermite 插值法每增加一阶的已知导数就可能增加 n n n 阶的多项式分量，因此通常也不会直接用于多个数据的插值。但是，如果考虑 n = 2 , p = 1 n=2, p=1 n=2,p=1 的情况，即仅有两个点的以及相应一阶导数，我们可得到一个最高阶数为 3 的多项式，称为三次 Hermite 插值多项式 (Cubic Hermite Interpolation Polynomial, CHIP)。由于其阶数较低，不容易产生较大的超调量，同时又能保持一定平滑性，所以被广泛用于数据的分段插值，该多项式也被称为分段三次 Hermite 插值多项式 (Piecewise CHIP, PCHIP)，具体在后面会有详细介绍。

1.2.1 CHIP 的具体形式

下面我们来推导 CHIP 的具体形式，假设有两点 ( x 0 , y 0 ) (x_0, y_0) (x0,y0) 和 ( x 1 , y 1 ) (x_1, y_1) (x1,y1)，其中 x 0 < x 1 x_0<x_1 x0<x1，假设其一阶导数存在，且分别为 d 0 d_0 d0 和 d 1 d_1 d1，令

P ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 , x ∈ [ x 0 , x 1 ] . (1-4) P(x) = {a_0} + {a_1}x + {a_2}{x^2} + {a_3}{x^3},{\rm{ }}x \in [{x_0},{x_1}].\tag{1-4} P(x)=a0+a1x+a2x2+a3x3,x∈[x0,x1].(1-4)

则有

P ( x 0 ) = a 0 + a 1 x 0 + a 2 x 0 2 + a 3 x 0 3 = y 0 , P ( x 1 ) = a 0 + a 1 x 1 + a 2 x 1 2 + a 3 x 1 3 = y 1 , P ′ ( x 0 ) = a 1 + 2 a 2 x 0 + 3 a 3 x 0 2 = d 0 , P ′ ( x 1 ) = a 1 + 2 a 2 x 1 + 3 a 3 x 1 2 = d 1 . (1-5) \begin{array}{l} P({x_0}) = {a_0} + {a_1}{x_0} + {a_2}x_0^2 + {a_3}x_0^3 = {y_0},\\ \\ P({x_1}) = {a_0} + {a_1}{x_1} + {a_2}x_1^2 + {a_3}x_1^3 = {y_1},\\ \\ P'({x_0}) = {a_1} + 2{a_2}{x_0} + 3{a_3}x_0^2 = {d_0},\\ \\ P'({x_1}) = {a_1} + 2{a_2}{x_1} + 3{a_3}x_1^2 = {d_1}. \end{array}\tag{1-5} P(x0)=a0+a1x0+a2x02+a3x03=y0,P(x1)=a0+a1x1+a2x12+a3x13=y1,P′(x0)=a1+2a2x0+3a3x02=d0,P′(x1)=a1+2a2x1+3a3x12=d1.(1-5)

通常我们不使用式 (1-4) 的形式，令

h = x 1 − x 0 , s = x − x 0 , δ = ( y 1 − y 0 ) / ( x 1 − x 0 ) , h = {x_1} - {x_0},{\rm{ }}s = x - {x_0},{\rm{ }}\delta = ({y_1} - {y_0})/({x_1} - {x_0}), h=x1−x0,s=x−x0,δ=(y1−y0)/(x1−x0),

可得

P ( x ) = H 1 ( x ) y 0 + H 2 ( x ) y 1 + H 3 ( x ) d 0 + H 4 ( x ) d 1 = h 3 − 3 h s 2 + 2 s 3 h 3 y 0 + 3 h s 2 − 2 s 3 h 3 y 1 + s ( s − h ) 2 h 2 d 0 + s 2 ( s − h ) h 2 d 1 = y 0 + s d 0 + s 2 ( 3 δ − 2 d 0 − d 1 ) / h + s 3 ( d 0 − 2 δ + d 1 ) / h 2 = y 0 + s d 0 + s 2 c 0 + s 3 b 0 . (1-6) \begin{aligned} P(x) &= {H_1}(x){y_0} + {H_2}(x){y_1} + {H_3}(x){d_0} + {H_4}(x){d_1}\\ &= \frac{{{h^3} - 3h{s^2} + 2{s^3}}}{{{h^3}}}{y_0} + \frac{{3h{s^2} - 2{s^3}}}{{{h^3}}}{y_1} + \frac{{s{{\left( {s - h} \right)}^2}}}{{{h^2}}}{d_0} + \frac{{{s^2}\left( {s - h} \right)}}{{{h^2}}}{d_1}\\ &= {y_0} + s{d_0} + {s^2}\left( {3\delta - 2{d_0} - {d_1}} \right)/h + {s^3}\left( {{d_0} - 2\delta + {d_1}} \right)/{h^2}\\ &= {y_0} + s{d_0} + {s^2}{c_0} + {s^3}{b_0}. \end{aligned}\tag{1-6} P(x)=H1(x)y0+H2(x)y1+H3(x)d0+H4(x)d1=h3h3−3hs2+2s3y0+h33hs2−2s3y1+h2s(s−h)2d0+h2s2(s−h)d1=y0+sd0+s2(3δ−2d0−d1)/h+s3(d0−2δ+d1)/h2=y0+sd0+s2c0+s3b0.(1-6)

其中 H k ( x ) H_k(x) Hk(x)称为三次 Herimite 的基函数。同样，我们很容易可推导其一阶、二阶导数以及定积分的具体形式，即

P ′ ( x ) = d 0 + 2 s c 0 + 3 s 2 b 0 , P ′ ′ ( x ) = 2 c 0 + 6 s b 0 , ∫ m n P ( x ) d x = ( s y 0 + 1 2 s 2 d 0 + 1 3 s 3 c 0 + 1 4 s 4 b 0 + C ) ∣ m − x 0 n − x 0 . (1-7) \begin{aligned} P'(x) &= {d_0} + 2s{c_0} + 3{s^2}{b_0},{\rm{ }}\\ P''(x) &= 2{c_0} + 6s{b_0},\\ \int_m^n {P(x)dx} &= \left. {\left( {s{y_0} + {\textstyle{1 \over 2}}{s^2}{d_0} + {\textstyle{1 \over 3}}{s^3}{c_0} + {\textstyle{1 \over 4}}{s^4}{b_0} + C} \right)} \right|_{m - {x_0}}^{n - {x_0}}. \end{aligned}\tag{1-7} P′(x)P′′(x)∫mnP(x)dx=d0+2sc0+3s2b0,=2c0+6sb0,=(sy0+21s2d0+31s3c0+41s4b0+C) m−x0n−x0.(1-7)

1.2.2 CHIP 的单调性

在一些应用中，我们希望在区间 [ x 0 , x 1 ] [x_0, x_1] [x0,x1] 上的插值是单调的，这样可以尽量避免数值的剧烈变化，减少超调量，使插值的结果符合人眼的直觉。在式 (1-7) 中，我们知道 CHIP 的一阶导数是一个二次函数，且与 d 0 d_0 d0 和 d 1 d_1 d1 的取值密切相关，不当的选择很有可能造成极值点的产生，破坏插值的单调性。因此我们有必要对 CHIP 的单调条件进行讨论。

首先，我们很容易能够得到 P ( x ) P(x) P(x) 在区间 [ x 0 , x 1 ] [x_0, x_1] [x0,x1] 上单调的必要条件为

s i g n ( d 0 ) = s i g n ( d 1 ) = s i g n ( δ ) , s i g n ( x ) = { 1 , x > 0 − 1 , x < 0 (1-8) sign({d_0}) = sign({d_1}) = sign(\delta ),\\ \\ sign(x) = \left\{ {\begin{array}{c} {1,}&{x > 0}\\ { - 1,}&{x < 0} \end{array}} \right.\tag{1-8} sign(d0)=sign(d1)=sign(δ),sign(x)={1,−1,x>0x<0(1-8)

其中 0 的符号位一般为 0，但有时与 1 或者 -1 比较时会不那么严格地认为是相等的。如果 δ = 0 \delta=0 δ=0，那么 P ( x ) P(x) P(x) 单调的充分必要条件为 d 0 = d 1 = 0 d_0=d_1=0 d0=d1=0，也就是说 P ( x ) P(x) P(x) 是常数。在下面的讨论中我们假设 δ ≠ 0 \delta \ne 0 δ=0。

P ( x ) P(x) P(x) 的导数的具体形式为

P ′ ( x ) = d 0 + 2 s ( 3 δ − 2 d 0 − d 1 ) / h + 3 s 2 ( d 0 − 2 δ + d 1 ) / h 2 . (1-9) P'(x) = {d_0} + 2s\left( {3\delta - 2{d_0} - {d_1}} \right)/h + 3{s^2}\left( {{d_0} - 2\delta + {d_1}} \right)/{h^2}.\tag{1-9} P′(x)=d0+2s(3δ−2d0−d1)/h+3s2(d0−2δ+d1)/h2.(1-9)

那么，

1）如果 d 0 + d 1 − 2 δ = 0 d_0+d_1-2\delta=0 d0+d1−2δ=0，则 P ( x ) P(x) P(x) 为二次(或线性)函数， P ′ ( x ) P'(x) P′(x) 则是线性的 (或常数)，于是有 min ⁡ { d 0 , d 1 } ≤ P ′ ( x ) ≤ max ⁡ { d 0 , d 1 } . \min \{ {d_0},{d_1}\} \le P'(x) \le \max \{ {d_0},{d_1}\} . min{d0,d1}≤P′(x)≤max{d0,d1}. 结合 (1-8) 的必要条件， P ′ ( x ) P'(x) P′(x) 总是大于 0 或小于 0，可知 d 0 + d 1 − 2 δ = 0 d_0+d_1-2\delta=0 d0+d1−2δ=0 是 P ( x ) P(x) P(x) 单调的一个充分条件。

2）如果 d 0 + d 1 − 2 δ ≠ 0 d_0+d_1-2\delta \ne0 d0+d1−2δ=0，则 P ′ ( x ) P'(x) P′(x) 为二次函数。当 δ > 0 \delta>0 δ>0 且 d 0 + d 1 − 2 δ < 0 d_0+d_1-2\delta<0 d0+d1−2δ<0 时， P ′ ( x ) P'(x) P′(x) 开口向下，且在区间 [ x 0 , x 1 ] [x_0, x_1] [x0,x1] 上有 0 ≤ min ⁡ { d 0 , d 1 } ≤ P ′ ( x ) 0 \le \min \{ {d_0},{d_1}\} \le P'(x) 0≤min{d0,d1}≤P′(x)，因此 P ( x ) P(x) P(x) 是单调递增的。同样地，当 δ < 0 \delta<0 δ<0 且 d 0 + d 1 − 2 δ > 0 d_0+d_1-2\delta>0 d0+d1−2δ>0 时， P ′ ( x ) P'(x) P′(x) 开口向上，且有 P ′ ( x ) ≤ max ⁡ { d 0 , d 1 } ≤ 0 P'(x) \le \max \{ {d_0},{d_1}\} \le 0 P′(x)≤max{d0,d1}≤0，因此 P ( x ) P(x) P(x) 是单调递减的。

为了整合 (2) 的两种情况，定义 α = d 0 / δ , β = d 1 / δ \alpha=d_0/\delta, \beta=d_1/\delta α=d0/δ,β=d1/δ，注意当 (1-8) 满足时， α , β \alpha,\beta α,β 都是非负的。那么当 α + β − 2 < 0 \alpha+\beta-2<0 α+β−2<0 时， P ( x ) P(x) P(x) 是单调的。因此，由以上的分析可得，当满足 (1-8) 的必要条件时，如果 α + β − 2 ≤ 0 \alpha+\beta-2 \le 0 α+β−2≤0 时，则 P ( x ) P(x) P(x) 是单调的。

α + β − 2 ≤ 0 \alpha+\beta-2 \le 0 α+β−2≤0 的条件要求 α , β \alpha, \beta α,β 不能太大，那么当 α + β − 2 > 0 \alpha+\beta-2 > 0 α+β−2>0 时我们需要做进一步讨论。由于 P ′ ( x ) P'(x) P′(x) 为二次函数，因此不考虑区间时有唯一一个极值点

x ∗ = x 0 + h 3 [ 2 α + β − 3 α + β − 2 ] . (1-10) {x^ * } = {x_0} + \frac{h}{3}\left[ {\frac{{2\alpha + \beta - 3}}{{\alpha + \beta - 2}}} \right].\tag{1-10} x∗=x0+3h[α+β−22α+β−3].(1-10)

且有

P ′ ( x ∗ ) = φ ( α , β ) δ , φ ( α , β ) = α − 1 3 ( 2 α + β − 3 ) 2 ( α + β − 2 ) . (1-11) P'({x^ * }) = \varphi \left( {\alpha ,\beta } \right)\delta ,{\text{ }}\varphi \left( {\alpha ,\beta } \right) = \alpha - \frac{1}{3}\frac{{{{(2\alpha + \beta - 3)}^2}}}{{(\alpha + \beta - 2)}}.\tag{1-11} P′(x∗)=φ(α,β)δ, φ(α,β)=α−31(α+β−2)(2α+β−3)2.(1-11)

那么，当满足 α + β − 2 > 0 \alpha+\beta-2 > 0 α+β−2>0 和式 (1-8) 时，如果以下任意一个条件成立， P ( x ) P(x) P(x) 是单调的：

1 ) x ∗ ∉ ( x 0 , x 1 ) ; 2 ) x ∗ ∈ ( x 0 , x 1 ) , a n d s i g n ( P ′ ( x ∗ ) ) = s i g n ( δ ) . \begin{aligned} &{\rm{1) }}{x^ * } \notin \left( {{x_0},{x_1}} \right);\\ &{\rm{2) }}{x^ * } \in \left( {{x_0},{x_1}} \right),{\rm{ and }}sign\left( {P'({x^ * })} \right) = {\rm{sign}}(\delta ). \end{aligned} 1)x∗∈/(x0,x1);2)x∗∈(x0,x1),andsign(P′(x∗))=sign(δ).

对于 (1)，当 x ∗ ≤ x 0 x^* \le x_0 x∗≤x0 时，相当于 2 α + β − 3 ≤ 0 2\alpha+\beta-3 \le 0 2α+β−3≤0； x ∗ ≥ x 1 x^* \ge x_1 x∗≥x1 时，相当于 α + 2 β − 3 ≤ 0 \alpha+2\beta-3 \le 0 α+2β−3≤0。而对于 (2)，也就是 φ ( α , β ) ≥ 0 \varphi(\alpha, \beta) \ge 0 φ(α,β)≥0。因此，当满足 α + β − 2 > 0 \alpha+\beta-2>0 α+β−2>0 和式 (1-8) 时，如果以下任意一个条件成立， P ( x ) P(x) P(x) 是单调的：

1 ) 2 α + β − 3 ≤ 0 ; 2 ) α + 2 β − 3 ≤ 0 ; 3 ) φ ( α , β ) ≥ 0. \begin{aligned} &1){\rm{ }}2\alpha + \beta - 3 \le 0;\\ &2){\rm{ }}\alpha + 2\beta - 3 \le 0;\\ &3){\rm{ }}\varphi \left( {\alpha ,\beta } \right) \ge 0. \end{aligned} 1)2α+β−3≤0;2)α+2β−3≤0;3)φ(α,β)≥0.

经过以上分析，我们可以得到使 P ( x ) P(x) P(x) 在区间 [ x 0 , x 1 ] [x_0, x_1] [x0,x1] 上单调的 α , β \alpha,\beta α,β 的取值范围，其图像如图 2 所示。其中 φ ( α , β ) = 0 \varphi(\alpha, \beta) = 0 φ(α,β)=0 是一个椭圆，与坐标轴在点 (0, 3) 和 (3, 0) 相切，其表达式为

( α − 1 ) 2 + ( α − 1 ) ( β − 1 ) + ( β − 1 ) 2 − 3 ( α + β − 2 ) = 0. {\left( {\alpha - 1} \right)^2} + \left( {\alpha - 1} \right)\left( {\beta - 1} \right) + {\left( {\beta - 1} \right)^2} - 3\left( {\alpha + \beta - 2} \right) = 0. (α−1)2+(α−1)(β−1)+(β−1)2−3(α+β−2)=0.

在图 2 中，对角阴影线代表的是 α + β − 2 ≤ 0 \alpha+\beta-2 \le 0 α+β−2≤0，垂直阴影线代表的是 α + β − 2 > 0 \alpha+\beta-2>0 α+β−2>0 且 2 α + β − 3 ≤ 0 2\alpha+\beta-3 \le 0 2α+β−3≤0，水平阴影线代表的是 α + β − 2 > 0 \alpha+\beta-2>0 α+β−2>0 且 α + 2 β − 3 ≤ 0 \alpha+2\beta-3 \le 0 α+2β−3≤0，点阴影代表的是 φ ( α , β ) ≥ 0 \varphi(\alpha, \beta) \ge 0 φ(α,β)≥0。其他空白区域则会使得 P ( x ) P(x) P(x) 非单调。注意 α + β − 2 > 0 \alpha+\beta-2>0 α+β−2>0 时有部分区域是重叠的。

图2 使 P(x) 单调的 α (横), β (纵) 的取值范围(阴影)

2 分段多项式插值

在前面的分析中我们已经知道，对于已知样点较多的情况，Lagrange 插值法会带来很高的多项式阶数，容易造成过冲现象，而且极不利于计算机计算，因此一个显而易见的解决方法是将插值区间分为若干个子区间，然后在各个区间内分别使用多项式插值，这样就可以极大地降低所用多项式的阶数。然而，分段低阶多项式则意味着某些样点处的高阶导数的不连续，这有可能带来一些主观上或者工业制造的问题，因此根据不同的需求也有很多不同的分段多项式插值方法。通常来说，综合考虑复杂度以及实用性，以两个样点作为一个子区间的分段插值方法比较常用，以下的讨论都基于这种分段方法。

2.1 分段线性插值

最简单的分段多项式插值莫过于分段线性插值，即直接使用一条直线连接相邻的两个样点。每个子区间内的直线实际上是 n = 2 n=2 n=2 的 Lagrange 多项式或者 n = 2 , p = 0 n=2, p=0 n=2,p=0 的 Hermite 多项式。给定一系列已知点 ( x k , y k ) (x_k, y_k) (xk,yk)，其中 k = 1 , 2 , . . . , n k=1, 2, ..., n k=1,2,...,n，且 x k x_k xk 各不相同，可得每个子区间内的函数表达式为

L ( x ) = y k + ( x − x k ) y k + 1 − y k x k + 1 − x k , x ∈ [ x k , x k + 1 ] , k = 1 , 2 , . . . , n − 1. (2-1) L(x) = {y_k} + (x - {x_k})\frac{{{y_{k + 1}} - {y_k}}}{{{x_{k + 1}} - {x_k}}},{\rm{ }}x \in \left[ {{x_k},{x_{k + 1}}} \right],{\rm{ }}k = 1,2,...,n - 1.\tag{2-1} L(x)=yk+(x−xk)xk+1−xkyk+1−yk,x∈[xk,xk+1],k=1,2,...,n−1.(2-1)

图3 分段线性插值结果

对于 1.1 节中的例子，可得其图像如图 3 所示。可以看到，分段线性插值具有完美的区间单调性，十分符合人眼对数据变化趋势的直觉。然而，在每个样点处，只有函数是连续的，而一阶导数则产生了跳变，即导数是不连续的，而曲线的平滑度与其导数的连续性相关，因此整体看起来非常不平滑。

实际上，在多项式插值问题中，区间的单调性和整体的平滑性是两个较难调和的问题，上述的分段线性插值和 1.1 节中的Lagrange 插值，则是区间单调性和整体平滑性的两个极端。分段线性插值具有很好的区间单调性，却具有十分糟糕的平滑性；而Lagrange 插值则保证了所具有的任何阶数的导数在样点处连续，是一个光滑函数，然而其具有很差的区间单调性，很容易出现过冲的现象。因此，为了调和两者，我们需要在分段的同时保证一定阶数的导数的连续性，这就需要用到分段 Hermite 插值，常用的是 n = 2 , p = 1 n=2, p=1 n=2,p=1 的分段三次插值多项式(PCHIP)，即至少保证各样点处的一阶导数是存在的。然而，不同的一阶导数的选择也就会导致不同的插值结果，常见的有两种，一种称为形状保持的 PCHIP (Shape-Preserving PCHIP, SPPCHIP)，其要求在一阶导数连续的同时，保证每个区间内是单调变化的，其单调性条件在 1.2.2 节中已经有详细讨论；另一种则是分段三次样条 (Spline) 插值，其不要求区间内的单调性，但要求各样点处的二阶导数是连续的，不过并不要求二阶导数有确定取值，不同的边界条件也会导致不同的插值结果。

2.2 形状保持的 PCHIP

SPPCHIP 是在保持区间单调性的同时，使得样点处的一阶导数连续，从而使得插值后的曲线更加平滑，从而在视觉更加符合人们对数据变化趋势的预期。在 1.2.2 中我们已经推导出在保证一阶导数存在的条件下使得各区间单调变化的一阶导数的取值范围，但我们没有给出具体的导数求解方法。SPPCHIP 最初是由 Fred Fritsch 等人在 1980 年左右提出的，其将各样点的一阶导数值设置为左右两侧割线斜率的加权调和平均值 (Weighted Harmonic Mean)，使得其刚好落在了图 2 中单调取值范围中，从而能够简单快速地求得各个样点的一阶导数取值，并且保证每个子区间是单调变化的。

具体来说，给定一系列已知点 ( x k , y k ) (x_k, y_k) (xk,yk)，其中 k = 1 , 2 , . . . , n k=1, 2, ..., n k=1,2,...,n，且 x k x_k xk 各不相同且从小到大排列，令

h k = x k + 1 − x k , s k = x − x k , δ k = ( y k + 1 − y k ) / ( x k + 1 − x k ) , (2-2) {h_k} = {x_{k + 1}} - {x_k},{\rm{ }}{s_k} = x - {x_k},{\rm{ }}{\delta _k} = ({y_{k + 1}} - {y_k})/({x_{k + 1}} - {x_k}),\tag{2-2} hk=xk+1−xk,sk=x−xk,δk=(yk+1−yk)/(xk+1−xk),(2-2)

那么非边界各个样点 d k d_k dk, k = 2 , 3 , … , n − 1 k=2,3,…,n-1 k=2,3,…,n−1 的一阶导数取值如下：

1）如果 δ k − 1 \delta_{k-1} δk−1 和 δ k \delta_k δk 的符号相反或者至少其中一个为 0，那么令 x k x_k xk 为局部极值点，即令 d k = 0 d_k =0 dk=0。

2）如果 δ k − 1 \delta_{k-1} δk−1 和 δ k \delta_k δk 具有相同的符号且都不为 0，则令 d k d_k dk 为两者的加权调和平均值，即

w 1 + w 2 d k = w 1 δ k − 1 + w 2 δ k , w 1 = 2 h k + h k − 1 , w 2 = h k + 2 h k − 1 . (2-3) \frac{{{w_1} + {w_2}}}{{{d_k}}} = \frac{{{w_1}}}{{{\delta _{k - 1}}}} + \frac{{{w_2}}}{{{\delta _k}}},{\rm{ }}{w_1} = 2{h_k} + {h_{k - 1}},{\rm{ }}{w_2} = {h_k} + 2{h_{k - 1}}.\tag{2-3} dkw1+w2=δk−1w1+δkw2,w1=2hk+hk−1,w2=hk+2hk−1.(2-3)

对于边界样点 d 1 d_1 d1 和 d n d_n dn，则有

d k = ( 2 h p + h q ) δ p − h p δ q h p + h q , if sign ( d k ) ≠ s i g n ( δ p ) , d k = 0 ; elseif sign ( δ p ) ≠ s i g n ( δ q ) a n d ∣ d k ∣ > ∣ 3 δ p ∣ , d k = 3 δ p . (2-4) \begin{aligned} &{d_k} = \frac{{\left( {2{h_p} + {h_q}} \right){\delta _p} - {h_p}{\delta _q}}}{{{h_p} + {h_q}}},{\rm{ }}\\ &{\text{if sign}}\left( {{d_k}} \right) \ne {\rm{sign}}\left( {{\delta _p}} \right),{d_k} = 0;\\ &{\text{elseif sign}}\left( {{\delta _p}} \right) \ne {\rm{sign}}\left( {{\delta _q}} \right){\rm{ and }}\left| {{d_k}} \right| > \left| {3{\delta _p}} \right|,{d_k} = 3{\delta _p}. \end{aligned}\tag{2-4} dk=hp+hq(2hp+hq)δp−hpδq,if sign(dk)=sign(δp),dk=0;elseif sign(δp)=sign(δq)and∣dk∣>∣3δp∣,dk=3δp.(2-4)

其中，当 k = 1 k=1 k=1 时， p = 1 , q = 2 p=1, q=2 p=1,q=2；当 k = n k=n k=n 时， p = n , q = n − 1 p=n, q=n-1 p=n,q=n−1。

从上面可知，各样点的一阶导数的取值符合式 (1-8) 的必要条件，下面证明其充分性。令

α k = d k / δ k , β k = d k + 1 / δ k , {\alpha _k} = {d_k}/{\delta _k},{\rm{ }}{\beta _k} = {d_{k + 1}}/{\delta _k}, αk=dk/δk,βk=dk+1/δk,

对于非边界区间 [ x k , x k + 1 ] , k = 2 , 3 , … , n − 2 [x_k, x_k+1], k=2,3,…,n-2 [xk,xk+1],k=2,3,…,n−2， d k d_k dk 要么为 0，要么为式 (2-3) 所定义。对于 d k = 0 d_k=0 dk=0，则有 α k = 0 \alpha_k=0 αk=0 以及 β k − 1 = 0 \beta_{k-1}=0 βk−1=0。对于式 (2-3)，有

α k = d k δ k = w 1 + w 2 ( w 1 δ k − 1 + w 2 δ k ) δ k < w 1 + w 2 w 2 = 1 + 2 h k + h k − 1 h k + 1 2 h k − 1 + 3 2 h k − 1 < 3 , β k − 1 = d k δ k − 1 = w 1 + w 2 ( w 1 δ k − 1 + w 2 δ k ) δ k − 1 < w 1 + w 2 w 1 = 3 ( h k + h k − 1 ) h k + h k − 1 + h k < 3. (2-5) \begin{aligned} {\alpha _k} = \frac{{{d_k}}}{{{\delta _k}}} = \frac{{{w_1} + {w_2}}}{{\left( {\frac{{{w_1}}}{{{\delta _{k - 1}}}} + \frac{{{w_2}}}{{{\delta _k}}}} \right){\delta _k}}} < \frac{{{w_1} + {w_2}}}{{{w_2}}} = 1 + \frac{{2{h_k} + {h_{k - 1}}}}{{{h_k} + {\textstyle{1 \over 2}}{h_{k - 1}} + {\textstyle{3 \over 2}}{h_{k - 1}}}} < 3,\\ \\ {\beta _{k - 1}} = \frac{{{d_k}}}{{{\delta _{k - 1}}}} = \frac{{{w_1} + {w_2}}}{{\left( {\frac{{{w_1}}}{{{\delta _{k - 1}}}} + \frac{{{w_2}}}{{{\delta _k}}}} \right){\delta _{k - 1}}}} < \frac{{{w_1} + {w_2}}}{{{w_1}}} = \frac{{3\left( {{h_k} + {h_{k - 1}}} \right)}}{{{h_k} + {h_{k - 1}} + {h_k}}} < 3. \end{aligned}\tag{2-5} αk=δkdk=(δk−1w1+δkw2)δkw1+w2<w2w1+w2=1+hk+21hk−1+23hk−12hk+hk−1<3,βk−1=δk−1dk=(δk−1w1+δkw2)δk−1w1+w2<w1w1+w2=hk+hk−1+hk3(hk+hk−1)<3.(2-5)

对于左边界区间 [ x 1 , x 2 ] [x_1, x_2] [x1,x2]，根据式 (2-4)，当 δ 1 \delta_1 δ1 和 δ 2 \delta_2 δ2 符号相同且都不为 0 时，有

α 1 = d 1 δ 1 = ( 2 h 1 + h 2 ) δ 1 − h 1 δ 2 ( h 1 + h 2 ) δ 1 < 2 h 1 + h 2 h 1 + h 2 < 2. (2-6) {\alpha _1} = \frac{{{d_1}}}{{{\delta _1}}} = \frac{{\left( {2{h_1} + {h_2}} \right){\delta _1} - {h_1}{\delta _2}}}{{\left( {{h_1} + {h_2}} \right){\delta _1}}} < \frac{{2{h_1} + {h_2}}}{{{h_1} + {h_2}}} < 2.\tag{2-6} α1=δ1d1=(h1+h2)δ1(2h1+h2)δ1−h1δ2<h1+h22h1+h2<2.(2-6)

这时 α 1 \alpha_1 α1 可能小于 0，但是根据式 (2-4)，如果 d 1 d_1 d1 的符号不等于 δ 1 \delta_1 δ1，那么 d 1 = 0 d_1=0 d1=0，从而 α 1 = 0 α_1=0 α1=0。另外。如果 δ 1 \delta_1 δ1 和 δ 2 \delta_2 δ2 符号不同，式 (2-6) 的不等式就不成立了，但最终 d 1 d_1 d1 被限制在 3 δ 1 3\delta_1 3δ1 的范围内。对于右边界区间 [ x n − 1 , x n ] [x_{n-1}, x_n] [xn−1,xn] 对 β n − 1 \beta_{n-1} βn−1 的分析同理，这里不赘述。

从上面分析可知，各个区间内的 α \alpha α 和 β \beta β 都是在图 2 的阴影范围内的，所以经过插值后每个区间内都是单调变化的。因此，SPPCHIP 给出了一种简便的求解各样点一阶导数的方法，可以使插值的结果在各个区间内单调变化，或者说视觉上令人舒适的 (visually pleasing)。对于 1.1 节的例子，可得其经过 SPPCHIP 插值的结果如图 4 所示。可以看出每个区间内的插值结果是单调变化的，其变化趋势比较符合人的直觉，另外分段函数的连接处也相对平滑。注意，SPPCHIP 并不能保证区间外的插值结果具有单调性。

图4 SPPCHIP 插值结果

2.3 分段三次样条插值

分段三次样条插值也属于 PCHIP，其要求每个样点处的二阶导数也是存在的，由于样点左右两侧的分段函数发生了改变，这些样点也被称为扭结点 (knot)。形状保持的 PCHIP 虽然能获得视觉上比较满意的插值结果，但其只具有一阶导数的连续性，所以分段曲线的连接处还是不够平滑，而在很多工业应用场景中，例如汽车零件制造，其平滑性要求分段曲线至少保证二阶导数连续，这就是三次样条插值在 PCHIP 家族中应用最为广泛的原因。

由于每个区间内的分段函数为三次多项式，经过两次求导后我们可以得到一个线性函数，根据式 (1-6) 和 (1-7) 可得

P ′ ′ ( x ) = ( 6 h k − 12 s k ) δ k + ( 6 s k − 2 h k ) d k + 1 + ( 6 s k − 4 h k ) d k h k 2 . (2-7) P''\left( x \right) = \frac{{\left( {6{h_k} - 12{s_k}} \right){\delta _k} + \left( {6{s_k} - 2{h_k}} \right){d_{k + 1}} + \left( {6{s_k} - 4{h_k}} \right){d_k}}}{{h_k^2}}.\tag{2-7} P′′(x)=hk2(6hk−12sk)δk+(6sk−2hk)dk+1+(6sk−4hk)dk.(2-7)

那么，当 x = x k x=x_k x=xk 时， s k = 0 s_k=0 sk=0，有

P ′ ′ ( x k + ) = 6 δ k − 2 d k + 1 − 4 d k h k . (2-8) P''\left( {{x_k} + } \right) = \frac{{6{\delta _k} - 2{d_{k + 1}} - 4{d_k}}}{{{h_k}}}.\tag{2-8} P′′(xk+)=hk6δk−2dk+1−4dk.(2-8)

同理，当x=xk+1时，sk=hk，有

P ′ ′ ( x k + 1 − ) = − 6 δ k + 4 d k + 1 + 2 d k h k . (2-9) P''\left( {{x_{k + 1}} - } \right) = \frac{{ - 6{\delta _k} + 4{d_{k + 1}} + 2{d_k}}}{{{h_k}}}.\tag{2-9} P′′(xk+1−)=hk−6δk+4dk+1+2dk.(2-9)

为了使二阶各样点处的导数连续，那么对于内点（即非边界点）有，

P ′ ′ ( x k − ) = P ′ ′ ( x k + ) , k = 2 , 3 , . . . , n − 1. (2-10) P''\left( {{x_k} - } \right) = P''\left( {{x_k} + } \right),{\rm{ }}k = 2,3,...,n - 1.\tag{2-10} P′′(xk−)=P′′(xk+),k=2,3,...,n−1.(2-10)

代入式 (2-8) 和 (2-9) 可得

h k d k − 1 + 2 ( h k − 1 + h k ) d k + h k − 1 d k + 1 = 3 ( h k δ k − 1 + h k − 1 δ k ) . (2-11) {h_k}{d_{k - 1}} + 2\left( {{h_{k - 1}} + {h_k}} \right){d_k} + {h_{k - 1}}{d_{k + 1}} = 3\left( {{h_k}{\delta _{k - 1}} + {h_{k - 1}}{\delta _k}} \right).\tag{2-11} hkdk−1+2(hk−1+hk)dk+hk−1dk+1=3(hkδk−1+hk−1δk).(2-11)

由此，我们就可得到关于 d 1 , d 2 , … , d n d_1, d_2, …, d_n d1,d2,…,dn 的 n − 2 n-2 n−2 个线性方程，很明显这有无穷多个解。但是这时我们还未对两个边界点添加约束，因此不同的边界条件将会得到不同的插值结果。

一种常用的边界约束是称为 not-a-knot 的方法，也就是区间 [ x 1 , x 3 ] [x_{1}, x_{3}] [x1,x3] 和 [ x n − 2 , x n ] [x_{n-2}, x_n] [xn−2,xn] 实际上各属于一条三次多项式曲线，这样 x 2 x_2 x2 和 x n − 2 x_{n-2} xn−2 两侧的函数并没有发生改变，也就不属于扭结点 (knot)。这样， x 2 x_2 x2 和 x n − 2 x_{n-2} xn−2 上的三次导数也是连续的，即有

P ′ ′ ′ ( x 2 − ) = P ′ ′ ′ ( x 2 + ) , P ′ ′ ′ ( x n − 2 − ) = P ′ ′ ′ ( x n − 2 + ) . (2-12) P'''\left( {{x_2} - } \right) = P'''\left( {{x_2} + } \right),{\rm{ }}P'''\left( {{x_{n - 2}} - } \right) = P'''\left( {{x_{n - 2}} + } \right). \tag{2-12} P′′′(x2−)=P′′′(x2+),P′′′(xn−2−)=P′′′(xn−2+).(2-12)

先考虑 k = 2 k=2 k=2的情况，可得

d 1 − 2 δ 1 + d 2 h 1 2 = d 2 − 2 δ 2 + d 3 h 2 2 ⇒ h 2 2 d 1 − ( h 1 2 − h 2 2 ) d 2 − h 1 2 d 3 = 2 h 2 2 δ 1 − 2 h 1 2 δ 2 . (2-13) \begin{aligned} &\frac{{{d_1} - 2{\delta _1} + {d_2}}}{{h_1^2}} = \frac{{{d_2} - 2{\delta _2} + {d_3}}}{{h_2^2}} \Rightarrow \\ &h_2^2{d_1} - \left( {h_1^2 - h_2^2} \right){d_2} - h_1^2{d_3} = 2h_2^2{\delta _1} - 2h_1^2{\delta _2}. \end{aligned}\tag{2-13} h12d1−2δ1+d2=h22d2−2δ2+d3⇒h22d1−(h12−h22)d2−h12d3=2h22δ1−2h12δ2.(2-13)

虽然以式 (2-13) 的形式加入方程组中并没有什么问题，但我们仔细分析式 (2-11)，可知其组成的矩阵是一个三对角矩阵，即只有主对角线和相邻两条对角线上的值非零，可以采用一些优化算法进行求解。而对于式 (2-13)，其实际上包含了三个未知数，而三对角矩阵第一行和最后一行应该只包含前两个未知数，因此为了构造一个三对角矩阵，需要对式 (2-13) 进行处理。综合式 (2-11) 和 (2-13)，可得

{ h 2 2 d 1 − ( h 1 2 − h 2 2 ) d 2 − h 1 2 d 3 = 2 h 2 2 δ 1 − 2 h 1 2 δ 2 h 2 d 1 + 2 ( h 1 + h 2 ) d 2 + h 1 d 3 = 3 ( h 2 δ 1 + h 1 δ 2 ) ⇒ h 2 d 1 + ( h 1 + h 2 ) d 2 = ( 3 h 1 + 2 h 2 ) h 2 δ 1 + h 1 2 δ 2 h 1 + h 2 = r 1 . (2-14) \begin{aligned} \left\{ {\begin{array}{c} {h_2^2{d_1} - \left( {h_1^2 - h_2^2} \right){d_2} - h_1^2{d_3} = 2h_2^2{\delta _1} - 2h_1^2{\delta _2}}\\ {{h_2}{d_1} + 2\left( {{h_1} + {h_2}} \right){d_2} + {h_1}{d_3} = 3\left( {{h_2}{\delta _1} + {h_1}{\delta _2}} \right)} \end{array}} \right. \Rightarrow \\ {h_2}{d_1} + \left( {{h_1} + {h_2}} \right){d_2} = \frac{{\left( {3{h_1} + 2{h_2}} \right){h_2}{\delta _1} + h_1^2{\delta _2}}}{{{h_1} + {h_2}}} = {r_1}. \end{aligned}\tag{2-14} {h22d1−(h12−h22)d2−h12d3=2h22δ1−2h12δ2h2d1+2(h1+h2)d2+h1d3=3(h2δ1+h1δ2)⇒h2d1+(h1+h2)d2=h1+h2(3h1+2h2)h2δ1+h12δ2=r1.(2-14)

同理可得

( h n − 1 + h n − 2 ) d n − 1 + h n − 2 d n = h n − 1 2 δ n − 2 + ( 3 h n − 1 + 2 h n − 2 ) h n − 2 δ n − 1 h n − 2 + h n − 1 = r n . (2-15) \left( {{h_{n - 1}} + {h_{n - 2}}} \right){d_{n - 1}} + {h_{n - 2}}{d_n} = \frac{{h_{n - 1}^2{\delta _{n - 2}} + \left( {3{h_{n - 1}} + 2{h_{n - 2}}} \right){h_{n - 2}}{\delta _{n - 1}}}}{{{h_{n - 2}} + {h_{n - 1}}}} = {r_n}.\tag{2-15} (hn−1+hn−2)dn−1+hn−2dn=hn−2+hn−1hn−12δn−2+(3hn−1+2hn−2)hn−2δn−1=rn.(2-15)

这时，我们总共得到 n n n 个约束方程，对应 n n n 个未知数即各个样点的一阶导数。写成矩阵形式可得

A d = r . (2-16) {\bf{Ad}} = {\bf{r}}.\tag{2-16} Ad=r.(2-16)

其中

A = ( h 2 h 1 + h 2 h 2 2 ( h 1 + h 2 ) h 1 h 3 2 ( h 2 + h 3 ) h 2 ⋱ ⋱ ⋱ h n − 1 2 ( h n − 2 + h n − 1 ) h n − 2 h n − 1 + h n − 2 h n − 2 ) , (2-17) {\bf{A}} = \left( {\begin{array}{c} {{h_2}}&{{h_1} + {h_2}}&{}&{}&{}&{}\\ {{h_2}}&{2\left( {{h_1} + {h_2}} \right)}&{{h_1}}&{}&{}&{}\\ {}&{{h_3}}&{2\left( {{h_2} + {h_3}} \right)}&{{h_2}}&{}&{}\\ {}&{}& \ddots & \ddots & \ddots &{}\\ {}&{}&{}&{{h_{n - 1}}}&{2\left( {{h_{n - 2}} + {h_{n - 1}}} \right)}&{{h_{n - 2}}}\\ {}&{}&{}&{}&{{h_{n - 1}} + {h_{n - 2}}}&{{h_{n - 2}}} \end{array}} \right),\tag{2-17} A= h2h2h1+h22(h1+h2)h3h12(h2+h3)⋱h2⋱hn−1⋱2(hn−2+hn−1)hn−1+hn−2hn−2hn−2 ,(2-17)

d = ( d 1 d 2 ⋮ d n ) , r = ( r 1 3 ( h 2 δ 1 + h 1 δ 2 ) 3 ( h 3 δ 2 + h 2 δ 3 ) ⋮ 3 ( h n − 1 δ n − 2 + h n − 2 δ n − 1 ) r n ) . (2-18) {\bf{d}} = \left( {\begin{array}{c} {{d_1}}\\ {{d_2}}\\ \vdots \\ {{d_n}} \end{array}} \right),{\rm{ }}{\bf{r}} = \left( {\begin{array}{c} {{r_1}}\\ {3\left( {{h_2}{\delta _1} + {h_1}{\delta _2}} \right)}\\ {3\left( {{h_3}{\delta _2} + {h_2}{\delta _3}} \right)}\\ \vdots \\ {3\left( {{h_{n - 1}}{\delta _{n - 2}} + {h_{n - 2}}{\delta _{n - 1}}} \right)}\\ {{r_n}} \end{array}} \right).\tag{2-18} d= d1d2⋮dn ,r= r13(h2δ1+h1δ2)3(h3δ2+h2δ3)⋮3(hn−1δn−2+hn−2δn−1)rn .(2-18)

可以看到 A A A 是一个三对角 (tridiagonal) 矩阵，区别于传统矩阵有一些特殊的优化求解算法，如追赶法或者 Thomas 法等，在 MATLAB 中可用 tridisolve() 函数，在 Python 中可用 scipy.linalg.solve_banded() 函数进行求解。

对于 1.1 节的例子，可得其经过分段三次样条插值的结果如图 5 所示。可以看到，由于没有限制各样点一阶导数的取值在图 2 的阴影范围内，各区间内的插值结果有可能不是单调变化的，但其相比于 SPPCHIP 插值的结果显得更加平滑，但是不像 Lagrange 插值那样产生很严重的 overshoot 现象。因此 spline 插值在一些不严格要求区间单调性的问题中应用十分广泛。

图5 分段三次 Spline 插值结果

2.4 小结

在前面的内容中，我们分析了几种常用的多项式插值算法，分别从其区间单调性以及整体的平滑性进行了详细的讨论，图 6 展示了几种方法插值的结果。由于这是插值问题，所以几种方法插值的结果都会准确地经过已知样点。但在非样点部分，它们的表现各不相同。

图6 几种多项式插值方法的结果比较

Lagrange 插值得到的结果是一个全局的多项式，因此具有无限阶连续的导数，所以非常光滑，但在某些区间，特别是靠近边界的区间，其非常容易产生较大的 overshoot 现象，使插值的结果难以让人信服。作为另外一种极端，分段线性插值的一阶导数在各样点处发生跳变，因此看上去极不光滑，在行业生产制造中不适合应用，但其形式简单，并且能够保证各个区间内的插值结果是单调变化的，因此在图像处理算法中也经常有使用。作为前两者两个极端的调和，PCHIP 增加了各样点出一阶导数连续的条件，因此相比于分段线性插值显得更加平滑，同时其保持了较低的多项式阶数，相比于 Lagrange 插值又不那么容易产生较大的超调量。在 PCHIP 中，SPPCHIP 保证了各区间的单调性，因此其与分段线性插值又有相似之处，在图 6 也可看到两者的结果比较相近，但 SPPCHIP 明显更加平滑；而 Spline 则不保证区间的单调性，但增加二阶导数连续的条件，所以相比于 SPPCHIP 又显得更加平滑。因此，哪种多项式插值算法更加适合，需要根据具体的应用需求具体分析。

3 Bjøntegaard Delta

视频编解码算法中简单而常用的两个客观评价指标是平均码率 (Rate) 和峰值信噪比 (PSNR)。由于视频编码过程中涉及很多的优化方法，修改某个参数就有可能影响整个率失真优化 (RDO) 的结果，因此难以直接反映 Rate-PSNR 的函数关系。目前通常的做法是固定其他参数配置，仅修改初始 QP 值，分别测试得到各初始 QP 配置下的 (Rate, PSNR) 对，然后对各个 (Rate, PSNR) 对之间的部分采用某些插值算法进行插值，这样就可得到该模式下的 Rate-PSNR 曲线，其也可以作为 R-D 曲线的一种形式。但由于测试的工作量比较大，通常只选用 4 个典型的 QP 值，H.26x 中常采用的是 22，27，32，37，其他编解码算法由于量化方法的不同可能会选择不同的 QP 值，但为了方便比较通常要求其码率或者 PSNR 变化范围相近。图 7 展示了两种编解码算法的结果。

图7 两种编解码算法的 Rate-PSNR 曲线
图8 两种编解码算法的 log(Rate)-PSNR 曲线

在得到两种编解码算法 Rate-PSNR 曲线后，从直观上我们可以对两者的性能进行比较，如果其中一种算法的 Rate-PSNR 曲线在另外一种算法曲线的上方，如图 7 中的 Plot1，那么很明显该算法是更优的，因为其在相同码率消耗的同时能够获得更高的 PSNR 即更好的客观图像质量，或者说在保证同样的 PSNR 的同时需要更低的码率消耗。但是，从图 7 中我们无法知道两者具体的性能差异，这就无法制定一个客观的评价标准来指导我们改善算法。对此，Gisle Bjøntegaard 等人在 H.264 标准开发过程中提出了 Bjøntegaard Delta 指标，这也是目前视频编解码算法中常用的客观评价标准。

根据不同的角度，Bjøntegaard Delta 可分为 BD-PSNR 和 BD-Rate 两种。BD-PSNR 反映了相同码率下平均的 PSNR(dB) 差值，而 BD-Rate 反映了相同 PSNR 下平均的码率变化比例 (%)。因为人眼对一定范围内 PSNR 的变化并不敏感，但码率的变化却会切实影响到硬盘存储或者网络传输，因此 BD-Rate 往往更受关注。由于 Bjøntegaard Delta 反映的是整体的平均变化量，我们需要通过区间积分来进行估计，这就要求 Rate-PSNR 曲线是连续可积的。另外，从图 7 中我们可以看到，由于 PSNR 计算的是峰值信噪比，当码率较高而失真较小时，较小的失真变化都有可能导致较大的 PSNR 变化，如果直接在图 7 中计算 BD-PSNR 或者 BD-Rate 有可能导致高码率部分占据主导因素而不能够准确反映低码率的情况，因此 Bjøntegaard 提出在 log(Rate) 尺度上进行计算，如图 8 所示。这时，两条曲线接近于直线（但不一定）且近乎平行，这样就能够更好地反映低码率的情况，并且近似直线通常也更加利于插值的准确性和稳定性。另外，采用 log(Rate) 尺度也能够更加方便 BD-Rate 的计算，下面会具体分析。

为了计算 BD-PSNR 和 BD-Rate，我们首先要对图 8 的曲线进行插值。从直观上我们可以发现 log(Rate)-PSNR 曲线应该是单调变化的，因此在选择插值算法时应该尽量保持插值结果的单调性。由于只有四个样点，如果我们采用 Lagrange 插值，那么我们获得的是一条全局的三次多项式曲线。因为阶数较低，虽然我们在前面的内容中已经讨论过 Lagrange 插值有可能出现严重的 overshoot 现象，但由于从图 8 中可以看出样点排列相对整齐，所以插值的结果还是比较可信的，因此在最初的版本中，Bjøntegaard Delta 使用的也是 Lagrange 插值方法。但在目前的版本中，为了保证其单调性，通常采用的是形状保持的分段三次 Hermite 插值法 (SPPCHIP)，因此在求解定积分时需要分段计算。

为了简便，令 R R R 为 Rate， D D D 为 PSNR， r = l o g ( R ) r=log(R) r=log(R)，其中 l o g ( ) log() log() 函数以 10 为底。为了计算 BD-PSNR，假设通过插值算法可得 D = f ( r ) D=f(r) D=f(r)，那么可得

Δ D = E ( D 2 − D 1 ) ≈ 1 r H − r L ∫ r L r H [ f 2 ( r ) − f 1 ( r ) ] d r , r L = max ⁡ { min ⁡ { r 1 } , min ⁡ { r 2 } } , r H = min ⁡ { max ⁡ { r 1 } , max ⁡ { r 2 } } . (3-1) \begin{aligned} &\Delta D = E\left( {{D_2} - {D_1}} \right) \approx \frac{1}{{{r_H} - {r_L}}}\int_{{r_L}}^{{r_H}} {\left[ {{f_2}\left( r \right) - {f_1}\left( r \right)} \right]dr} ,\\ &{r_L} = \max \left\{ {\min \left\{ {{{\bf{r}}_1}} \right\},{\rm{ }}\min \left\{ {{{\bf{r}}_2}} \right\}} \right\},{\rm{ }}{r_H} = \min \left\{ {\max \left\{ {{{\bf{r}}_1}} \right\},{\rm{ }}\max \left\{ {{{\bf{r}}_2}} \right\}} \right\}. \end{aligned}\tag{3-1} ΔD=E(D2−D1)≈rH−rL1∫rLrH[f2(r)−f1(r)]dr,rL=max{min{r1},min{r2}},rH=min{max{r1},max{r2}}.(3-1)

同理，为了计算 BD-Rate，假设通过插值算法可得 r = g ( D ) r=g(D) r=g(D)，注意 f f f 和 g g g 并不是反函数关系，那么可得

Δ R = E ( R 2 − R 1 R 1 ) = E ( R 2 R 1 ) − 1 = E ( 10 r 2 − r 1 ) − 1 ≈ 1 0 E ( r 2 − r 1 ) − 1 ≈ 1 0 1 D H − D L ∫ D L D H [ g 2 ( D ) − g 1 ( D ) ] d D − 1 , D L = max ⁡ { min ⁡ { D 1 } , min ⁡ { D 2 } } , D H = min ⁡ { max ⁡ { D 1 } , max ⁡ { D 2 } } . (3-2) \begin{aligned} \Delta R &= E\left( {\frac{{{R_2} - {R_1}}}{{{R_1}}}} \right) = E\left( {\frac{{{R_2}}}{{{R_1}}}} \right) - 1 = E\left( {{{10}^{{r_2} - {r_1}}}} \right) - 1\\ &\approx {10^{E\left( {{r_2} - {r_1}} \right)}} - 1 \approx {10^{\frac{1}{{{D_H} - {D_L}}}\int_{{D_L}}^{{D_H}} {\left[ {{g_2}\left( D \right) - {g_1}\left( D \right)} \right]dD} }} - 1,\\ {D_L} &= \max \left\{ {\min \left\{ {{{\bf{D}}_1}} \right\},{\rm{ }}\min \left\{ {{{\bf{D}}_2}} \right\}} \right\},{\rm{ }}{D_H} = \min \left\{ {\max \left\{ {{{\bf{D}}_1}} \right\},{\rm{ }}\max \left\{ {{{\bf{D}}_2}} \right\}} \right\}. \end{aligned}\tag{3-2} ΔRDL=E(R1R2−R1)=E(R1R2)−1=E(10r2−r1)−1≈10E(r2−r1)−1≈10DH−DL1∫DLDH[g2(D)−g1(D)]dD−1,=max{min{D1},min{D2}},DH=min{max{D1},max{D2}}.(3-2)

其中， r , D r, D r,D 代表测试所得的几个样点。从式 (3-1) 和 (3-2) 可以看到，一般当两条曲线完全平行时约等号才可以替换为等号，而根据图 8 两条曲线在高码率时的差距明显要比低码率要大一些，因此 Bjøntegaard Delta 在两个编码器高低码率表现相差较大时可能不那么准确。另外，为了使积分区间包含大部分曲线，两条曲线的区间范围应该尽可能相近。

参考资料

https://www.mathworks.com/content/dam/mathworks/mathworks-dot-com/moler/interp.pdf
Monotone Piecewise Cubic Interpolation.
https://zero.sci-hub.tw/1536/cb55e0aa6f5eba4f009e7e812f4fc953/fritsch1980.pdf#view=FitH
https://coast.nd.edu/jjwteach/www/www/30125/pdfnotes/lecture5_9v14.pdf
https://blogs.mathworks.com/cleve/2012/07/16/splines-and-pchips/
Calculation of average coding efficiency based on subjective quality scores.
https://infoscience.epfl.ch/record/192802/files/Elsevier2013_preprint.pdf
http://www.uta.edu/faculty/krrao/dip/Courses/EE5351/BDPSNR.doc

Python脚本

# -*- coding: utf-8 -*-from __future__ import print_function, division
import numpy as np
import scipydef pchip_slopes(h, delta):d = np.zeros(len(h) + 1)k = np.argwhere(np.sign(delta[:-1]) * np.sign(delta[1:]) > 0).reshape(-1) + 1w1 = 2*h[k] + h[k-1]w2 = h[k] + 2*h[k-1]d[k] = (w1 + w2) / (w1 / delta[k-1] + w2 / delta[k])d[0] = pchip_end(h[0], h[1], delta[0], delta[1])d[-1] = pchip_end(h[-1], h[-2], delta[-1], delta[-2])return ddef pchip_end(h1, h2, del1, del2):d = ((2*h1 + h2)*del1 - h1*del2) / (h1 + h2)if np.sign(d) != np.sign(del1):d = 0elif np.sign(del1) != np.sign(del2) and np.abs(d) > np.abs(3*del1):d = 3 * del1return ddef spline_slopes(h, delta):a, r = np.zeros([3, len(h)+1]), np.zeros(len(h)+1)a[0, 1], a[0, 2:] = h[0] + h[1], h[:-1]a[1, 0], a[1, 1:-1], a[1, -1] = h[1], 2*(h[:-1] + h[1:]), h[-2]a[2, :-2], a[2, -2] = h[1:], h[-1] + h[-2]r[0] = ((h[0] + 2*a[0, 1])*h[1]*delta[0] + h[0]**2*delta[1]) / a[0, 1]r[1:-1] = 3*(h[1:] * delta[:-1] + h[:-1] * delta[1:])r[-1] = (h[-1]**2*delta[-2] + (2*a[2, -2] + h[-1])*h[-2]*delta[-1]) / a[2, -2]d = scipy.linalg.solve_banded((1, 1), a, r)return dclass PCHIP:def __init__(self, x, y, use_spline=False):assert len(np.unique(x)) == len(x)order = np.argsort(x)self.xi, self.yi = x[order], y[order]h = np.diff(self.xi)delta = np.diff(self.yi) / hself.d = spline_slopes(h, delta) if use_spline else pchip_slopes(h, delta)self.c = (3*delta - 2*self.d[:-1] - self.d[1:]) / hself.b = (self.d[:-1] - 2*delta + self.d[1:]) / h**2"""The piecewise function is like p(x) = y_k + s*d_k + s*s*c_k + s*s*s*b_kwhere s = x - xi_k, k is the interval includeing x.So the original function of p(x) is P(x) = s*y_k + 1/2*s*s*d_k + 1/3*s*s*s*c_k + 1/4*s*s*s*s*b_k + C."""self.interval_int_coeff = []self.interval_int = np.zeros(len(x)-1)for i in range(len(x)-1):self.interval_int_coeff.append(np.polyint([self.b[i], self.c[i], self.d[i], self.yi[i]]))self.interval_int[i] = np.polyval(self.interval_int_coeff[-1], h[i]) - np.polyval(self.interval_int_coeff[-1], 0)def interp(self, x):if len(x[x < np.min(self.xi)]) > 0 or len(x[x > np.max(self.xi)]) > 0:print('Warning: Some samples are out of the interval and the results may be strange!')"""find the intervals the xs belong to"""k = np.zeros(len(x), dtype='int')for i in range(1, len(self.xi)-1):idx = np.argwhere(self.xi[i] <= x).reshape(-1)k[idx] = i * np.ones(len(idx), dtype='int')s = x - self.xi[k]y = self.yi[k] + s*(self.d[k] + s*(self.c[k] + s*self.b[k]))return ydef _integral(self, lower, upper):assert lower <= upperif lower < np.min(self.xi):lower = np.min(self.xi)print('Warning: The lower bound is less than the interval and clipped!')elif lower > np.max(self.xi):print('Warning: The lower bound is greater than the interval!')return 0if upper > np.max(self.xi):upper = np.max(self.xi)print('Warning: The upper bound is greater than the interval and clipped!')elif upper < np.min(self.xi):print('Warning: The lower bound is less than the interval!')return 0left = np.arange(len(self.xi))[self.xi - lower > -1e-6][0]right = np.arange(len(self.xi))[self.xi - upper < 1e-6][-1]inte = np.sum(self.interval_int[left:right])if self.xi[left] - lower > 1e-6:inte += (np.polyval(self.interval_int_coeff[left-1], self.xi[left]-self.xi[left-1]) - np.polyval(self.interval_int_coeff[left-1], lower-self.xi[left-1]))if self.xi[right] - upper < -1e-6:inte += (np.polyval(self.interval_int_coeff[right], upper-self.xi[right]) - np.polyval(self.interval_int_coeff[right], 0))return intedef integral(self, lower, upper):if lower > upper:return -self._integral(upper, lower)else:return self._integral(lower, upper)def BD_Rate(R1, PSNR1, R2, PSNR2, piecewise=True):lR1, lR2 = np.log10(R1), np.log10(R2)min_int = np.max((np.min(PSNR1), np.min(PSNR2)))max_int = np.min((np.max(PSNR1), np.max(PSNR2)))if piecewise == True:int1 = PCHIP(PSNR1, lR1, use_spline=False).integral(min_int, max_int)int2 = PCHIP(PSNR2, lR2, use_spline=False).integral(min_int, max_int)else:p1 = np.polyfit(PSNR1, lR1, len(lR1)-1)p2 = np.polyfit(PSNR2, lR2, len(lR1)-1)p_int1 = np.polyint(p1)p_int2 = np.polyint(p2)int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)avg_exp_diff = (int2 - int1) / (max_int - min_int)avg_diff = (np.power(10, avg_exp_diff) - 1) * 100.return avg_diffdef BD_PSNR(R1, PSNR1, R2, PSNR2, piecewise=True):lR1, lR2 = np.log10(R1), np.log10(R2)min_int = np.max((np.min(lR1), np.min(lR2)))max_int = np.min((np.max(lR1), np.max(lR2)))if piecewise == True:int1 = PCHIP(lR1, PSNR1, use_spline=False).integral(min_int, max_int)int2 = PCHIP(lR2, PSNR2, use_spline=False).integral(min_int, max_int)else:p1 = np.polyfit(lR1, PSNR1, len(lR1)-1)p2 = np.polyfit(lR2, PSNR2, len(lR1)-1)p_int1 = np.polyint(p1)p_int2 = np.polyint(p2)int1 = np.polyval(p_int1, max_int) - np.polyval(p_int1, min_int)int2 = np.polyval(p_int2, max_int) - np.polyval(p_int2, min_int)avg_diff = (int2 - int1) / (max_int - min_int)return avg_diffif __name__ == '__main__':'''# x 假设是单调的，内部会进行排序，不然就不是函数了x = np.array([1, 2, 3, 4, 5, 6])y = np.array([16., 18., 21., 17., 15., 12.])xx = np.linspace(x[0]-0.25, x[-1]+0.25, 101)interp = PCHIP(x, y, use_spline=True)spline_y = interp.interp(xx)p = np.polyfit(x, y, 5)poly_y = np.polyval(p, xx)interp2 = PCHIP(x, y, use_spline=False)pchip_y = interp2.interp(xx)plt.plot(x, y, 'ro')plt.plot(x, y, '--', linewidth=1, label='linear')plt.plot(xx, poly_y, label='Lagrange')plt.plot(xx, pchip_y, label='sppchip')plt.plot(xx, spline_y, label='spline')plt.legend()plt.show()'''rate1 = np.array([514.2076, 337.1483, 225.6044, 148.7819])psnr1 = np.array([51.2271, 47.5517, 43.9461, 40.1898])rate2 = np.array([513.5298, 337.4435, 225.9397, 148.7625])psnr2 = np.array([51.2174, 47.5473, 43.9479, 40.1754])bd_psnr = BD_PSNR(rate1, psnr1, rate2, psnr2, piecewise=True)bd_rate = BD_Rate(rate1, psnr1, rate2, psnr2, piecewise=True)print(bd_psnr, bd_rate)bd_psnr1 = BD_PSNR(rate1, psnr1, rate2, psnr2, piecewise=False)bd_rate1 = BD_Rate(rate1, psnr1, rate2, psnr2, piecewise=False)print(bd_psnr1, bd_rate1)

多项式插值之Lagrange、PCHIP与Spline以及BD-Rate和BD-PSNR的计算相关推荐

解决龙格现象matlab,matlab实现Lagrange多项式插值观察龙格现象
Matlab进行Lagrange多项式插值拉格朗日插值法对函数y=1./(1+25*x.^2)在区间[-1,1]进行5次.10次.15次插值观察龙格现象主程序 1．拉格朗日 function [c ...
数学建模准备插值（拉格朗日多项式插值，牛顿多项式插值，分段线性插值，分段三次样条插值，分段三次Hermite插值）
文章目录摘要(必看) 0 基础概念什么是插值插值用途什么是拟合插值和拟合的相同点插值和拟合的不同点 1 常用的基本插值方法 1.1 多项式插值法 1.1.1 拉格朗日多项式插值法多项式插 ...
MATLAB运用——多项式插值
文章目录实验1-误差传播与算法稳定性问题提出实验内容实验要求程序运行程序理论实验2-多项式插值的震荡和样条插值的收敛问题提出实验内容实验要求程序运行程序总结实验1-误差传 ...
插值与拟合 (一) ：拉格朗日多项式插值、Newton插值、分段线性插值、Hermite插值、样条插值、 B 样条函数插值、二维插值
插值:求过已知有限个数据点的近似函数. 拟合:已知有限个数据点,求近似函数,不要求过已知数据点,只要求在某种意义下它在这些点上的总偏差最小. 插值和拟合都是要根据一组数据构造一个函数作为近似,由于近似 ...
吉布斯现象与插值优化（上）Matlab实现多项式插值
前言傅里叶变换是贯穿一系列课程的应用,从起初的不知其所以然到学习电路理论时的轮廓逐渐清晰. 在电路理论的结课复习阶段,我在图书馆找到国外教材翻阅,关注到行文中提到的吉布斯现象(Gibbs Pheno ...
c语言构造插值多项式,拉格朗日多项式插值(C语言).docx
拉格朗日多项式插值(C语言) #include #include #include float lagrange(float *x,float *y,float xx,int n)/*拉¤-格?朗¤¨ ...
matlab 平滑曲线连接_平滑轨迹插值方法之多项式插值（附代码）
前言今天我们来聊聊轨迹插值,在机器人的运动规划和控制领域,参考轨迹的生成是一个历史悠久的问题,已经发展出了一系列的方法.今天我们就来聊一聊轨迹插值领域中最常见的轨迹插值方法:多项式插值. 说明:本文 ...
matlab全域基函数,多项式函数插值：全域多项式插值（一）单项式基插值、拉格朗日插值、牛顿插值 [MATLAB]...
全域多项式插值指的是在整个插值区域内形成一个多项式函数作为插值函数.关于多项式插值的基本知识,见"计算基本理论". 在单项式基插值和牛顿插值形成的表达式中,求该表达式在某一点处的值 ...
基于神经网络多项式插值的图像超分辨重构研究-附Matlab代码
⭕⭕ 目录 ⭕⭕ ✳️ 一.引言 ✳️ 二.基于单帧图像的超分辨率重构技术 ✳️ 2.1 最近邻域插值法 ✳️ 2.2 双线性插值法 ✳️ 2.3 双三次插值法(Keys'插值) ✳️ 三.神经网络 ...

多项式插值之Lagrange、PCHIP与Spline以及BD-Rate和BD-PSNR的计算