4.1.1 频率的直觉#

想象你站在海边，看一波波浪花打上来：

• 低频浪花：起伏缓慢，对应图像中的平坦区域 / 整体亮度
• 高频浪花：起伏密集，对应图像中的边缘 / 纹理 / 噪声

傅里叶定理（一句话）：任何信号都可以分解成无数个”不同频率、不同相位、不同幅度的正弦波”的叠加。

对图像而言： [ \text{图像} = \sum_{\text{频率}} \text{某频率的正弦波} \times \text{对应系数} ]

频域处理的本质：

1. 把图像”拆”成一堆正弦波的系数
2. 修改某些系数（比如削弱高频）
3. 把它们”拼”回图像

4.1.2 频域的工程价值#

4.1.3 频谱的”几何语言”#

一张图的频谱（DFT 结果）是一张复数图。通常画幅度谱 (|F(u, v)|)：

• 中心：低频（DC 分量 = 全图均值）
• 向外：频率升高
• 水平亮线：图中有垂直边缘/纹理
• 垂直亮线：图中有水平边缘/纹理
• 对角亮线：图中有对角纹理

1
幅度谱示意:
2
┌─────────────┐
3
│             │
4
│    ·        │   ← 高频边缘噪声散布
5
│  ·  ·       │
6
│  · ★ ·      │   ← 中心亮点 = DC 分量（全图平均）
7
│  ·  ·       │
8
│    ·        │
9
│             │
10
└─────────────┘

看懂频谱是重要基本功。

4.2.1 连续傅里叶变换（CFT）#

定义 4.1： [ F(u) = \mathcal{F}{f(t)} = \int_{-\infty}^{+\infty} f(t) e^{-j2\pi u t} dt ] 逆变换： [ f(t) = \mathcal{F}^{-1}{F(u)} = \int_{-\infty}^{+\infty} F(u) e^{j2\pi u t} du ]

其中 (u) 的单位是”循环/秒”或”循环/米”，不是角频率（注意：有的教材用 (\omega = 2\pi u)）。

4.2.2 欧拉公式回顾#

[ e^{j\theta} = \cos\theta + j\sin\theta ]

因此 (F(u)) 的实部对应余弦成分，虚部对应正弦成分。分解复振荡的利器。

4.2.3 离散傅里叶变换（DFT）#

对长度 (M) 的序列 (f(0), f(1), \ldots, f(M-1))： [ F(u) = \sum_{x=0}^{M-1} f(x) e^{-j 2\pi u x / M}, \quad u = 0, 1, \ldots, M-1 ] [ f(x) = \frac{1}{M} \sum_{u=0}^{M-1} F(u) e^{j 2\pi u x / M} ]

基函数： [ \phi_u(x) = e^{-j 2\pi u x / M}, \quad u = 0, 1, \ldots, M-1 ] 这 (M) 个基向量正交： [ \sum_{x=0}^{M-1} \phi_u(x) \overline{\phi_v(x)} = M \delta(u - v) ]

DFT 可以看作 (\mathbb{C}^M) 上一组标准正交基的坐标变换。

4.2.4 快速傅里叶变换（FFT）#

DFT 的朴素实现 (O(M^2))。FFT（Cooley-Tukey 1965）利用”分治 + 对称性”降到 (O(M \log M))。

加速比：

没有 FFT，图像处理就没有频域一章。

FFT 尺寸建议：2 的幂最快，其他组合使用 “Bluestein 变换” 或补零到最近 2 的幂。cv2.getOptimalDFTSize(N) 给出推荐尺寸。

4.3.1 定义#

定义 4.2（2D-DFT）：对 (M \times N) 图像 (f(x, y))： [ F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) e^{-j 2\pi (ux/M + vy/N)} ]

逆变换： [ f(x, y) = \frac{1}{MN} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{j 2\pi (ux/M + vy/N)} ]

4.3.2 可分离性（重要）#

指数函数可分离： [ e^{-j2\pi(ux/M + vy/N)} = e^{-j2\pi ux/M} \cdot e^{-j2\pi vy/N} ]

因此 2D-DFT 可以先对每行做 1D-DFT，再对每列做 1D-DFT： [ F(u, v) = \sum_{y=0}^{N-1} \left[ \sum_{x=0}^{M-1} f(x, y) e^{-j 2\pi u x / M} \right] e^{-j 2\pi v y / N} ]

复杂度：(M \cdot N \cdot (M + N)) → 用 FFT 后为 (O(MN \log(MN)))。

4.3.3 频率的物理意义#

对一幅 (M \times N) 图：

• (u = 0)：水平方向常数（DC）
• (u = M/2)：水平方向最高频（每 2 像素振荡一次 = 奈奎斯特频率）
• (v) 同理

周期性： [ F(u, v) = F(u + M, v + N) ]

所以 (u \in [0, M-1]) 其实覆盖了频率范围 ([-M/2, M/2))。

4.3.4 频谱中心化（必会操作）#

DFT 直接输出的频谱：

• 低频在四角
• 高频在中间

这对人眼不直观。解决办法：

方法 1：fftshift

1
Fshift = np.fft.fftshift(F)  # 把左上角的 DC 移到中心

方法 2：预乘 ((-1)^{x+y})（冈萨雷斯教材做法）

利用 DFT 平移性质： [ f(x, y)(-1)^{x+y} \leftrightarrow F(u - M/2, v - N/2) ]

在空间域乘 ((-1)^{x+y})（即棋盘格式的±1），等价于在频域把频谱中心移到 ((M/2, N/2))。

两种方法等价，fftshift 简单得多。

4.3.5 频谱可视化#

1
F = np.fft.fft2(f)
2
Fshift = np.fft.fftshift(F)
3
magnitude = np.log1p(np.abs(Fshift))   # 对数压缩
4
spectrum = cv2.normalize(magnitude, None, 0, 255, cv2.NORM_MINMAX).astype(np.uint8)

为什么要取对数？

• 频谱动态范围大（DC 系数 (\sim 10^8)，高频 (\sim 10^{-2})）
• 直接显示只看到中心亮点
• 取对数压缩后可见细节（回顾 3.2.3）

4.3.6 幅度 vs 相位（反直觉但重要）#

把 (F(u, v)) 写成极形式： [ F(u, v) = |F(u, v)| e^{j \phi(u, v)} ]

• 幅度谱 (|F|)：各频率有多强
• 相位谱 (\phi)：各频率如何对齐

Oppenheim 经典实验：用图像 A 的幅度 + 图像 B 的相位重构，结果看起来更像 B——

1
|A| + φ(B) → 看起来像 B
2
|B| + φ(A) → 看起来像 A

结论：相位携带图像结构（物体位置/轮廓），幅度只管能量。

这个性质在相位相关（Phase Correlation）中用来做平移估计：两图只差一个平移，幅度谱不变，相位差含有平移量。

4.4.1 卷积定理的严格证明（一维）#

定理 4.1：对周期延拓的信号， [ \mathcal{F}{f * h} = F \cdot H ]

证明： [ (f * h)(x) = \sum_m f(m) h(x - m) ] 代入 DFT： [ \begin{aligned} \mathcal{F}{f * h}(u) &= \sum_x \left[\sum_m f(m) h(x-m)\right] e^{-j2\pi ux/M} \ &= \sum_m f(m) \sum_x h(x-m) e^{-j2\pi ux/M} \end{aligned} ] 令 (k = x - m)： [ = \sum_m f(m) e^{-j2\pi um/M} \sum_k h(k) e^{-j2\pi uk/M} = F(u) H(u) ]

二维情形类似。证毕。

4.4.2 卷积定理的工程意义#

1
空间域             频率域
2
  f                  F
3
  *   ≡            ·
4
  h                  H
5
  ▼                  ▼
6
 f*h                F·H

三种用法：

1. 加速大核卷积：(k > 30) 时，FFT 方案更快
2. 设计滤波器：直接在频域画 (H)，反变换即得空间核
3. 分析图像：看 (F) 就知道 (f) 的频率成分

4.4.3 循环卷积 vs 线性卷积（坑点）#

DFT 的周期性意味着：DFT 对应的卷积是循环卷积！

循环卷积：把信号首尾相接，周期延拓。 线性卷积（空间域真实做的）：无延拓，边界外为 0。

两者结果不等——循环卷积会让滤波器从一边”卷”到另一边，产生边缘错位。

解决方案：零填充

要用 FFT 计算线性卷积 (f \in \mathbb{R}^M, h \in \mathbb{R}^k)：

1. 对 (f, h) 各自补零到长度 (P \ge M + k - 1)
2. 各自做 DFT
3. 点乘
4. 逆 DFT
5. 取前 (M + k - 1) 个样本

二维情形同理：补零到 ((M + m - 1) \times (N + n - 1))。

1
# 正确做法
2
P = cv2.getOptimalDFTSize(M + m - 1)  # 优化尺寸
3
Q = cv2.getOptimalDFTSize(N + n - 1)
4
f_pad = np.pad(f, ((0, P-M), (0, Q-N)), mode='constant')
5
h_pad = np.pad(h, ((0, P-m), (0, Q-n)), mode='constant')
6
result = np.real(np.fft.ifft2(np.fft.fft2(f_pad) * np.fft.fft2(h_pad)))[:M, :N]

4.5.1 七步法#

1
(1) 输入图像 f(x, y)       大小 M × N
2
                 │
3
                 ▼ 零填充
4
(2) f_pad           大小 P × Q  (P ≥ 2M, Q ≥ 2N)
5
                 │
6
                 ▼ 预乘 (-1)^(x+y) 或后续 fftshift
7
(3)
8
                 │
9
                 ▼ DFT
10
(4) F(u, v)
11
                 │  ×
12
                 │
13
(5) 滤波器 H(u, v)       (需与 F 等尺寸)
14
                 │
15
                 ▼
16
(6) G = F · H
17
                 │
18
                 ▼ IDFT, 取实部
19
(7) g'
20
                 │
21
                 ▼ 反预乘 / ifftshift
22
(8) 裁剪到 M × N → g

4.5.2 设计滤波器的两种路线#

1. 频域直接画：设计 (H(u, v)) 的形状（低通、高通、带阻、陷波），反变换给出空间核
2. 空间域核 → 频域：已知 (h(x, y))，H = FFT(h)

两种方法在工程中都常用。

4.6.1 距离函数#

所有低通滤波器都用到”到频域中心的距离”： [ D(u, v) = \sqrt{(u - M/2)^2 + (v - N/2)^2} ]

截止频率 (D_0)：决定低通的”宽度”。(D_0) 越小，保留的频率越少，图像越模糊。

4.6.2 理想低通 ILPF#

[ H_{\text{ILPF}}(u, v) = \begin{cases} 1, & D(u, v) \le D_0 \ 0, & D(u, v) > D_0 \end{cases}

1
**频域图像**：一个半径 \(D_0\) 的实心圆盘，内 1 外 0。
2

3
#### 空间域核（反变换）
4

5
\[
6
h(x, y) \propto \frac{J_1(2\pi D_0 r)}{r}, \quad r = \sqrt{x^2 + y^2}
7
\]
8
（\(J_1\) 为 1 阶贝塞尔函数，对应一维 `sinc`）
9

10
这是一个**振荡、无限支集**的函数！
11

12
#### 严重缺陷：振铃（Gibbs 现象）
13

14
ILPF 的空间核有无穷振荡 → 卷积后边缘附近产生"波纹"：

原边缘: ILPF 滤波后: ───── ───── │ ╲ ← 过冲 (overshoot) └──── └──↘ ╲↗↘ ← 振铃 ╲↗

1
**理论原因**：频域剧烈不连续（1→0 的阶跃）⇄ 空间域无限延伸。
2

3
ILPF 不实用，但是教学上重要。
4

5
### 4.6.3 Butterworth 低通 BLPF
6

7
\[
8
H_{\text{BLPF}}(u, v) = \frac{1}{1 + [D(u, v) / D_0]^{2n}}
9
\]
10

11
**参数 \(n\)**（阶数）：
12
- \(n = 1\)：过渡平缓，无振铃，但保留高频多
13
- \(n = 2\)：工程常用
14
- \(n \to \infty\)：退化为 ILPF（出现振铃）
15

16
\(D = D_0\) 处 \(H = 0.5\)，恰是"半功率"点。
17

18
### 4.6.4 高斯低通 GLPF
19

20
\[
21
H_{\text{GLPF}}(u, v) = e^{-D^2(u, v) / (2 D_0^2)}
22
\]
23

24
**优点**：
25
- 高斯的 FT 还是高斯 → **无振铃**
26
- 频域空间都**紧致平滑**，工程中最推荐
27

28
**GLPF 与空间域 Gaussian Blur 的关系**：两者等价，只差一个 \(\sigma\) 的对应。
29

30
### 4.6.5 三种低通的对比
31

32
| 滤波器 | 振铃 | 过渡带 | 推荐度 |
33
|--------|------|-------|--------|
34
| ILPF | 严重 | 无 | 理论分析 |
35
| BLPF n=2 | 轻微 | 平滑 | 工程折中 |
36
| GLPF | 无 | 平滑 | **实际首选** |

三者频率响应 (D0 相同):

H(D) 1 ┤───── (ILPF) │ │ │ │╲ (BLPF) ½ ┤ │ ╲╲ │ │ ╲_ │ ╲ ╲╲╲__ (GLPF) 0 ┼─────┴────────→ D D0

1
---
2

3
## 4.7 频域高通滤波（锐化）
4

5
### 4.7.1 高通由低通派生

H_HP(u, v) = 1 - H_LP(u, v)

1
于是：
2
- **IHPF** = 1 - ILPF
3
- **BHPF** = 1 - BLPF
4
- **GHPF** = 1 - GLPF
5

6
**直觉**：低通保留低频，1 - 低通 保留高频。
7

8
**注意**：高通滤波会**去除 DC 分量**（\(H(0, 0) = 0\)）→ 输出图像均值为 0，负值多，需要加一个偏置或取绝对值显示。
9

10
### 4.7.2 拉普拉斯的频域表示
11

12
拉普拉斯算子 \(\nabla^2\)，其频域响应为：
13
\[
14
\mathcal{F}\{\nabla^2 f\} = -4\pi^2 (u^2 + v^2) F(u, v)
15
\]
16

17
中心化后：
18
\[
19
H_{\text{Laplacian}}(u, v) = -4\pi^2 [(u - M/2)^2 + (v - N/2)^2]
20
\]
21

22
这是一个**二次曲面**，中心 0，越往外越负（绝对值越大）。
23

24
**锐化公式**（频域等价于 \(g = f - c \nabla^2 f\)）：
25
\[
26
G(u, v) = F(u, v) - c \cdot H_{\text{Lap}}(u, v) F(u, v) = [1 - c H_{\text{Lap}}] F
27
\]
28

29
归一化后即一个高通滤波器。
30

31
### 4.7.3 High-Boost 滤波（频域版）
32

33
High-Boost = 强化版的 Unsharp Masking：
34
\[
35
H_{\text{HB}}(u, v) = (A - 1) + H_{\text{HP}}(u, v), \quad A \ge 1
36
\]
37

38
- \(A = 1\)：普通高通
39
- \(A > 1\)：既保留原图能量又增强高频（= 锐化但不过度）
40

41
---
42

43
## 4.8 同态滤波 Homomorphic Filtering
44

45
### 4.8.1 动机
46

47
回忆 1.1.1：图像 = 照度 × 反射率
48
\[
49
f(x, y) = i(x, y) \cdot r(x, y)
50
\]
51
- \(i\)：缓变（低频），不均匀光照的元凶
52
- \(r\)：反映细节（高频）
53

54
若场景光照不均（一半亮一半暗），想**抑制光照变化、保留物体细节**，就要**分别处理 \(i\) 和 \(r\)**。
55

56
### 4.8.2 关键技巧：取对数
57

58
\(f = i \cdot r\)（非线性，乘法不能直接 FFT）。取对数：
59
\[
60
\ln f = \ln i + \ln r
61
\]
62

63
现在是**加性**！可以：
64
1. 对两边 FFT：\(\mathcal{F}\{\ln f\} = \mathcal{F}\{\ln i\} + \mathcal{F}\{\ln r\}\)
65
2. 设计一个滤波器 \(H(u, v)\) **压低低频、提升高频**
66
3. IFFT，再 `exp`
67

68
### 4.8.3 典型同态滤波器
69

70
\[
71
H(u, v) = (\gamma_H - \gamma_L) \left[1 - e^{-c D^2(u,v) / D_0^2}\right] + \gamma_L
72
\]
73

74
- \(\gamma_L < 1\)：低频衰减倍数（压暗部光照）
75
- \(\gamma_H > 1\)：高频增益（提升细节对比度）
76
- \(c\)：过渡斜率

H(D) γ_H ┤─────── │ ╲_ │ ╲_ ← 过渡区 │ ╲___ γ_L ┤───────────── 0 D0 D

1
### 4.8.4 完整流程

f → ln(1+f) → FFT → H(u,v) 相乘 → IFFT → exp - 1 → 归一化 → g

1
### 4.8.5 应用
2

3
- 医学 X 光**背光不均**校正
4
- 卫星图光照不均校正
5
- 低光照手机拍照预处理
6

7
---
8

9
## 4.9 选择性滤波：带阻、带通、陷波
10

11
### 4.9.1 带阻 Bandreject
12

13
去除某**环形**频率带：
14
\[
15
H_{\text{BR}}(u, v) = \frac{1}{1 + \left[\frac{D W}{D^2 - D_0^2}\right]^{2n}}
16
\]
17
其中 \(D_0\) 中心半径，\(W\) 带宽。
18

19
**用途**：消除具有**特定频率**的周期噪声（如旧电视的行干扰）。
20

21
### 4.9.2 带通 Bandpass
22

23
\[
24
H_{\text{BP}} = 1 - H_{\text{BR}}
25
\]
26

27
提取特定尺度的纹理。
28

29
### 4.9.3 陷波 Notch（必会）
30

31
更精确——去掉频谱上**某几个对称点**（而非整个环）。
32

33
设要抑制的噪声点为 \((u_k, v_k)\)（频谱上的亮斑），由 DFT 共轭对称性，同时存在对称点 \((-u_k, -v_k)\)。两个点都要抑制。
34

35
Butterworth 陷波：
36
\[
37
H(u, v) = \prod_{k=1}^{Q} \frac{1}{1 + [D_0 / D_k(u, v)]^{2n}} \cdot \frac{1}{1 + [D_0 / D_{-k}(u, v)]^{2n}}
38
\]
39
其中
40
\[
41
D_k(u, v) = \sqrt{(u - M/2 - u_k)^2 + (v - N/2 - v_k)^2}
42
\]
43

44
**直觉**：\(D_k\) 是到陷波中心的距离。当 \(D_k \to 0\)（靠近要抑制的点），\(H \to 0\)；远离则 \(H \to 1\)。
45

46
### 4.9.4 陷波使用流程（真实工作流）

1. 读取含周期干扰的图 f
2. 计算 F = fftshift(fft2(f))，显示 log|F|
3. 肉眼找出噪声对应的亮斑对 (u_k, v_k)
4. 构造 notch reject 滤波器 H
5. G = F · H
6. g = real(ifft2(ifftshift(G)))
7. 对比恢复效果

1
这是图像处理中极实用的技能（扫描文档去条纹、遥感图像去条带等）。
2

3
### 4.9.5 最优陷波（Gonzalez 推荐）
4

5
实际工程中，噪声估计 \(\hat{\eta}\) 与原图有相关，直接减去会残留。**最优陷波**：
6
\[
7
\hat{f}(x, y) = g(x, y) - w(x, y) \hat{\eta}(x, y)
8
\]
9
\(w\) 通过**最小化局部方差**选取。复杂但效果最好。
10

11
---
12

13
## 4.10 FFT 的实现与优化
14

15
### 4.10.1 主流库
16

17
```python
18
# NumPy
19
import numpy as np
20
F = np.fft.fft2(f)             # 完整复数输出 (实数可用 rfft2 省一半)
21
F_shift = np.fft.fftshift(F)
22

23
# OpenCV
24
import cv2
25
F = cv2.dft(f.astype(np.float32), flags=cv2.DFT_COMPLEX_OUTPUT)
26
# cv2 的 DFT 比 np 快，但接口稍复杂
27

28
# SciPy
29
from scipy.fft import fft2, ifft2, fftshift

4.10.2 实信号优化：RFFT#

图像是实数。实信号的 DFT 具有共轭对称性： [ F(-u, -v) = F^*(u, v) ]

只需存储一半即可。numpy.fft.rfft2 输出尺寸 ((M, N/2 + 1))，省一半内存、一半时间。

4.10.3 尺寸优化#

FFT 在 “2 的幂” 上最快。现代库支持 (2^a 3^b 5^c) 型尺寸：

1
N_opt = cv2.getOptimalDFTSize(N)

一张 1023×1023 图 → 补零到 1024×1024，速度快几倍。

4.10.4 GPU FFT#

• cuFFT（CUDA）
• OpenCV CUDA 模块：cv2.cuda.dft
• PyTorch torch.fft

大图像的实时处理必备。

✅ 必须掌握#

• 2D-DFT 定义与可分离实现
• 频谱中心化与对数可视化
• 幅度与相位的分工
• 卷积定理及其证明
• 循环卷积 vs 线性卷积，零填充的必要性
• ILPF/BLPF/GLPF 的振铃对比及其数学原因
• 由低通派生高通
• 同态滤波流程（log → FFT → H → IFFT → exp）
• 陷波滤波去周期噪声

💡 高频面试题#

Q1. 卷积定理的物理意义？

答：空间域的滤波（卷积）等价于频率域的每个频率分量乘以滤波器在该频率的响应。这意味着：

1. 空间核越简单，频域响应越平滑（反之亦然）
2. 大核卷积可以用 FFT 加速到 (O(N\log N))
3. 设计滤波器可以直接在频率域画图

Q2. 为什么直接在频域乘后反变换会有环绕伪影？

答：DFT 隐含周期延拓，对应循环卷积。滤波器在空间中等效为 (h(x, y) \ast f(x, y)) 的循环版本——从左边卷出去会从右边卷回来。解决办法：将 (f) 与 (h) 都零填充到 (P \ge M + k - 1)，强行变成线性卷积。

Q3. 高斯低通为什么无振铃？

答：高斯的 FT 仍是高斯。空域与频域都紧致连续、无跳变，不会产生 Gibbs 振荡。与之相反，理想低通在频域是矩形（跳变），其反变换为 sinc（有无穷振荡），卷积后边缘处震荡即振铃。

Q4. 如何从频谱判断图像内容？

答：

• 中心亮斑：DC 分量（图的平均亮度）
• 水平亮线：图像中有显著垂直边缘（梳子纹、栅栏）
• 径向对称衰减：自然图像（1/f 统计）
• 离散亮斑（非中心）：周期性纹理或干扰

Q5. FFT 为什么比 DFT 快？

答：利用 “(e^{-j2\pi/N}) 的周期性和对称性” 做分治。(N = 2^k) 下通过”偶 + 奇”递归拆分，复杂度从 (O(N^2)) 降到 (O(N \log N))。

Q6. 相位谱为什么比幅度谱更重要？

答：Oppenheim 实验表明，图像的结构信息（轮廓、位置） 主要编码在相位中，幅度只管能量。幅度与另一图置换后，重构图像看起来更像”相位那张图”。这也解释了为什么相位相关（Phase Correlation）可用作图像配准 —— 平移只改变相位。

Q7. 同态滤波为什么要取对数？

答：图像模型 (f = i \cdot r)（乘性）。对数把乘法变加法：(\ln f = \ln i + \ln r)。加性之后就可以线性滤波（分别处理低频照度和高频反射率）。IDFT 后再做 exp 回到原空间。

Q8. 一张含横条纹的图，如何去噪？

答：

1. FFT 观察频谱 → 会出现纵向的对称亮斑对（因为横条纹在空间横向是低频、纵向是高频）
2. 用陷波滤波器在那几对亮斑处 (H \to 0)
3. 其他频率 (H = 1)
4. IFFT 回空间域，条纹几乎消失