11  频域处理与 FFT

前几章里我们一直在空间域（spatial domain）工作：图像是像素的阵列，滤波是邻域内的加权运算，所有问题都围绕”哪里的灰度是多少”展开。本章换一个视角——频域（frequency domain）关心的不是”哪里”，而是”多快”：图像的灰度沿空间变化得有多快？缓慢起伏的背景是低频，锐利的边缘与细密的纹理是高频。这个视角并非数学游戏，它能解决空间域束手无策的实际问题。设想产线图像被一种周期性的电磁干扰污染——整幅图叠上了斜向的细密条纹。在空间域里，这种干扰无处不在、与图像内容纠缠在一起，均值滤波抹不掉它，中值滤波也无能为力；可是到了频域，这一整幅图的干扰坍缩成两个孤零零的亮点，用一个小圆把它们挖掉，干扰就近乎无损地消失了。本章将完整演示这个”杀手级”应用，而在此之前，先把频域的语言学会。

我们的实验场景（图 11.1）是一幅 512×512 的合成图像：一个大暗矩形和一个大圆提供低频的块状结构，中部是一条周期恰好为 8 px 的竖条纹带（已知频率的高频成分），下方是点阵字母 “FFT”（类文本细节）。每一种成分都会在频谱上留下可辨认的签名。

图 11.1: 512×512 实验场景：大矩形与大圆（低频块状结构）、周期 8 px 的竖条纹带（已知高频）、点阵文字 “FFT”（细节成分）。

11.1 二维傅里叶变换

傅里叶变换（Fourier transform）的核心思想是：任何信号都可以分解为一组不同频率的正弦波的叠加。对 \(W\times H\) 的离散图像 \(f[n,m]\)（\(n\) 为行、\(m\) 为列），二维离散傅里叶变换（2D Discrete Fourier Transform, DFT）定义为

\[ F[u,v] = \sum_{n=0}^{H-1}\sum_{m=0}^{W-1} f[n,m]\, e^{-j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)}, \]

其中 \(u,v\) 是水平与垂直方向的频率序号。逆变换（inverse DFT）把频率成分重新叠加回图像：

\[ f[n,m] = \frac{1}{WH}\sum_{u=0}^{W-1}\sum_{v=0}^{H-1} F[u,v]\, e^{+j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)}. \]

直接按定义计算 DFT 需要 \(O(N^2)\) 次运算（\(N=WH\)），512×512 的图像约要 \(7\times 10^{10}\) 次——不可接受。快速傅里叶变换（Fast Fourier Transform, FFT）利用旋转因子的对称性把复杂度降到 \(O(N\log N)\)，这正是频域方法得以实用的前提，也是本章实验把图像尺寸取为 2 的整数次幂的原因。

\(F[u,v]\) 是复数，工程上几乎总是把它拆成两个更直观的量来看。幅度谱（magnitude spectrum）\(|F[u,v]| = \sqrt{\mathrm{Re}^2 + \mathrm{Im}^2}\) 回答”频率为 \((u,v)\) 的正弦成分有多少”；相位谱（phase spectrum）\(\arg F[u,v]\) 回答”这个成分摆在哪里”——同样一组正弦波，相位不同，叠出来的图像就完全不同。滤波操作主要在幅度上做文章，但重建图像时相位必须原样保留，否则图像结构会面目全非。

有一个频率成分值得单独点名：\(u=v=0\) 时指数项恒为 1，故 \(F[0,0] = \sum f[n,m]\)，恰好是全图灰度之和，即 \(WH\) 倍的平均灰度。这就是直流分量（DC component）——章节 6 中”核权重之和等于直流增益”的说法正源于此：直流增益为 1 的滤波器不改变 \(F[0,0]\)，也就不改变图像的平均亮度。

直接把 \(|F|\) 画成图像几乎什么也看不见：直流分量的能量比绝大多数高频成分大几个数量级，线性灰度映射下整幅谱除中心一点外近乎全黑。标准做法是显示对数谱 \(\log(1+|F|)\)，把巨大的动态范围压缩到肉眼可辨的区间——本章所有谱图都按此绘制。

11.2 读谱

图 11.2 是实验场景的对数幅度谱（直流分量已移到图像中心）。读谱是一项可以速成的技能，对照场景逐一辨认。

图 11.2: 实验场景的对数幅度谱（DC 居中）。中心亮斑为低频能量，水平/垂直亮十字来自矩形与条纹带的横平竖直边缘，水平轴上 ±64 处的对称点对正是 8 px 竖条纹的签名。
  • 中心亮斑：低频能量。大矩形、大圆这类缓慢变化的块状结构，能量几乎全部集中在中心附近。
  • 水平与垂直的亮十字：矩形的水平/垂直边缘以及条纹带的水平边界都是”沿某一方向的突变”，突变在频域里沿垂直于边缘的方向铺开成一条线。
  • 水平轴上一对对称亮点：纯粹的周期条纹在频域里就是一对脉冲——这是本章最重要的对应关系。

这对亮点的位置可以先预测、再验证。周期为 \(T\) 像素的条纹，在宽 \(W\) 的图像里恰好重复 \(W/T\) 次，因此其谱峰出现在水平频率序号

\[ u = \frac{W}{T} = \frac{512}{8} = 64 \]

处，即中心两侧 \(\pm 64\) 的位置。实验中我们排除中心低频区后搜索全谱最大幅值，实测峰位为 \((\Delta u, \Delta v) = (-64, 0)\)——与预测精确命中。这种”已知周期 → 谱峰位置”的换算在排查产线周期性干扰时是直接可用的武器：量出干扰条纹的像素周期，就知道该去谱上哪里找它。

共轭对称性（conjugate symmetry）：实值图像的频谱满足 \(F[-u,-v] = F^*[u,v]\)，幅度谱因此关于中心对称——这就是谱峰总是成对出现的原因。后文做陷波时，每个干扰峰的共轭峰必须一并处理，否则逆变换的结果不再是实图像。

11.3 频域滤波

频域滤波的理论基石是卷积定理（convolution theorem）：空间域的卷积等价于频域的逐点相乘，

\[ f * h \;\longleftrightarrow\; F \cdot H. \]

这意味着 章节 6 的全部线性滤波都可以搬到频域执行：先对图像做 FFT，把谱与滤波器的频率响应 \(H\) 逐点相乘，再做逆变换。两条路径在数学上完全等价，工程上的取舍是计算量：\(K\times K\) 的核在空间域卷积每像素需 \(O(K^2)\) 次乘加，而频域路径的开销与核大小无关——核小用空间域，核大（数十像素以上）用频域。

最简单的频域滤波器是理想低通滤波器（ideal lowpass filter）：以直流为中心、半径 \(D_0\) 以内的 \(H=1\)，以外的 \(H=0\)，把高频”一刀切”地清零；理想高通（ideal highpass）则恰好相反。我们用截止半径 30 px（归一化频率 \(30/256 \approx 0.1172\)）做实验，结果见 图 11.3。

(a) 理想低通（截止半径 30）
(b) 理想高通（截止半径 30）
图 11.3: 理想低通与高通滤波结果。(a) 低通：8 px 条纹（频率 64 > 30）被整带抹平成均匀灰块，“FFT” 文字模糊，所有边缘附近出现一圈圈波纹状的振铃；(b) 高通：块状结构的内部被清空，只剩边缘轮廓，而条纹带因频率高于截止而完整保留。

低通结果验证了预期：条纹频率 64 远在截止半径 30 之外，整条条纹带被抹成均匀的灰色。但请注意矩形、圆和文字周围那一圈圈同心波纹——这是吉布斯振铃（Gibbs ringing），理想滤波器必须付出的代价。成因很明确：频域的矩形硬截断，其空间域对应物是 sinc 函数，而 sinc 有绵延不绝的振荡旁瓣；按卷积定理，频域相乘等于空间域同这个带拖尾的核做卷积，每条边缘都被拖出一串波纹。工程上因此很少使用理想滤波器，而是改用高斯或巴特沃斯（Butterworth）滤波器——频率响应平滑过渡，空间域的核不再振荡，代价只是截止边缘不那么”锋利”。

振铃的教训可以一句话记住：频域里切得越狠，空间域里拖得越长。这与 章节 6 中”高斯核频率响应无旁瓣、平滑效果干净”是同一枚硬币的两面。

高通结果同样符合谱的语言：半径 30 以内的低频被清零后，块状结构只剩下边缘轮廓（边缘是高频）；而条纹带几乎原样保留——它的频率 64 高于截止 30，在高通的通带之内。高通滤波因此常被用作边缘增强与背景抑制的手段。

11.4 周期噪声陷波

现在进入本章的核心演示。给场景叠加一个对角正弦干扰 \(40\sin\!\big(2\pi(r+c)/8\big)\)（沿对角方向周期约 5.7 px），模拟电磁干扰或机械振动造成的周期性条纹污染。含噪图像与干净图像的 RMSE 为 28.13——与理论值 \(40/\sqrt{2} \approx 28.3\)（正弦波的均方根幅值）吻合。空间域滤波对它无能为力：干扰的空间频率与场景里有用的 8 px 条纹相当接近，任何强到能压制干扰的平滑都会把有用条纹一起抹掉。

但在频域里，这个干扰只是两个点。图 11.4 展示完整流程。

(a) 含周期噪声
(b) 含噪图像的谱
(c) 陷波去噪结果
图 11.4: 周期噪声的陷波去除。(a) 全图被斜向细条纹覆盖，RMSE 28.13；(b) 干扰在谱上坍缩为对角线上 \((\pm 64, \pm 64)\) 处的一对共轭亮点；(c) 以半径 5 的小圆挖除两个亮点后逆变换，条纹干扰消失，RMSE 降至 0.61，有用的 8 px 竖条纹完好无损。

先做预测：干扰沿行、列方向的周期都是 8 px，谱峰应出现在 \((\pm 512/8, \pm 512/8) = (\pm 64, \pm 64)\)。实测搜索谱峰得到 \((-64, -64)\) 及其共轭峰 \((+64, +64)\)——再次精确命中。接下来是陷波（notch filtering）四步流程：

  1. 对含噪图像做 FFT（DC 居中）；
  2. 排除中心低频区，搜索幅值最大的干扰峰；
  3. 以峰位为圆心、半径 5 px 的小圆把复数谱置零，共轭峰同样处理；
  4. 逆变换回空间域，钳位到 [0, 255]。

结果（图 11.4 (c)）几乎与原图无异：RMSE 从 28.13 降到 0.61，干扰能量被消除了 98% 以上，而有用的竖条纹（峰位在水平轴 \(\pm 64\)，与干扰峰相距甚远）毫发无伤。这正是频域方法的不可替代之处：两种成分在空间域的”快慢”可以很接近，但只要方向或周期略有差异，它们在谱上就是分离的点——空间域滤波按”快慢”一刀切，频域陷波却可以指哪打哪。残余的 0.61 来自被挖掉的小圆内同时损失的少量图像本征能量——这提示陷波半径宁小勿大。

11.5 采样定理回看

有了频域语言，章节 1 中”每个周期至少采样两次”的混叠（aliasing）直觉终于可以升级为正式表述。奈奎斯特采样定理（Nyquist sampling theorem）：只要信号的最高频率成分低于采样频率的一半，就能从采样值完美重建原信号。频域里看得最清楚——采样使原信号的频谱以采样频率为间隔周期性复制，一旦信号带宽超过采样频率之半，相邻的频谱副本就会搬移重叠，高频成分伪装成低频混入信号，且事后无法分离。这也解释了 章节 10 中缩小图像前必须先做平滑的规矩：降采样等于降低采样频率，预滤波就是先把频谱裁剪到新的奈奎斯特限以内，宁可损失细节，也不能让细节变成虚假的低频伪影。

11.6 SciVision 实现

本章全部实验由 SCIMV::SciSvFFT 类完成，核心调用如下：

SCIMV::SciSvFFT fft;
SciImage fftImg;
// mode=0：直流分量移到谱中心；fftImg 为 2 通道 32F 图像，交错存储实部/虚部
fft.ApplyFFTGeneric(src, 0, &fftImg, NULL);

// 理想低通/高通：frequency 为归一化截止（相对半宽 W/2），0.1172 → 实测半径 30.0 px
float cutoff = 30.0f / (512 / 2);
SciImage lp;
fft.GenerateLowpass(cutoff, 0, 512, 512, &lp);

// 卷积定理的直接体现：谱与滤波器逐点相乘
SciImage fftLP, out32;
fft.ApplyFFTConvolution(fftImg, lp, &fftLP);
// 逆变换：务必请求 SCI_IMAGE_32F，再自行钳位到 [0,255]
fft.ApplyIFFTGeneric(fftLP, 0, &out32, true, SCI_IMAGE_32F);

ApplyFFTGeneric 输出的复数谱可以经 ImageData() 与 Step() 直接读写——第 \(r\) 行第 \(c\) 列的实部、虚部分别位于浮点偏移 r*stride + 2*c 与 r*stride + 2*c + 1。小节 11.4 的陷波正是这样手动把两个小圆内的复数置零实现的，SDK 不提供现成陷波接口也无妨。

使用这套 API 有三个实测确认的陷阱，如实记录：

  • IFFT 输出 8U 时会被 min-max 归一化：直接请求 8 位输出，灰度被整体拉伸，往返（FFT→IFFT）的 RMSE 高达 33.62，看似算法有错，实为归一化所致。请求 SCI_IMAGE_32F 拿到原始浮点结果、自行钳位后，往返 RMSE 为 0.0000——变换本身是无损的。
  • ApplyFFTGeneric 自带的幅度图近乎全黑：其 imageMagnitude 输出做的是线性归一化，直流独大，其余频率不可见。日志谱 \(\log(1+|F|)\) 需要自己从复数数据计算。
  • RemovePeriodicPatternsByFFT 不可用：该接口在头文件中明确标注”未实现”，周期噪声去除请按本章流程手动陷波。

生成本章全部图像的完整工程位于 code/fourier/，谱峰预测值与实测值的对照会直接打印在控制台。

工业案例：编织物上的网纹干扰

某纺织表面缺陷检测项目中，织物本身的经纬织纹是强烈的周期性结构，灰度起伏远大于断纱、油污等缺陷信号，缺陷在原图上几乎不可见。空间域滤波陷入两难：核小了织纹滤不干净，核大了缺陷连同织纹一起被糊掉。改用频域方案后问题迎刃而解——织纹的基频及其谐波在谱上是一组位置稳定的亮点对，对每个点对做小半径陷波，再逆变换回空间域，织纹被整体移除，缺陷在残差图上一目了然。现场调试的关键经验是陷波半径宁小勿大：缺陷能量同样分布在织纹峰附近，半径取大固然把织纹滤得更净，却会吞掉缺陷能量、降低检出对比度——以刚好覆盖谱峰为度。

11.7 小结

  • 频域回答”多快”，空间域回答”哪里”：幅度谱给出各频率成分的多少，相位谱给出它们的摆放位置；直流分量 \(F[0,0]\) 就是全图灰度和，对数显示 \(\log(1+|F|)\) 是看谱的标准姿势。
  • 谱是可以读、更可以预测的：周期 \(T\) 的条纹在宽 \(W\) 的图像中产生 \(u=W/T\) 处的共轭峰对——8 px 条纹预测 ±64，实测精确命中。
  • 卷积定理 \(f*h \leftrightarrow F\cdot H\) 打通两个域：小核用空间域卷积，大核走频域更快；频域硬截断会招致吉布斯振铃（sinc 拖尾），工程上用高斯/巴特沃斯软过渡替代理想滤波器。
  • 陷波是周期干扰的特效药：整幅图的周期污染在谱上只是几个点，小圆置零后逆变换即可——实验中 RMSE 从 28.13 降到 0.61，这是任何空间域滤波都做不到的。
  • 奈奎斯特定理的频域表述：信号最高频率必须低于采样频率之半，否则频谱副本搬移重叠产生混叠——这正是降采样前必须预滤波的根本原因。

使 FFT 实用化的快速算法源头是 Cooley 与 Tukey 1965 年的经典论文 (Cooley 和 Tukey 1965)，它把 DFT 的复杂度从 \(O(N^2)\) 降到 \(O(N\log N)\)；傅里叶变换在数字图像处理中的系统论述（含频域滤波与陷波）见 Gonzalez 与 Woods 的教材 (Gonzalez 和 Woods 2018)。关于频域方法在工业图像处理中的系统论述（包括最优滤波与频域特征），可进一步阅读 Steger 等人的著作 (Steger, Ulrich, 和 Wiedemann 2018)。