DFT、循环卷积与差分算子

傅里叶变换通常被解释为将一个信号分解到不同频率的正弦和余弦基上。这个理解很直观，不过在处理离散信号、卷积以及差分算子时，我更习惯从线性代数的角度来看傅里叶变换：离散傅里叶变换可以写成一次矩阵乘法，而卷积、差分等很多离散线性算子，在傅里叶基下恰好具有非常简单的表示。

顺着这个思路，离散傅里叶变换、循环卷积以及差分算子其实可以自然地串在一起。首先将 DFT 写成矩阵形式 $\hat{x}=Fx$，然后把循环卷积表示为循环矩阵与向量的乘法，最后利用循环矩阵与傅里叶变换之间的关系，将卷积和差分运算转换到频域。

本文主要整理这几部分内容，同时利用 NumPy 和 SciPy 对其中的一些结论进行简单验证。

1. 离散傅里叶变换的矩阵表示

假设有一个长度为 $N$ 的离散序列

$$x=\begin{bmatrix}x_0 & x_1 & \cdots & x_{N-1}\end{bmatrix}^T.$$

用 $\hat{x}$ 表示其离散傅里叶变换。按照 DFT 的定义，第 $m$ 个频率分量可以写成

$$\hat{x}_m=\mathcal{F}(x)_m=\sum_{n=0}^{N-1}x_ne^{-\frac{2\pi mn}{N}i}. \tag{1}$$

定义

$$w=e^{-\frac{2\pi i}{N}}=\cos\frac{2\pi}{N}-i\sin\frac{2\pi}{N},$$

那么式（1）可以写成

$$\hat{x}_m=\sum_{n=0}^{N-1}x_nw^{nm}.$$

对于固定的 $m$，这实际上就是向量 $x$ 与一组复指数基的乘积：

$$\hat{x}_m=\begin{bmatrix}w^{0m} & w^{1m} & \cdots & w^{(N-1)m}\end{bmatrix}x.$$

将所有 $m=0,1,\ldots,N-1$ 的结果组合起来，可以得到

$$\hat{x}= \begin{bmatrix} w^{0\times0} & w^{1\times0} & \cdots & w^{(N-1)\times0}\\ w^{0\times1} & w^{1\times1} & \cdots & w^{(N-1)\times1}\\ \vdots & \vdots & \ddots & \vdots\\ w^{0\times(N-1)} & w^{1\times(N-1)} & \cdots & w^{(N-1)\times(N-1)} \end{bmatrix} \begin{bmatrix} x_0\\ x_1\\ \vdots\\ x_{N-1} \end{bmatrix}.$$

定义傅里叶矩阵

$$F= \begin{bmatrix} 1 & 1 & \cdots & 1\\ 1 & w & \cdots & w^{N-1}\\ 1 & w^2 & \cdots & w^{2(N-1)}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & w^{N-1} & \cdots & w^{(N-1)^2} \end{bmatrix},$$

那么 DFT 就可以非常紧凑地写成

$$\hat{x}=Fx.$$

因此，从矩阵的角度来看，离散傅里叶变换可以理解为使用复指数函数构成的基矩阵，对原始序列进行一次线性变换。

按照这里使用的傅里叶矩阵定义，可以得到

$$F^T=F,$$ $$(F^*)^T=F^*,$$ $$FF^*=NI,$$ $$FF^H=NI,$$

以及

$$F^{-1}=\frac{1}{N}F^*.$$

这里 $F^*$ 表示元素逐个取复共轭，$F^H$ 表示共轭转置。由于这里构造的傅里叶矩阵满足 $F^T=F$，因此 $F^H=(F^*)^T=F^*$。

这些性质可以直接通过 NumPy 进行验证：

import numpy as np

N = 6
x = np.array([-1, 0, 0, 0, 0, 1])

theta = np.arange(N)[:, None] * np.arange(N)[None, :]
F = np.exp(theta * (-2j * np.pi) / N)

Fx = F @ x
xhat = np.fft.fft(x)

assert np.allclose(xhat, Fx, rtol=1e-9)
assert np.allclose(np.linalg.inv(F), F.conj() / N)
assert np.allclose(F, F.T)
assert np.allclose(F @ F.conj().T, N * np.eye(N))

这里直接构造出来的 $F$ 与 np.fft.fft 所采用的 DFT 定义一致，因此 F @ x 与 np.fft.fft(x) 可以得到相同的结果。

不过，如果只是把傅里叶变换写成矩阵乘法，意义还不是特别明显。真正有意思的是，当卷积也被写成矩阵形式以后，可以看到两者之间非常直接的关系。

2. 离散傅里叶变换与循环卷积

傅里叶变换和卷积定理通常会一起出现。最常见的一句话是：

时域中的卷积，对应频域中的逐元素乘法。

对于有限长度的离散序列，更准确地说，DFT 对应的是循环卷积。

设 $x$ 和 $z$ 是两个长度均为 $N$ 的离散序列，它们的循环卷积记为 $x*z$。循环卷积定义为

$$(x*z)_m=\sum_{n=0}^{N-1}x_{(m-n)\bmod N}z_n.$$

其中 $(m-n)\bmod N$ 表示当索引超出范围时，会重新绕回序列的另一端。因此，循环卷积实际上隐含了周期性的边界处理。

循环卷积定理可以写成

$$F(x*z)=\hat{x}\odot\hat{z},$$

其中 $\odot$ 表示逐元素乘法。

在直接讨论这个结论之前，可以先把循环卷积写成矩阵形式。

假设 $N=4$，那么循环卷积可以展开成

$$x*z= \begin{bmatrix} x_0 & x_3 & x_2 & x_1\\ x_1 & x_0 & x_3 & x_2\\ x_2 & x_1 & x_0 & x_3\\ x_3 & x_2 & x_1 & x_0 \end{bmatrix} \begin{bmatrix} z_0\\ z_1\\ z_2\\ z_3 \end{bmatrix}.$$

可以看到，卷积实际上对应一个具有循环结构的矩阵。

对于序列 $x$，定义循环矩阵

$$C(x)= \begin{bmatrix} x_0 & x_1 & x_2 & \cdots & x_{N-1}\\ x_{N-1} & x_0 & x_1 & \cdots & x_{N-2}\\ x_{N-2} & x_{N-1} & x_0 & \cdots & x_{N-3}\\ \vdots & \vdots & \ddots & \vdots & \vdots\\ x_1 & x_2 & \cdots & x_{N-1} & x_0 \end{bmatrix}.$$

按照这里采用的定义，循环卷积可以写成

$$x*z=C^T(x)z.$$

这样一来，卷积就从一个“不断平移序列并求和”的过程变成了普通的矩阵乘法问题。

3. 循环矩阵与傅里叶变换

循环矩阵和傅里叶变换之间存在非常紧密的联系。原来的推导可以从一个循环移位置换矩阵开始。

定义

$$P= \begin{bmatrix} 0 & 0 & 0 & \cdots & 1\\ 1 & 0 & 0 & \cdots & 0\\ 0 & 1 & 0 & \cdots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 & 0 \end{bmatrix}.$$

矩阵 $P$ 作用于序列时，相当于将序列循环移动一个位置。利用 $w^N=1$，可以写成

$$P=F\operatorname{diag}(1,w^{-1},w^{-2},\ldots,w^{-(N-1)})F^{-1}.$$

进一步有

$$P^n=F\operatorname{diag}(1,w^{-n},w^{-2n},\ldots,w^{-(N-1)n})F^{-1}.$$

另一方面，循环矩阵可以通过不同次数的 $P$ 组合得到：

$$C^T(x)=x_0I+x_1P+x_2P^2+\cdots+x_{N-1}P^{N-1}.$$

将 $P$ 的傅里叶表示代入，可以得到

$$\begin{aligned} C^T(x) &=F\operatorname{diag}\left( \sum_{n=0}^{N-1}x_n, \sum_{n=0}^{N-1}x_nw^{-n}, \sum_{n=0}^{N-1}x_nw^{-2n}, \ldots, \sum_{n=0}^{N-1}x_nw^{-(N-1)n} \right)F^{-1}. \end{aligned}$$

也就是说，循环矩阵经过傅里叶变换以后，会转换为对角矩阵。按照前面的符号，可以写成

$$C(x)=F\operatorname{diag}(\hat{x})F^{-1},$$

以及

$$C^T(x)=F\operatorname{diag}(\hat{x}^*)F^{-1}.$$

由循环卷积 $x*z=C^T(x)z$，可以进一步得到

$$C^T(x)z=F^{-1}(\hat{x}\odot\hat{z}),$$

也就是

$$x*z=F^{-1}(\hat{x}\odot\hat{z}).$$

两边进行傅里叶变换，最终得到

$$F(x*z)=\hat{x}\odot\hat{z}.$$

这就是前面提到的循环卷积定理。

从这个角度来看，“卷积在频域中变成乘法”并不是一个孤立的公式。其背后真正发生的事情，是循环卷积所对应的循环矩阵在傅里叶基下变成了对角形式。

原本需要进行矩阵与向量乘法的运算，在傅里叶域中只需要对对应的频率分量逐元素相乘。

4. 循环矩阵的一些性质

有了傅里叶表示以后，循环矩阵的一些性质也比较容易整理出来。例如两个循环矩阵满足

$$C(x)C(z)=C(z)C(x),$$

以及线性关系

$$\alpha C(x)+\beta C(z)=C(\alpha x+\beta z).$$

另外，

$$C(x)C(z)=C\left(F^{-1}(\hat{x}\odot\hat{z})\right).$$

如果 $C(x)$ 非奇异，那么还可以写成

$$C^{-1}(x)=C\left(F^{-1}(1/\hat{x})\right).$$

这些关系背后的原因基本相同：循环矩阵经过傅里叶变换以后，对应的是一个对角矩阵，而对角矩阵的乘法、求逆等操作都可以直接通过其对角元素完成。

因此，很多原本表现为矩阵运算的问题，到了傅里叶域以后都会转化为更加简单的逐元素运算。

5. 用代码验证循环卷积

循环卷积在实际代码中的一个重要问题是卷积核的位置。

考虑一个长度为 3 的卷积核：

x = np.array([1, 2, 3])

如果使用 scipy.ndimage.convolve，卷积核的中心位于中间位置。对于长度为 3 的卷积核，其半径为 r = len(x) // 2。

而使用 FFT 实现循环卷积时，需要先将卷积核补零到与信号相同的长度，并通过循环移动将卷积核的中心移动到数组原点。

对于周期边界，可以写成：

import numpy as np
from scipy.ndimage import convolve

ks = 3
N = 5
r = ks // 2

x = np.array([1, 2, 3])
y = np.random.randn(N)

z = convolve(y, x, mode="wrap")

x_periodic = np.roll(np.pad(x, (0, N - ks)), shift=-r)
xhat = np.fft.fft(x_periodic)
yhat = np.fft.fft(y)

z_ = np.fft.ifft(xhat * yhat)

assert np.allclose(z_, z, rtol=1e-9)

mode="wrap" 使用循环边界，因此能够和 DFT 所对应的循环卷积直接对应。

如果采用 mode="constant"，则需要先对信号进行补零，将问题嵌入一个更大的序列中。例如：

ks = 3
N = 5
r = ks // 2
L = N + 2 * r

x = np.array([1, 2, 3])
y = np.random.randn(N)

z = convolve(y, x, mode="constant")

x_periodic = np.roll(np.pad(x, (0, L - ks)), shift=-r)
y_periodic = np.pad(y, (r, r))

xhat = np.fft.fft(x_periodic)
yhat = np.fft.fft(y_periodic)

z_ = np.fft.ifft(xhat * yhat)[r:r + N]

assert np.allclose(z_, z, rtol=1e-9)

这里实际做的是先扩展计算区域，再利用扩展区域中的循环卷积模拟原始边界条件下的卷积。

因此，在实际利用傅里叶变换处理卷积时，除了卷积定理本身，还需要特别注意两个问题：边界条件是否一致，以及卷积核中心在离散序列中的位置是否正确。

如果这两个地方处理不一致，即使频域中的逐元素乘法没有任何问题，最终结果也可能出现平移或者边界上的差异。

6. 差分算子也可以写成卷积

理解卷积以后，差分算子的傅里叶表示就比较自然了。

对于一维离散序列，一阶差分可以通过卷积核

$$[1,-1]$$

实现。二阶差分则可以通过

$$[1,-2,1]$$

实现，因此可以写成

$$\nabla^2x=[1,-2,1]*x.$$

在周期边界条件下，对应的矩阵形式为

$$\nabla^2x= \begin{bmatrix} -2 & 1 & 0 & \cdots & 1\\ 1 & -2 & 1 & \cdots & 0\\ 0 & 1 & -2 & \cdots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 1 & 0 & \cdots & 1 & -2 \end{bmatrix}x.$$

可以看到，这仍然是一个循环矩阵，因此可以直接使用前面的傅里叶分析方法。

按照循环卷积中的卷积核排列方式，与二阶差分对应的周期序列可以写成

$$\begin{bmatrix}-2 & 1 & 0 & \cdots & 0 & 1\end{bmatrix}^T.$$

对这个序列进行 DFT，第 $k$ 个频率分量为

$$w^k+w^{-k}-2.$$

由于 $w^k=e^{-2\pi ik/N}$，有

$$w^k+w^{-k}=2\cos\frac{2\pi k}{N}.$$

因此，二阶差分算子在傅里叶域中的表示为

$$\lambda_k=2\cos\frac{2\pi k}{N}-2.$$

于是空间域中的二阶差分运算可以写成

$$\nabla^2x=F^{-1}(\lambda\odot Fx).$$

也就是说，一维二阶差分在空间域中需要访问相邻位置，而在傅里叶域中只需要将每一个频率分量乘以对应的 $\lambda_k$。

可以用下面的代码进行验证：

import numpy as np
from scipy.ndimage import convolve

N = 7
x = np.random.randn(N)

kernel = np.array([1, -2, 1])
A = 2 * np.cos(2 * np.pi * np.arange(N) / N) - 2

z = convolve(x, kernel, mode="wrap")
z_ = np.fft.ifft(np.fft.fft(x) * A)

assert np.allclose(z, z_, rtol=1e-9)

这里的 A 就是差分算子在每一个离散频率位置上的响应。

7. 二维差分算子的频域形式

同样的思路可以自然推广到二维。

对于二维数据，常见的四邻域 Laplacian 卷积核为

$$K= \begin{bmatrix} 0 & 1 & 0\\ 1 & -4 & 1\\ 0 & 1 & 0 \end{bmatrix}.$$

假设输入数据大小为 $H\times W$。在周期边界条件下，其傅里叶域表示可以写成

$$A_{p,q}=2\left(\cos\frac{2\pi p}{H}+\cos\frac{2\pi q}{W}\right)-4.$$

因此二维 Laplacian 可以表示为

$$\nabla^2x=\mathcal{F}^{-1}(A\odot\mathcal{F}(x)).$$

代码验证如下：

import numpy as np
from scipy.ndimage import convolve

H, W = 10, 10
kernel = np.array([[0., 1., 0.], [1., -4., 1.], [0., 1., 0.]])
x = np.random.randn(H, W)

H_ind = np.arange(H)[:, None]
W_ind = np.arange(W)[None, :]
A = 2 * (np.cos(2 * np.pi * H_ind / H) + np.cos(2 * np.pi * W_ind / W)) - 4

z = convolve(x, kernel, mode="wrap")
z_ = np.fft.ifft2(np.fft.fft2(x) * A)

assert np.allclose(z, z_, rtol=1e-5)

二维情况下与一维情况并没有本质区别，只是原来的一维频率索引 $k$ 变成了二维索引 $(p,q)$。每增加一个空间维度，只需要把对应方向上的频率响应继续加入即可。

8. 三维差分算子的频域形式

对于三维数据，使用六邻域 Laplacian 时，卷积核中心位置为 $-6$，并在前后、上下、左右六个相邻位置取值 1。

假设三维数据尺寸为 $H\times W\times D$，那么对应的频率响应为

$$A_{p,q,r}=2\left(\cos\frac{2\pi p}{H}+\cos\frac{2\pi q}{W}+\cos\frac{2\pi r}{D}\right)-6.$$

因此三维 Laplacian 同样可以写成

$$\nabla^2x=\mathcal{F}^{-1}(A\odot\mathcal{F}(x)).$$

代码验证如下：

import numpy as np
from scipy.ndimage import convolve

H, W, D = 30, 100, 50

kernel = np.zeros((3, 3, 3))
kernel[1, 1, 1] = -6
kernel[0, 1, 1] = kernel[2, 1, 1] = 1
kernel[1, 0, 1] = kernel[1, 2, 1] = 1
kernel[1, 1, 0] = kernel[1, 1, 2] = 1

x = np.random.randn(H, W, D)

H_ind = np.arange(H)[:, None, None]
W_ind = np.arange(W)[None, :, None]
D_ind = np.arange(D)[None, None, :]
A = 2 * (np.cos(2 * np.pi * H_ind / H) + np.cos(2 * np.pi * W_ind / W) + np.cos(2 * np.pi * D_ind / D)) - 6

z = convolve(x, kernel, mode="wrap")
z_ = np.fft.ifftn(np.fft.fftn(x) * A)

assert np.allclose(z, z_, rtol=1e-5)

从一维到二维再到三维，整个推导过程实际上没有发生本质变化。空间域中的差分算子仍然对应一个卷积，而周期卷积具有循环结构，因此可以通过傅里叶变换将其转换成逐频率的乘法。

9. 从 DFT 再回头看卷积和差分

回过头来看，本文实际上只用了一个非常简单的起点：

$$\hat{x}=Fx.$$

但是将傅里叶变换放在线性代数框架下以后，可以逐步得到很多有用的结果。

首先，循环卷积可以写成循环矩阵与向量的乘法：

$$x*z=C^T(x)z.$$

循环矩阵可以通过傅里叶变换转化成对角形式，因此循环卷积满足

$$F(x*z)=\hat{x}\odot\hat{z}.$$

差分算子本身又只是特殊的卷积核，因此同样可以转换到傅里叶域。一维二阶差分对应

$$A_k=2\cos\frac{2\pi k}{N}-2,$$

二维四邻域 Laplacian 对应

$$A_{p,q}=2\left(\cos\frac{2\pi p}{H}+\cos\frac{2\pi q}{W}\right)-4,$$

三维六邻域 Laplacian 则对应

$$A_{p,q,r}=2\left(\cos\frac{2\pi p}{H}+\cos\frac{2\pi q}{W}+\cos\frac{2\pi r}{D}\right)-6.$$

因此，在空间域中表现为局部邻域计算的卷积和差分运算，进入傅里叶域以后，都可以转化成对应频率分量上的逐元素乘法。

这也是我觉得从矩阵角度理解傅里叶变换比较有意思的地方。很多看起来彼此独立的公式，其实都可以归结到同一个结构上：DFT 提供了一个新的基底，而循环卷积和周期差分算子在这个基底下具有更加简单的表示。

真正理解这一点以后，“时域卷积对应频域乘法”就不再只是一条需要记忆的结论，而是能够从循环矩阵与傅里叶变换之间的关系自然推导出来。

在实际计算中，还需要特别留意边界条件和卷积核位置。DFT 天然对应周期性的离散序列，因此 mode="wrap" 可以直接与循环卷积对应；如果使用常数填充或其他边界处理，则需要先进行相应的扩展和对齐，之后才能正确地利用 FFT 实现对应运算。