频率域与采样

频率域的视角

上一章讨论了空间域滤波——直接对像素的邻域进行加权求和。本章从一个全新的角度来理解图像和滤波:频率域(Frequency Domain)。

在空间域中,我们关注「每个像素的值是多少」;在频率域中,我们关注「图像包含哪些频率成分,各成分有多强」。傅里叶变换在这两种等价表示之间架起桥梁。频率域的视角能回答空间域分析难以解答的问题:为什么高斯滤波比均值滤波效果更好?为什么降低图像分辨率时会出现摩尔纹?为什么同一张混合图像在近处和远处看起来不同?

傅里叶的核心思想

1807 年,Fourier 提出了一个大胆的想法:任何函数都可以表示为不同频率的正弦和余弦函数的加权和。这一想法在当时备受争议——Lagrange、Laplace 等数学大家都对其持怀疑态度,直到 1878 年才被翻译为英文。但这一想法(在满足一定技术条件下)基本上是正确的,后来被称为傅里叶级数(Fourier Series)。

直觉上,正弦波 Asin(ωx+ϕ)A\sin(\omega x + \phi) 是基本构建模块,其中 AA 控制振幅,ω\omega 控制频率,ϕ\phi 控制相位。将足够多的正弦波叠加,就可以逼近任意信号。

方波的傅里叶级数

方波可以表示为奇数次谐波的叠加:

f(t)=k=1,3,5,1ksin(kt)f(t) = \sum_{\substack{k=1,3,5,\ldots}}^{\infty} \frac{1}{k} \sin(kt)

只用基频 k=1k=1 时是一条简单的正弦曲线;叠加三次谐波 k=3k=3 后波形开始出现平顶;继续叠加更高次谐波,波形逐步逼近方波。每增加一个分量,近似精度就提高一些。

重要信号

在傅里叶分析中,有几类信号扮演着基础角色。

脉冲信号(Impulse)δ[n]\delta[n]n=0n = 0 处取值 1,其余位置为 0。脉冲是卷积运算的单位元:g[n]δ[n]=g[n]g[n] * \delta[n] = g[n]。与平移脉冲 δ[nn0]\delta[n - n_0] 卷积则得到平移后的信号 g[nn0]g[n - n_0]

余弦和正弦波 cos(2πkn/N)\cos(2\pi k n / N)sin(2πkn/N)\sin(2\pi k n / N) 是傅里叶分析的基函数。参数 k[1,N/2]k \in [1, N/2] 表示在长度为 NN 的区间内完成的完整周期数——kk 越大,频率越高。

复指数(Complex Exponential)将余弦和正弦统一在一个框架中。由欧拉公式(Euler's Formula):

eiθ=cosθ+isinθ\e^{\mathrm{i}\theta} = \cos\theta + \mathrm{i}\sin\theta

复指数在复平面上对应单位圆上的一个旋转点——实部是余弦,虚部是正弦。离散时间的复指数信号为 ei2πkn/N\e^{\mathrm{i} 2\pi k n / N}。使用复指数可以将傅里叶变换写成紧凑的形式。推广到二维:

ei2π(k1n/N+k2m/M)\e^{\mathrm{i} 2\pi (k_1 n / N + k_2 m / M)}

其中 (k1,k2)(k_1, k_2) 分别控制水平和垂直方向的频率。

频谱

将信号分解为频率分量后,每个频率上的振幅构成该信号的频谱(Frequency Spectrum)。例如信号 g(t)=sin(2πft)+13sin(2π3ft)g(t) = \sin(2\pi f t) + \frac{1}{3}\sin(2\pi \cdot 3f \cdot t) 由基波(频率 ff,振幅 1)和三次谐波(频率 3f3f,振幅 1/31/3)两个分量组成,其频谱在 ff3f3f 处有两个峰。

频谱提供了一种互补的分析视角——不再关注信号在每个时刻的值,而是关注信号由哪些频率成分构成。

一维离散傅里叶变换

基变换

从数学角度看,傅里叶变换是一种基变换(Change of Basis)。将图像的所有像素排成向量 f\bm{f},变换可以写成矩阵乘法:

f^=Uf\hat{\bm{f}} = \bm{U} \bm{f}

其中 U\bm{U} 是变换矩阵,f^\hat{\bm{f}} 是频率域表示。除了傅里叶变换,小波变换(Wavelet Transform)和可控金字塔变换(Steerable Pyramid Transform)等也是常用的基变换,各有适用场景。

傅里叶变换的变换矩阵是酉矩阵(Unitary Matrix),满足 U1=U\bm{U}^{-1} = \bm{U}^*(逆等于共轭转置)。直觉上,酉矩阵保持向量的「长度」不变(Uf=f\|\bm{U}\bm{f}\| = \|\bm{f}\|),因此变换是完全可逆的,不丢失也不创造信息。

DFT 的定义

对长度为 NN 的离散信号 f[n]f[n],其离散傅里叶变换(Discrete Fourier Transform, DFT)为:

DFT 与逆 DFT

正变换(空间域 → 频率域):

F[u]=n=0N1f[n]ei2πun/N,u=0,1,,N1F[u] = \sum_{n=0}^{N-1} f[n] \, \e^{-\mathrm{i} 2\pi u n / N}, \quad u = 0, 1, \ldots, N-1

逆变换(频率域 → 空间域):

f[n]=1Nu=0N1F[u]e+i2πun/N,n=0,1,,N1f[n] = \frac{1}{N} \sum_{u=0}^{N-1} F[u] \, \e^{+\mathrm{i} 2\pi u n / N}, \quad n = 0, 1, \ldots, N-1

两者的区别仅在于指数的符号(- vs ++)和归一化因子 1/N1/N

uu 是频率索引(对应于前面的 kk,但范围扩展到 [0,N1][0, N-1]),F[u]F[u] 是复数,编码了频率分量 uu 的振幅和相位。直觉上,DFT 计算的是信号 ff 与频率为 uu 的复指数波的内积(逐点相乘再求和)——内积越大,该频率分量越显著。

一个简单的例子

N=4N = 4,信号 f=[1,2,3,4]f = [1, 2, 3, 4]DFTF[0]F[0] 是所有元素之和 1+2+3+4=101 + 2 + 3 + 4 = 10,除以 NN 就是信号的平均值 2.52.5——这就是直流分量。F[1]F[1] 衡量信号中「每 4 个采样点完成一个周期」的频率成分有多强。

DFT 可以写成矩阵形式 f^=Wf\hat{\bm{f}} = \bm{W}\bm{f},其中 W\bm{W} 的第 (u,n)(u, n) 个元素为 ei2πun/N\e^{-\mathrm{i} 2\pi un / N}。矩阵的每一行对应一个频率的基函数:第 0 行是常数(直流分量,即信号的平均值),后续各行的振荡频率依次递增。

幅度与相位

DFT 的输出 F[u]F[u] 是复数,由实部 R[u]R[u] 和虚部 I[u]I[u] 组成。从中可以提取两个物理量:

幅度(Magnitude)衡量频率 uu 的分量有多强:

F[u]=R[u]2+I[u]2|F[u]| = \sqrt{R[u]^2 + I[u]^2}

相位(Phase)编码该频率分量在空间上的偏移:

ϕ[u]=arctanI[u]R[u]\phi[u] = \arctan\frac{I[u]}{R[u]}

幅度和相位完整描述了 F[u]F[u]F[u]=F[u]eiϕ[u]F[u] = |F[u]| \e^{\mathrm{i}\phi[u]}

常见变换对

一些信号与其傅里叶变换之间存在优雅的对应关系。掌握这些变换对有助于快速判断滤波器的频率特性,也是理解后续采样理论的基础:

空间域 频率域 说明
脉冲 δ\delta 常数 1 脉冲包含所有频率且振幅相等
常数 1 脉冲 δ\delta 常数信号仅有零频率分量
方波(Box) sinc 函数 sinc 有无穷振荡的旁瓣
高斯函数 高斯函数 高斯在变换下保持形状
三角波(Tent) sinc2\text{sinc}^2 两个方波的卷积;卷积对应频率域乘法
梳状函数(Shah) 梳状函数 采样函数的变换仍是采样函数

高斯函数在傅里叶变换下仍为高斯——空间域宽高斯对应频率域窄高斯,反之亦然。这种「无旁瓣」的特性是高斯滤波器在图像处理中广受青睐的重要原因。

二维图像傅里叶变换

二维 DFT 的定义

将一维 DFT 推广到二维,对大小为 N×MN \times M 的图像 f[n,m]f[n, m]

F[u,v]=n=0N1m=0M1f[n,m]ei2π(un/N+vm/M)F[u, v] = \sum_{n=0}^{N-1} \sum_{m=0}^{M-1} f[n, m] \, \e^{-\mathrm{i} 2\pi (un/N + vm/M)}

逆变换为:

f[n,m]=1NMu=0N1v=0M1F[u,v]e+i2π(un/N+vm/M)f[n, m] = \frac{1}{NM} \sum_{u=0}^{N-1} \sum_{v=0}^{M-1} F[u, v] \, \e^{+\mathrm{i} 2\pi (un/N + vm/M)}

(u,v)(u, v) 是二维频率坐标,uu 对应垂直方向频率,vv 对应水平方向频率。

二维波

二维复指数 ei2π(un/N+vm/M)\e^{\mathrm{i} 2\pi (un/N + vm/M)} 在空间中表现为特定方向的平面波。(u,v)(u, v) 共同决定波的方向和频率:

  • (u,v)=(2,0)(u, v) = (2, 0):纯垂直方向的波,有 2 个完整周期
  • (u,v)=(3,1)(u, v) = (3, 1):斜方向的波
  • (u,v)=(7,5)(u, v) = (7, -5):更高频率的斜向波

(u,v)(u, v) 的模 u2+v2\sqrt{u^2 + v^2} 决定总频率的高低,而 arctan(v/u)\arctan(v/u) 决定波的传播方向。

幅度谱与相位谱

二维 DFT 的结果通常可视化为幅度谱相位谱两张图像。幅度谱的标准显示约定是将零频率(直流分量,即图像的平均值)置于中心,频率向外递增。由于直流分量通常远大于其他频率,一般取对数 log(1+F[u,v])\log(1 + |F[u,v]|) 来增强可视性。

1
2
3
4
5
6
7
8
import numpy as np
import matplotlib.pyplot as plt

img = plt.imread('image.png')[:, :, 0]      # 灰度化
F = np.fft.fft2(img)                         # 二维 DFT
F_shift = np.fft.fftshift(F)                 # 将零频率移到中心
magnitude = np.log(1 + np.abs(F_shift))      # 对数幅度谱
phase = np.angle(F_shift)                    # 相位谱

观察幅度谱可以得到以下直觉:

平移不改变幅度谱。移动图像只改变各频率分量的相位,不改变振幅——图像的频率组成不因空间位置变化而变化。

空间尺度与频率互逆。图像中尺度越小的细节(如精细纹理、锐利边缘),对应频率域中越高的频率。

方向性。含有强烈水平边缘的图像,其幅度谱在垂直轴上能量集中——水平边缘意味着亮度沿垂直方向快速变化。

一个有趣的实验是在频率域中选择性地将某些方向的频率分量置零再逆变换:如果将所有水平频率分量置零,逆变换后图像中的垂直边缘消失;反之亦然。这直观地展示了不同方向的频率分量各自编码了图像中对应方向的结构信息。

相位比幅度更重要

如果将图像 A 的幅度谱与图像 B 的相位谱组合并逆变换,得到的图像在视觉上几乎完全像图像 B。这说明相位编码了物体的空间结构和轮廓——人类视觉最依赖的信息。幅度更多影响整体纹理和对比度,对「像什么」的贡献较小。

傅里叶变换的性质

基本性质

线性F[αf+βg]=αF[f]+βF[g]\mathcal{F}[\alpha f + \beta g] = \alpha \mathcal{F}[f] + \beta \mathcal{F}[g]

实信号的对称性:实值信号的傅里叶变换关于原点共轭对称:F[u,v]=F[u,v]F[-u, -v] = F^*[u, v],因此幅度谱关于原点对称。

帕塞瓦尔定理(Parseval's Theorem):信号在空间域和频率域的能量相等:

nf[n]2=1NuF[u]2\sum_n |f[n]|^2 = \frac{1}{N} \sum_u |F[u]|^2

傅里叶变换不创造也不消灭能量,只是在不同频率分量间重新分配。

卷积定理

卷积定理(Convolution Theorem)是连接空间域与频率域的核心桥梁:

卷积定理

空间域的卷积等价于频率域的逐元素乘法

F[gh]=F[g]F[h]\mathcal{F}[g * h] = \mathcal{F}[g] \cdot \mathcal{F}[h]

反过来,频率域的乘法等价于空间域的卷积:

F1[GH]=F1[G]F1[H]\mathcal{F}^{-1}[G \cdot H] = \mathcal{F}^{-1}[G] * \mathcal{F}^{-1}[H]

卷积定理有两个重要的实际意义。

计算加速。直接卷积的复杂度为 O(N2M2)O(N^2 M^2)NN 为图像大小,MM 为核大小)。利用快速傅里叶变换(Fast Fourier Transform, FFT),先将图像和核变换到频率域(O(N2logN)O(N^2 \log N)),逐元素相乘,再逆变换回来。当核较大时(M>logNM > \log N),FFT 方法更快。

频率域滤波。在空间域应用滤波器等价于在频率域将图像频谱乘以滤波器的频率响应(Frequency Response)。这提供了设计滤波器的新思路——直接指定要保留或抑制的频率范围。

flowchart LR
    subgraph 空间域
        A["图像 f"] -- "卷积 *" --> C["输出"]
        B["滤波器 h"] --> C
    end
    subgraph 频率域
        D["F = FFT(f)"] -- "逐元素相乘" --> F["G = F · H"]
        E["H = FFT(h)"] --> F
        F -- "IFFT" --> G["输出"]
    end

    A -.-> D
    B -.-> E

    classDef spatial fill:#e3f2fd,stroke:#1565c0,stroke-width:2px
    classDef freq fill:#e8f5e8,stroke:#2e7d32,stroke-width:2px
    classDef result fill:#fff3e0,stroke:#ef6c00,stroke-width:2px

    class A,B spatial
    class D,E freq
    class C,G result
    class F freq

频率域视角下的滤波器

卷积定理为上一章观察到的一个现象提供了频率域的解释:为什么高斯滤波平滑自然,而均值滤波会产生方块状伪影

高斯滤波器的傅里叶变换仍然是高斯函数——一条平滑的钟形曲线,从低频到高频逐渐衰减,没有任何振荡。高频分量被平滑地抑制。

均值滤波器的傅里叶变换是 sinc 函数——主瓣之外有无穷个振荡的旁瓣。这些旁瓣意味着某些高频分量不仅未被抑制,反而可能被放大,导致振铃效应(Ringing Artifact)——在边缘附近出现波纹状伪影。

高斯函数的独特之处在于它是唯一在空间域和频率域都没有旁瓣的函数。此外,空间域的标准差 σ\sigma 与频率域的宽度成反比关系:σ\sigma 越大,频率域的高斯越窄,低通滤波效果越强——这与空间域中「σ\sigma 越大,平滑越强」的直觉一致。

采样

降采样与混叠

降采样(Downsampling)是减小图像分辨率的操作。最简单的 2 倍降采样就是丢弃每隔一行一列的像素,得到原图 1/41/4 大小的图像。但直接降采样可能导致严重的视觉伪影——混叠(Aliasing)。

从频率域理解采样过程需要分几步来看:

  1. 原始信号的频谱占据一定的带宽——假设最高频率为 fmaxf_{\max},频谱分布在 [fmax,fmax][-f_{\max}, f_{\max}] 范围内
  2. 采样操作等价于在空间域将信号乘以一个脉冲序列(梳状函数)。由卷积定理,空间域的乘法对应频率域的卷积。由于梳状函数的傅里叶变换仍是梳状函数,采样后的频谱变成原始频谱沿频率轴的周期性复制——以采样频率 fsf_s 为间隔,无限重复
  3. 采样率足够fs>2fmaxf_s > 2f_{\max})时,频谱副本之间有间隔,不会重叠。此时可以用一个理想低通滤波器「切出」中心的那一份频谱,完美重建原始信号
  4. 采样率不足fs<2fmaxf_s < 2f_{\max})时,频谱副本发生重叠——高频分量「伪装」成低频,混入原始频谱中,无法再被分离。这就是混叠,一种不可逆的信息损坏

混叠的常见表现

  • 电影中车轮似乎向后转(时间采样率不足以捕捉旋转频率)
  • 光线追踪中棋盘格图案崩解(空间采样率不足以表示高频纹理)
  • 彩色电视上条纹衬衫出现奇怪花纹(空间频率超过显示器的采样能力)

奈奎斯特-香农采样定理

奈奎斯特-香农采样定理(Nyquist-Shannon Sampling Theorem)给出了避免混叠的精确条件:

采样定理

以离散间隔对连续信号采样时,采样频率必须大于信号最高频率 fmaxf_{\max} 的两倍:

fs2fmaxf_s \ge 2 f_{\max}

满足此条件即可从采样信号完美重建原始信号。临界频率 2fmax2f_{\max} 称为奈奎斯特频率(Nyquist Frequency)。

直觉上,捕获频率为 ff 的正弦波至少需要每周期两个采样点(一峰一谷),否则该正弦波与一个更低频率的正弦波不可区分——这正是混叠的本质。想象对一个高频正弦波以过低的频率采样:采样点恰好可以被另一条频率完全不同的正弦曲线完美穿过,于是原始的高频信号在采样后「看起来」变成了低频信号。

抗混叠

混叠源于信号中存在超出采样能力的高频分量。解决方案有两条路径:

提高采样率。采集更多样本使采样频率满足奈奎斯特条件,但往往受硬件限制。

预滤波。在降采样前用低通滤波器去除超过新采样频率一半的频率分量。虽然丢失了部分高频信息,但有控制地丢失高频远好于让高频混叠到低频产生伪影。这就是抗混叠(Anti-aliasing)。

2 倍降采样的标准算法:

  1. 对原图像应用高斯低通滤波进行平滑
  2. 每隔一行一列取一个像素,得到 H/2×W/2H/2 \times W/2 的输出
1
2
3
4
from scipy.ndimage import gaussian_filter

blurred = gaussian_filter(img, sigma=1.0)    # 高斯预滤波(抗混叠)
downsampled = blurred[::2, ::2]              # 每隔一行一列取样

对比两组结果可以清楚看到:不做预滤波时,图像在低分辨率下出现严重的锯齿和摩尔纹;做了高斯预滤波后,图像虽然更模糊,但视觉上自然且无伪影。

混合图像

原理

混合图像(Hybrid Image)利用人类视觉系统在不同观看距离下对频率的敏感度差异,将两张图像融合为一张——近看是一张图像,远看是另一张。

制作方法是将图像 I1I_1 的低频分量与图像 I2I_2 的高频分量叠加:

Ihybrid=LowPass(I1)+HighPass(I2)I_{\text{hybrid}} = \text{LowPass}(I_1) + \text{HighPass}(I_2)

低通滤波用高斯模糊实现,高通滤波通过原图减去其低通版本得到:

HighPass(I)=ILowPass(I)\text{HighPass}(I) = I - \text{LowPass}(I)

在频率域中观察混合图像,可以清晰看到低频部分来自图像 I1I_1,高频部分来自图像 I2I_2

对比敏感度函数

混合图像有效的生理学基础是对比敏感度函数(Contrast Sensitivity Function, CSF)。经典实验揭示了人眼对空间频率的敏感度并不均匀:

  • 峰值灵敏度出现在约 6 个周期/度视角
  • 非常低的频率(对应视野中很大或很近的物体)不易被感知
  • 非常高的频率(对应很小或很远的细节)同样不易被感知

走近混合图像时,高频细节变得清晰可见(占据了灵敏度峰值附近的视角范围),于是我们看到高频图像;走远时,同样的物理频率对应更高的「每度视角周期数」,超出感知范围而变得不可见,此时低频分量主导感知。

Campbell-Robson 图表是一个直观体验 CSF 的工具:图表中水平方向频率递增,垂直方向对比度递减。在正常观看距离下,你会注意到中间频率处能看到的对比度最低的条纹(灵敏度最高),而两侧的条纹在更高的对比度下才可见——这正是 CSF 曲线的直接可视化。

人类视觉的多尺度处理

人类视觉系统的早期处理阶段对输入进行多尺度的边缘和斑块检测——相当于一组不同方向和尺度的带通滤波器。研究表明,感知线索主要由中高频信息主导:我们辨认物体依靠的是边缘和轮廓(中高频),而非大面积的亮度分布(低频)。Dalí 的画作和 Gilbert 的「All Is Vanity」等艺术作品正是利用了这种多尺度感知特性——从不同距离观看时,你会「看到」完全不同的图像。

这一生物学发现启发了计算机视觉中的多尺度图像表示方法(如图像金字塔),将在后续章节中讨论。