离散卷积定理

发表于 2024-05-15 更新于 2025-03-06 分类于 Fourier 分析阅读次数：

首先, 定义两种向量之间的运算: 乘积 \(f \circ g\) 和卷积 \(f*g\):

向量 \(f,g\) 之间的乘积计算 \(f \circ g: \mathbb{R}^n \times \mathbb{R}^n \mapsto \mathbb{R}^n\) 定义为 \[ (f \circ g)_k \triangleq f_k g_k, \quad k = 0,1,\ldots,n-1 \] 向量 \(f,g\) 之间的卷积计算 \(f*g: \mathbb{R}^n \times \mathbb{R}^n \mapsto \mathbb{R}^n\) 定义为 \[ (f*g)_k \triangleq \sum_{i=0}^{n-1} f_i g_{k-i}, \quad k = 0,1,\ldots,n-1 \] 我在这篇文章中希望说明的是, 向量的乘积/卷积运算, 总是可以通过卷积/乘积运算得到同样的结果. 将这两种运算联系到一起的就是离散卷积定理.

注:

按照上述定义, 计算卷积时会出现指标为负数的情况, 例如, 计算卷积的第一个分量时: \[ (f*g)_0 = f_0 g_0 + f_1 g_{-1} + f_2 g_{-2} + \cdots + f_{n-1} g_{-(n-1)} \] 但因为周期性 (或者称为 aliasing), \(g_{k} = g_{k+n}\), 所以上式又可以写成 \[ (f*g)_0 = f_0 g_0 + f_1 g_{n-1} + f_2 g_{n-2} + \cdots + f_{n-1} g_{1} \]

接下来我们将举一个数例, 说明卷积的计算过程

离散卷积计算

若 \(f=(1,2,3,4)\), \(g=(5,6,7,8)\), 那么 \(f,g\) 之间的卷积为

\((f*g)_0\)	\(f\)	\(g\)		\((f*g)_1\)	\(f\)	\(g\)
	1	5	5		1	6	6
	2	8	16		2	5	10
	3	7	21		3	8	24
sum = 66	4	6	24	sum = 68	4	7	28

\((f*g)_2\)	\(f\)	\(g\)		\((f*g)_3\)	\(f\)	\(g\)
	1	7	7		1	8	8
	2	6	12		2	7	14
	3	5	15		3	6	18
sum = 66	4	8	32	sum = 60	4	5	20

因此, \[ f*g = (66,68,66,60). \]

如果直接按照定义去计算卷积, 那么需要 (\(n\) 个分量 \(\times\) 每个分量需要 \(n\) 次乘法 \(=\) \(n^2\)) 计算量, 接下来我们将给出离散形式的卷积定理, 使用该方法, 并搭配上快速 Fourier 变换, 可以实现更高效率地计算卷积.

离散形式的卷积定理

若物理空间中的向量 \(f\) 和 \(g\) 分别对应相空间中的 \(c\) 和 \(d\) : \[ c \triangleq (F)^{-1} f, \quad d \triangleq (F)^{-1} g, \] 那么 \(f*g\) 等于相空间中向量之间的乘积 \(c \circ d\) 的变换, 并乘因子 \(n\) \[ f*g = n F(c \circ d) \tag{1} \]

在证明该定理之前, 我们首先通过一个算例验证一下该定理:

我们使用式 (1) 重新计算 \(f=(1,2,3,4)\), \(g=(5,6,7,8)\) 的卷积. 首先, 维度为 4 的 Fourier 变换矩阵以及逆矩阵为 \[ F = \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{pmatrix}, \quad (F)^{-1} = \frac{1}{4} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix} \] 对应的相空间中的 \(c\) 和 \(d\) 为 \[ c = (F)^{-1}f = \frac{1}{4} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix} \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \end{pmatrix} = \frac{1}{4} \begin{pmatrix} 10 \\ -2 + 2i \\ -2 \\ -2-2i \end{pmatrix} \]

\[ d = (F)^{-1}g = \frac{1}{4} \begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & -i & -1 & i \\ 1 & -1 & 1 & -1 \\ 1 & i & -1 & -i \end{pmatrix} \begin{pmatrix} 5 \\ 6 \\ 7 \\ 8 \end{pmatrix} = \frac{1}{4} \begin{pmatrix} 26 \\ -2 + 2i \\ -2 \\ -2-2i \end{pmatrix} \]

\(c\) 与 \(d\) 乘积为 \[ c \circ d = \frac{1}{16} \begin{pmatrix} 260 \\ -8i \\ 4 \\ 8i \end{pmatrix} \] 再对 \(c \circ d\) 作变换得到 \[ n F(c \circ d) = 4\begin{pmatrix} 1 & 1 & 1 & 1 \\ 1 & i & -1 & -i \\ 1 & -1 & 1 & -1 \\ 1 & -i & -1 & i \end{pmatrix} \frac{1}{16} \begin{pmatrix} 260 \\ -8i \\ 4 \\ 8i \end{pmatrix} = \begin{pmatrix} 66 \\ 68 \\ 66 \\ 60 \end{pmatrix} \]

使用离散卷积定理计算卷积的计算量为

物理空间到相空间的变换, \(2n^2\)
乘积 \(c \circ d\), \(n\)
相空间到物理空间的变换, \(n^2\)

总共需要 \(3n^2 + n\) 次乘法! 这比直接用定义计算还要多上 \(3\) 倍! 那么为什么说使用卷积定理计算速度更快呢? 这是因为在第 1 和第 3 步的矩阵乘法运算可以通过快速 Fourier 变换 (Fast Fourier Transformation, FFT) 替代, 这只需要 \(\mathcal{O}(n \log n)\) 次计算量.

接下来, 我们将用矩阵方法证明卷积定理. 首先将左侧的 \(f\) 向量拉成一个 \(n \times n\) 的矩阵, 矩阵的每一行对应上一行向右移动一个元素的位置, 右侧溢出的元素补缺到每一行的最左侧: \[ C_f \triangleq \begin{pmatrix} f_0 & f_{n-1} & \cdot & \cdot & f_{1} \\ f_{1} & f_{0} & f_{n-1} & \cdot & f_{2} \\ \cdot & f_{1} & f_{0} & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & f_{n-1} \\ f_{n-1} & f_{n-2} & \cdot & f_{1} & f_{0} \end{pmatrix} \] 可以看到, 卷积 \(f*g\) 结果就等于矩阵乘法 \(C_f g\) 得到的向量. 接下来我们考察矩阵 \(C_f\) 的特征值和特征向量. 立刻可以看到的是, 向量 \((1,1,\ldots,1)^T\) 是矩阵的一个特征向量, 并且对应的特征值为 \(\lambda_0 = \sum_{k=0}^{n-1} f_k\). 向量 \((1,w,\ldots,w^{n-1})^T\) 同样也是矩阵 \(C_f\) 的一个特征向量, 对应特征值为 \[ \lambda_1 = \sum_{k=0}^{n-1} f_{k} w^{n-k} = \sum_{k=0}^{n-1} f_{k} w^{-k} \] 我们可以验证第 \(j\) 行的计算结果, 左侧矩阵乘法得到的第 \(j\) 个分量为 \[ LHS = f_j \times 1 + f_{j-1} \times w + \cdots + f_0 w^j + f_{n-1} w^{j+1} + \ldots + f_{j+1} w^{n-1} \] 提取因子 \(w^j\) 得到 \[ LHS = w^j \left( \underbrace{f_j \times w^{-j} + f_{j-1} \times w^{-(j-1)} + \cdots + f_1 w^{-1} + f_0 w^0}_{\clubsuit} + \underbrace{f_{n-1} w^{1} + \ldots + f_{j+1} w^{n-j-1}}_{\spadesuit} \right) \] 对于式中 \(\spadesuit\) 部分, 我们可以除因子 \(w^n\) (并且不改变这一部分的大小, 因为 \(w^n=1\)) 得到 \[ f_{n-1} w^{1} + \ldots + f_{j+1} w^{n-j-1} = f_{n-1} w^{-(n-1)} + \ldots + f_{j+1} w^{-(j+1)} \] 重新排列就可以得到右端项向量第 \(j\) 行分量: \[ LHS = w^j \lambda_1 = RHS \] 一般的, 我们可以得到矩阵 \(C_f\) 的特征向量为 Fourier 变换矩阵 \(F\) 的列向量, 特征值为 \(f\) 的 Fourier 变换系数 \(c[f]\) 再乘 \(n\). 写成矩阵的形式为 \[ C_f F = F \Lambda_f,\quad \Lambda_f = n ~\mathrm{diag}~ (c[f]) \] 而这样的矩阵分解, 恰好对应于离散卷积定理. 我们现在用矩阵乘法将卷积 \(f*g\) 表示为周期矩阵 \(C_f\) 与向量 \(g\) 的矩阵乘积, \(C_f\) 可以分解为三个特殊的矩阵之间的乘积 \(F \Lambda_f (F)^{-1}\), 因此, 卷积计算过程也可以分解为对应地 “三步走”:

对向量 \(g\) 作变换 \((F)^{-1}g\), 得到关于 \(g\) 的 Fourier 变换系数 \(d\)
\(\Lambda d\) 得到向量 \(n(c \circ d)\), 也就是 \(n\) 倍的 \(c\), \(d\) 乘积
对乘积 \(n(c \circ d)\) 作变换 \(F\), 得到卷积 \(f*g\)

由此得到如下离散卷积定理的表达式: \[ (F)^{-1}(f*g) = (F)^{-1} F \Lambda_f (F)^{-1} g = n \left[ (F)^{-1}f \right] \circ \left[ (F)^{-1}g \right] \tag{2} \] 我们可以立刻看到, 式 (2) 与前文得到的式 (1) 其实说明的是一回事.

按照上面的矩阵证明方法, 同样也可以说明反方向用乘积计算卷积也是可行的. 如果我们将变换后的向量 \(c\) \(d\) 之间卷积写成矩阵乘积的形式, 就得到 \[ c*d = C_c d \] 矩阵 \(C_c\) 的特征向量为变换矩阵 \(\overline F\) 的列向量, 特征值为 \(f[c]\), 写成矩阵的形式就得到 \[ C_c \overline F = \overline F \Lambda_c,\quad \Lambda_c = \mathrm{diag}~ (f[c]) \] 所以 \[ c*d = \overline F \Lambda_c (\overline F)^{-1} d = \frac{1}{n} \overline F \Lambda_c g = F^{-1}(f \circ g) \] 综合上述讨论, 我们可以将离散卷积定理总结为如下对称且优美的形式: \[ \begin{aligned} (F)^{-1}(f*g) &= n \left[ (F)^{-1}f \right] \circ \left[ (F)^{-1}g \right] \\ F^{-1}(f \circ g) &= \left[ (F)^{-1}f \right]*\left[ (F)^{-1}f \right] \end{aligned} \]