B样条

在数学的子学科数值分析里，B-样条是样条曲线一种特殊的表示形式。它是B-样条基曲线的线性组合。B-样条是贝兹曲线的一种一般化，可以进一步推广为非均匀有理B样条（NURBS），使得我们能给更多一般的几何体建造精确的模型。

De Boor算法是一个数值上稳定的计算B样条的方法。

术语 B样条是Isaac Jacob Schoenberg创造的，B 是基（basis）样条的缩略。

## 定义

给定m+1 个节点ti，分布在0,1区间，满足

$$t_0 < t_1 < \ldots < t_m$$

一个n次B样条是一个参数曲线：

$$\mathbf{S}:[0,1] \to \mathbb{R}^2$$

它由n次B样条基(basis B-spline)组成

$\mathbf{S}(t)= \sum_{i=0}^{m} \mathbf{P}_{i} b_{i,n}(t) \mbox{ , } t \in [0,1]$.

Pi称为控制点或de Boor点.

m+1个n次B样条基可以用Cox-de Boor递归公式 定义

$$b_{j,0}(t) := \left\{\begin{matrix} 1 & \mathrm{} \quad t_j < t < t_{j+1} \\ 0 & \mathrm{...} \end{matrix} \right.$$

$$b_{j,n}(t) := \frac{t - t_{j}}{t_{j+n} - t_j} b_{j,n-1}(t) + \frac{t_{j+n+1} - t}{t_{j+n+1} - t_{j+1}} b_{j+1,n-1}(t).$$

当节点等距，称B样条为均匀(uniform)否则为非均匀(non-uniform)。

### 均匀B样条曲线

当B样条是均匀的时候，对于给定的n，每个B样条基是其他基的平移拷贝而已。一个可以作为替代的非递归定义是

$$b_{j,n}(t) := b_n(t + n - j) \qquad \mbox{ , } j = -1, \ldots m+1$$

满足

$$b_n(t) := (m+1) \sum_{i=0}^{m+1} \omega_i(t_i - t)_+^{m} \qquad \mbox{ , } t \in [0,1]$$

满足

$$\omega_i := \prod_{j=0, i \neq j}^{m+1} \frac{1}{t_i - t_k}$$

其中

$$(t_i - t)_+$$

是截断幂函数（truncated power function）

## 注解

当节点数和多项式次数相等时，B样条退化为贝兹曲线。即函数的形状由节点的位置决定。缩放或者平移节点向量不会改变基函数。

样条包含在它的控制点的凸包中

n次B样条的一个基

$$b_{i,n}(t)$$

仅当在区间ti, ti+n+1上非0。就是

$$b_{i,n}(t) = \left\{\begin{matrix} >0 & \mathrm{} \quad t_{i} \le t < t_{i+n+1} \\ 0 & \mathrm{...} \end{matrix} \right.$$

换句话说，如果我们操作一个控制点，我们只改变曲线在局部的行为，而不像Bezier曲线那样是全局行为。

## 例子

### 常数B样条

常数B样条是最简单的样条。只定义在一个节点距离上，而且不是节点的函数。它只是不同节点段（knot span）的指示函数。

$$b_{j,0}(t) = 1_{[t_j,t_{j+1})} = \left\{\begin{matrix} 1 & \mathrm{} \quad t_j \le t < t_{j+1} \\ 0 & \mathrm{...} \end{matrix} \right.$$

### 线性B样条

线性B样条定义在两个相邻的节点段上，在节点连续但不可微。

$$b_{j,1}(t) = \left\{\begin{matrix} \frac{t - t_j}{t_{j+1} - t_j} & \mathrm{if} \quad t_j \le t < t_{j+1} \\ \frac{t_{j+2} - t}{t_{j+2} - t_{j+1}} & \mathrm{ } \quad t_{j+1} \le t < t_{j+2} \\ 0 & \mathrm{... } \end{matrix} \right.$$

### 三次B样条

一个片断上的B样条的表达式可以写作：

$$S_{i} (t) = \sum_{k=0}^3 \mathbf{P}_{i-3+k} b_{i-3+k,3} (t) \qquad \mbox{ , } t \in [0,1]$$

其中Si是第i个B样条片断而P是一个控制点集，i和k是局部控制点索引。控制点的集合会是$P_i^w = ( w_i x_i, w_i y_i, w_i z_i, w_i)$的集合，其中$w_i$是比重，当它增加时曲线会被拉向控制点$P_i$，在减小时则把曲线远离该点。

片段的整个集合m-2条曲线（$S_3,S_4,...,S_m$）由m+1个控制点（$P_0,P_1,...,P_m, m \ge 3$）定义，作为t上的一个B样条可以定义为

$$S(t) = \sum_{i=0}^m \mathbf{P}_{i} b_{i,} (t)$$

其中i是控制点数，t是取节点值的全局参数。这个表达式把B样条表示为B样条基函数的线性组合，这也是这个名称的原因。

有两类B样条-均匀和非均匀。非均匀B样条相邻控制点间的距离不一定要相等。一个一般的形式是区间随着插入控制点逐步变小到0。

## B样条的程式指令

### Matlab

In Matlab，the command“spline” can be used for spline interpolation.

(Note： In the command, the cubic B-spline is used)

Cubic B-Spline Interpolation by Matlab

Generating a sine-like spline curve and samples it over a finer mesh:

x = 0:1:10; % original sampling points

y = sin(x);

xx = 0:0.1:10; % new sampling points

yy = spline(x,y,xx);

plot(x,y,'o',xx,yy)

### Python

事前安装模组

- pip install numpy
- pip install scipy
- pip install matplotlib

Cubic B-Spline Interpolation by Python

from scipy.interpolate import interp1d

import matplotlib.pyplot as plt

import numpy as np

x = np.arange(0, 11) # original sample points, 0, 1, 2, …, 9, 10

y = np.sin(x)

f = interp1d(x, y, kind=' cubic ') ) # Cubic means the cubic B-spline.

x_new = np.arange(0, 10.1, 0.1) # new sample points, 0, 0.1, 0.2, ….., 9.9, 10

y_new = f(x_new)

plt.plot(x,y,'o',x_new, y_new)

plt.show

## B样条曲线及曲面相关算法

关于此处涉及的算法，在著作中有针对Bézier、B样条（B-spline）以及非均匀有理B样条（Nurbs）的相关算法的详细数学表达和程序实现方法。

### 求导

在几何处理中，对参数曲线及曲面的求导是最基本的运算之一，由于参数表达的特性，在给定点的切线及法线可通过求导直接得到。

先来考察曲线的情形：采用本页定义中的B样条曲线表达式$$\mathbf{S}(t)= \sum_{i=0}^{m} \mathbf{P}_{i} b_{i,n}(t) \mbox{ , } t \in [0,1]$$ 对参数$t$进行求导：

$$\frac{d\mathbf{S}}{dt} = \sum_{i=0}^{m} b'_{i,n}(t) \mathbf{P}_{i}$$

## 应用

### 阵列信号处理

假设有多个讯号源，并考虑在远场 (far-field) 的条件下，以一个多天线接收机所接收，多天线接收机具有 N 个独立天线，且天线之间间距相同 𝑑 = 𝜆/2。

接收讯号$$x(t) = a(\theta_s(t))s(t) + a(\theta_i(t))i(t) + n(t), \quad t = 1, 2, \dots, L$$

其中，s(t) 为欲接收的讯号源，i(t) 为其他讯号源造成的干扰，n(t) 为噪声。

指向向量 (steering vector)为:

$$a(\theta_s(t)) = \begin{bmatrix} 1 & e^{j\pi \sin\theta} & e^{j2\pi \sin\theta} & \dots & e^{j(N-1)\pi \sin\theta} \end{bmatrix}^T$$

权重向量 (weight vector) 定义为:

$$w = \begin{bmatrix} w_0 & w_1 & \dots & w_{N-2} & w_{N-1} \end{bmatrix}^T$$

波束成型 (beamformer)的输出为:

$$y(t) = w^H x(t) = w^H \left( a(\theta_s(t))s(t) + a(\theta_i(t))i(t) + n(t) \right)$$

由于接收端不知道讯号源传送的讯号为何，利用分段平均、内插重建方法可以降低噪声 n(t) 的影响，还原低噪声的角度成分。

Step 1: 分段平均 (Segment Average)

$$\tilde{\theta}_s(\tau + 1) = \frac{1}{L} \sum_{i=0}^{L-1} \tilde{\theta}_s(L\tau + i), \; \tau \in \mathbb{N}$$

$$\tilde{\theta}_i(\tau + 1) = \frac{1}{L} \sum_{i=0}^{L-1} \tilde{\theta}_i(L\tau + i), \; \tau \in \mathbb{N}$$

Step 2: 利用 cubic B-Spline Interpolation B-Spline Interpolation 的关系式如下:

$$B_{n,0}(t) = \begin{cases} 1, & \text{for } t_n < t < t_{n+1} \\ 0, & \text{otherwise} \end{cases}$$

$$B_{n,m}(t) = \frac{t - t_n}{t_{n+m} - t_n} B_{n,m-1}(t) + \frac{t_{n+m+1} - t}{t_{n+m+1} - t_{n+1}} B_{n+1,m-1}(t) \quad \text{for } t_n < t < t_{n+1}$$

$$\text{where } B_{n,m}(t) \text{ is the } m\text{th order polynomial.}$$

因为我们是使用 cubic B-Spline interpolation，所以选择 m = 3 代入关系式。

借由上述两步骤，我们可以得到讯号源和干扰源的DOA:

$$\hat{\theta}_s(t) = \sum_{n=1}^N \tilde{\theta}_s(t_n) B_{n,3}(t)$$

$$\hat{\theta}_i(t) = \sum_{n=1}^N \tilde{\theta}_i(t_n) B_{n,3}(t)$$

优点:

用多项式内插函数重建平均化后的离散讯号可拟和出更接近连续的波型，同时降低高频成分造成的影响。

## 参考

本条目部分或全部内容出自以GFDL授权发布的《自由线上电脑词典》（FOLDOC）

## 外部链接

- Interactive java applets for B-splines
- Introduction to Computing with Geometry Notes (Dr.C.-K.Shene/Michigan Technological University)