从 O(n²) 暴力到 O(n·log n) 蝶形运算——OpenCV dft() 源码全链路拆解，附频域去条纹噪声实战

一、你的图像处理为什么需要频域？

做图像处理的程序员，十个里面有九个是从空间域起步的。模糊用 GaussianBlur()，锐化用拉普拉斯算子，去噪用中值滤波——这些操作直觉上很好理解，每个像素和它的邻居做一轮加权平均就完事了。

但总有一天你会碰到这样的场景：一张图片上出现了规律性的条纹干扰，间距固定、方向一致，像是被某种周期信号污染了。你试了高斯模糊，条纹确实淡了一点，但图像整体也糊了；你试了中值滤波，条纹纹丝不动——因为中值滤波擅长对付椒盐噪声，对这种有规律的周期性干扰束手无策。

这时候你需要换一个视角来看这张图像。空间域里看到的是像素值的起伏，而频域里看到的是"这张图包含了哪些频率的信号"。那些规律性的条纹，在频谱图上会变成几个刺眼的亮点——精准地挖掉这些亮点，再变换回去，条纹就消失了，而图像的其他细节几乎不受影响。

这就是频域滤波的威力：它能精准定位并消除特定频率的干扰，而空间域滤波做不到这么精细的频率选择性。

要进入频域，核心工具就是离散傅里叶变换（DFT）。OpenCV 用一个函数 cv::dft() 封装了整个变换过程，但这个函数背后是 4722 行的 dxt.cpp，包含了混合基 FFT、位反转排列、旋转因子预计算、SSE3 SIMD 加速、IPP 硬件加速、OpenCL GPU 加速等一整套工程优化。今天我们就从 DFT 的数学定义出发，一路拆到 OpenCV 源码的每一个关键环节，最后用频域陷波滤波器消除一张图像上的周期性条纹噪声。

二、离散傅里叶变换的物理意义——空间频率与图像纹理

2.1 什么是"空间频率"

在信号处理里，"频率"指的是信号每秒振荡多少次。但图像不是随时间变化的信号，图像是随空间位置变化的二维信号。所以图像处理里用的是"空间频率"——单位长度内灰度值变化的次数。

打个比方：一张棋盘格图像，黑白交替非常密集，那它的空间频率就高；一片均匀灰色的区域，灰度值几乎不变，空间频率就低。所以：

• 低频分量对应图像中大面积的、缓慢变化的区域——比如天空、背景色块

• 高频分量对应像素值剧烈变化的区域——边缘、纹理、噪声

• 特定频率的周期信号会在频谱上形成离散的亮点——这正是我们后面要定位和消除的目标

2.2 DFT 的数学定义

一维 DFT 的定义很简洁。给定一个长度为 N 的序列 x[n]，它的 DFT 为：

X[k] = Σ(n=0 to N-1) x[n] · exp(-j·2π·k·n/N), k = 0, 1, …, N-1

其中 exp(-j·2π·k·n/N) 是旋转因子（twiddle factor），本质上就是单位圆上等间隔采样的复数点。OpenCV 源码中这个公式对应的是 core.hpp 第 2146 行的注释：F(N)_jk = exp(-2πi·j·k/N)。

对于二维图像，DFT 可以分解为先对每一行做一维 DFT，再对每一列做一维 DFT。数学表达为：

Y = F(M) · X · F(N)

这也是 OpenCV 2D DFT 的实际实现策略——OcvDftImpl 类中 stages[0]=0（行变换）、stages[1]=1（列变换）的两阶段流程。

2.3 暴力 DFT 的计算量

直接按定义计算 DFT，对于每一个输出 X[k]，需要遍历所有 N 个输入 x[n]，做一次复数乘法和一次复数加法。N 个输出就是 N×N = O(N²) 次复数运算。

对于一个 512×512 的灰度图像来说，二维 DFT 相当于先做 512 次长度 512 的一维 DFT（行方向），再做 512 次长度 512 的一维 DFT（列方向），总共 2×512×512² ≈ 2.68 亿次复数运算。这个计算量在实时处理场景下是不可接受的。

快速傅里叶变换（FFT）把复杂度降到了 O(N·log N)。从量级上比较，暴力 DFT 的 N² = 262144，FFT 的 N·log₂N = 512×9 = 4608，量级上相差约 57 倍（实际蝶形运算的复数乘法次数约为 (N/2)·log₂N = 2304 次，但每次蝶形同时产生两个输出，总体工作量远小于暴力方法）。这就是为什么 OpenCV 的 dft() 不是"按定义算"，而是实现了一套精心优化的 FFT 算法。

三、从暴力 DFT 到 FFT：Cooley-Tukey 分治思想

3.1 核心思想：把大问题拆成小问题

FFT 的核心是分治法。最经典的 Cooley-Tukey 算法是这样的：把 N 点 DFT 分解成两个 N/2 点 DFT（一个处理偶数下标，一个处理奇数下标），然后用 O(N) 的蝶形运算合并结果，递归下去直到子问题大小为 1。

N 点 → 2 个 N/2 点 → 4 个 N/4 点 → … → N 个 1 点

每层合并的代价是 O(N)，总共 log₂N 层，所以总复杂度 O(N·log N)。

3.2 蝶形运算

所谓蝶形运算，就是 FFT 合并步骤中的基本操作单元。以 radix-2 为例，它的核心公式是：

a' = a + W·b
b' = a – W·b

其中 W 是旋转因子 exp(-j·2π·k/N)。一次蝶形运算只需要一次复数乘法和两次复数加法，就能同时算出两个输出——这比直接计算高效得多。

3.3 OpenCV 的混合基策略

经典的 Cooley-Tukey 要求 N 是 2 的幂。但实际应用中图像尺寸五花八门，不可能总是 2 的幂。OpenCV 采用了混合基（mixed-radix）FFT：把 N 分解为 2、3、5 的乘积，然后依次执行 radix-2、radix-3、radix-4、radix-5 的蝶形运算。

DFTFactorize() 函数（dxt.cpp 第 164 行）负责因子分解，我们来看这段源码：

// dxt.cpp 第164-206行
static int
DFTFactorize( int n, int* factors )
{
int nf = 0, f, i, j;

if( n <= 5 )
{
factors[0] = n;
return 1;
}

// 先提取最大的2的幂因子
f = (((n – 1)^n)+1) >> 1;
if( f > 1 )
{
factors[nf++] = f;
n = f == n ? 1 : n/f;
}

// 再提取3、5等奇数因子
for( f = 3; n > 1; )
{
int d = n/f;
if( d*f == n )
{
factors[nf++] = f;
n = d;
}
else
{
f += 2;
if( f*f > n )
break;
}
}

if( n > 1 )
factors[nf++] = n;

// 将最大的因子放前面（2的幂因子在最前）
f = (factors[0] & 1) == 0;
for( i = f; i < (nf+f)/2; i++ )
CV_SWAP( factors[i], factors[nf-i-1+f], j );

return nf;
}

这段代码有几个精妙之处值得注意：

第一，2 的幂因子的提取方式。(((n-1)^n)+1) >> 1 这个位运算技巧，能直接算出 n 中包含的最大 2 的幂因子。比如 n=360=8×45，二进制 360=101101000，359=101100111，XOR 得 000001111=15，加 1 得 16，右移 1 位得 8——正好是 360 的最大 2 的幂因子。一行代码搞定，不需要循环除 2。

第二，因子排列顺序。最终排好的 factors 数组，2 的幂因子放在最前面，其他奇数因子按照从大到小排列（源码中对奇数因子做了反转操作）。这个顺序不是随意的——它决定了后续蝶形运算的执行顺序：先做 radix-4/radix-2（处理 2 的幂因子），再做 radix-3、radix-5（处理奇数因子）。2 的幂因子用 radix-4 处理效率最高，所以放在前面优先执行。

第三，对大素数的处理。如果 n 包含大于 5 的素数因子（比如 n=7 或 n=11），那个素因子会被当做一个"通用基"处理，效率显著下降。这也是为什么 getOptimalDFTSize() 要把尺寸填充到只包含 2、3、5 因子的数——这一点我们后面会详细展开。

四、OpenCV dft() 源码全链路拆解

掌握了 FFT 的基本思想之后，我们来看 OpenCV 是怎么把这套理论变成高性能工程代码的。整个调用链从 cv::dft() 入口开始，经过多层加速分派，最终落到混合基蝶形运算的核心实现。

4.1 入口函数：四层加速分派

// dxt.cpp 第3505-3552行
void cv::dft( InputArray _src0, OutputArray _dst, int flags, int nonzero_rows )
{
CV_INSTRUMENT_REGION();

// 第一层：AMD clFFT GPU加速
#ifdef HAVE_CLAMDFFT
CV_OCL_RUN(ocl::haveAmdFft() && …, ocl_dft_amdfft(_src0, _dst, flags))
#endif

// 第二层：OpenCL GPU加速
#ifdef HAVE_OPENCL
CV_OCL_RUN(_dst.isUMat() && _src0.dims() <= 2,
ocl_dft(_src0, _dst, flags, nonzero_rows))
#endif

Mat src0 = _src0.getMat(), src = src0;
// … 类型检查、输出矩阵创建 …

// 第三层：HAL 硬件抽象层（可能是 IPP）
Ptr<hal::DFT2D> c = hal::DFT2D::create(src.cols, src.rows, depth,
src.channels(), dst.channels(), f, nonzero_rows);
// 第四层：软件回退（混合基FFT）
c->apply(src.data, src.step, dst.data, dst.step);
}

OpenCV 的 dft() 入口只有 47 行，但它体现了 OpenCV 一贯的工程哲学：多层加速分派，从最快的硬件加速开始尝试，逐层回退到软件实现。

加速优先级从高到低是：AMD clFFT（GPU）→ OpenCL（GPU）→ IPP（CPU SIMD）→ 软件回退（混合基 FFT + SSE3）。这意味着同样调用 cv::dft()，在不同的硬件环境下，实际执行的代码路径可能完全不同。

CV_OCL_RUN 宏是 OpenCV 的通用 OpenCL 分派机制——它会先检查目标 Mat 是否是 UMat（统一内存矩阵），如果是并且 OpenCL 可用，就走 GPU 路径。如果 GPU 路径执行失败或不可用，就自动回退到 CPU 路径，对调用者完全透明。

4.2 HAL 层与 OcvDftImpl：2D DFT 的行列分解

当 GPU 加速不可用时，执行流程进入 hal::DFT2D::create()，最终创建 OcvDftImpl 对象。这个类的 init() 方法实现了 2D DFT 的行列分解策略：

// dxt.cpp 第2888-2910行
DftDims dims = determineDims(height, width, isRowTransform, isContinuous);
if (dims == TwoDims)
{
stages.resize(2);
if (mode == InvCCSToReal || mode == InvComplexToReal)
{
stages[0] = 1; // 逆变换：先列后行
stages[1] = 0;
}
else
{
stages[0] = 0; // 正变换：先行后列
stages[1] = 1;
}
}

正变换先行后列，逆变换先列后行，这个顺序的选择跟 CCS 紧凑存储格式有关，保证了对称性填充的正确性。每个阶段创建一个 hal::DFT1D 对象来执行一维 DFT，最终 apply() 方法按照 stages 顺序依次调用它们。

4.3 旋转因子预计算：DFTInit()

蝶形运算中频繁使用旋转因子 W = exp(-j·2π·k/N)。如果每次蝶形运算都现算 cos/sin，性能会很差。OpenCV 在初始化阶段就把所有旋转因子预计算好，存到 wave[] 数组里。

// dxt.cpp 第343-353行
if( (n0 & (n0-1)) == 0 )
{
// N 是 2 的幂：直接查表
w.re = w1.re = DFTTab[m][0];
w.im = w1.im = -DFTTab[m][1];
}
else
{
// N 不是 2 的幂：用 cos/sin 计算
t = -CV_PI*2/n0;
w.im = w1.im = sin(t);
w.re = w1.re = std::sqrt(1. – w1.im*w1.im);
}

这里有一个细节值得关注：当 N 是 2 的幂时，OpenCV 不用调 cos()/sin()，而是直接查预计算好的 DFTTab[] 表（dxt.cpp 第 95-129 行），里面存了 32 个精度达 17 位有效数字的双精度值。用一次查表代替一次三角函数计算，在大量调用场景下节省了可观的时间。

然后通过递推关系 w_new = w * w1（复数乘法展开就是两次实数乘法两次实数加法）生成所有旋转因子，避免了对每个 k 都调用 cos(2πk/N) 和 sin(2πk/N)。

4.4 位反转排列：bitrevTab 查找表

FFT 的蝶形运算要求输入数据按照"位反转"顺序排列。比如 N=8 时，下标 0,1,2,3,4,5,6,7 的二进制表示是 000,001,010,011,100,101,110,111，位反转后变成 000,100,010,110,001,101,011,111，即下标 0,4,2,6,1,5,3,7。

暴力计算位反转需要对每个下标做 log₂N 次位操作。OpenCV 用一个 256 字节的查找表 bitrevTab[256] 一步到位：

// dxt.cpp 第75-93行
static unsigned char bitrevTab[] =
{
0x00,0x80,0x40,0xc0,0x20,0xa0,0x60,0xe0,
0x10,0x90,0x50,0xd0,0x30,0xb0,0x70,0xf0,
// … 共256个值
};

// dxt.cpp 第158-162行
#define BitRev(i,shift) \\
((int)((((unsigned)bitrevTab[(i)&255] << 24)+ \\
((unsigned)bitrevTab[((i)>>8)&255] << 16)+ \\
((unsigned)bitrevTab[((i)>>16)&255] << 8)+ \\
((unsigned)bitrevTab[((i)>>24)])) >> (shift)))

bitrevTab[i] 存的就是 8 位整数 i 的位反转结果。对于超过 8 位的整数，BitRev 宏把 32 位整数拆成 4 个字节，分别查表后拼接，再右移对齐——4 次查表代替了 32 次位操作。对于 N≤256 的情况（这在图像处理中非常常见，比如 256 点 FFT），甚至只需要一次查表。

4.5 混合基蝶形运算：radix-2/3/4/5

初始化完旋转因子和位反转表后，进入核心的蝶形运算阶段。DFT<T>() 模板函数（第 840 行）分三步执行：

第 0 步：数据洗牌（Shuffle）。根据位反转表重排输入数据，如果 src 和 dst 不同就直接拷贝到正确位置，如果是 in-place 变换就做原地交换。

第 1 步：处理 2 的幂因子。先尽可能多地执行 radix-4 蝶形运算（一次处理 4 个元素，效率最高），剩余的做 radix-2：

// dxt.cpp 第980-1064行
n = 1;
// 先做 radix-4
if( (c.factors[0] & 1) == 0 )
{
// 如果有 SSE3，用 SIMD 加速的 radix-4
if( c.factors[0] >= 4 && c.haveSSE3)
{
DFT_VecR4<T> vr4;
n = vr4(dst, c.factors[0], c.n, dw0, wave);
}
// 标量 radix-4 回退
for( ; n*4 <= c.factors[0]; )
{
nx = n;
n *= 4;
dw0 /= 4;
// … radix-4 蝶形运算核心循环 …
}
// 剩余的 radix-2
for( ; n < c.factors[0]; )
{
n *= 2;
dw0 /= 2;
// … radix-2 蝶形运算 …
}
}

为什么优先用 radix-4 而不是 radix-2？因为一个 radix-4 蝶形级可以替代两级 radix-2 蝶形。两级 radix-2 需要 2×(N/2) = N 次复数乘法，而一个 radix-4 级只需要 3N/4 次——其中基本 4 点 DFT 内部的旋转因子 W⁰=1 不需要乘法，加上 W²=-1 可以简化为取反，这些特殊值使得复数乘法次数比 radix-2 少约 25%。

第 2 步：处理奇数因子。按照 factors 数组中的顺序，依次执行 radix-3、radix-5 和通用基的蝶形运算：

// dxt.cpp 第1067-1174行
for( f_idx = (c.factors[0]&1) ? 0 : 1; f_idx < c.nf; f_idx++ )
{
int factor = c.factors[f_idx];
if( factor == 3 )
{
DFT_R3<T> or DFT_VecR3<T>(SSE3); // radix-3 蝶形
}
else if( factor == 5 )
{
DFT_R5<T>(); // radix-5 蝶形
}
else
{
// 通用基：效率较低，O(factor²) 操作
}
}

Radix-3 蝶形运算（DFT_R3<T>，第 438 行）利用了 sin(120°) = √3/2 ≈ 0.866 这个常数来避免三角函数计算。Radix-5 蝶形运算（DFT_R5<T>，第 479 行）更复杂，用了 5 个预计算常数（fft5_2 到 fft5_5）。如果遇到 7、11 等大素数因子，就走通用基路径——这个路径的内层循环是 O(factor²) 的，所以遇到大素数因子性能会急剧下降。

4.6 SSE3 SIMD 加速：复数乘法的指令优化

对于 float 类型的 DFT，OpenCV 提供了 SSE3 特化版本。核心优化在于复数乘法的 SIMD 实现：

// dxt.cpp 第562-574行
inline __m128 complexMul(const Complex<float>* const a,
const Complex<float>* const b) {
const __m128 z = _mm_setzero_ps();
const __m128 neg_elem0 = _mm_set_ps(0.0f,0.0f,0.0f,-0.0f);
const __m128 v_a = _mm_loadl_pi(z, (const __m64*)a);
const __m128 v_b = _mm_loadl_pi(z, (const __m64*)b);
// 交叉排列实部虚部
const __m128 v_a_riri = _mm_shuffle_ps(v_a, v_a, _MM_SHUFFLE(0,1,0,1));
const __m128 v_b_irri = _mm_shuffle_ps(v_b, v_b, _MM_SHUFFLE(1,0,0,1));
const __m128 mul = _mm_mul_ps(v_a_riri, v_b_irri);
const __m128 xored = _mm_xor_ps(mul, neg_elem0);
return _mm_hadd_ps(xored, z);
}

标量复数乘法 (a_re + j·a_im)·(b_re + j·b_im) 需要 4 次乘法和 2 次加减法。这段 SSE3 代码巧妙地利用了 _mm_shuffle_ps 交叉排列实部虚部，_mm_mul_ps 一次做 4 路并行乘法，_mm_xor_ps 加符号位翻转（替代减法），最后 _mm_hadd_ps 水平相加得到结果——整个复数乘法只用了 6 条 SIMD 指令，而且利用了 SSE3 独有的 _mm_hadd_ps（水平加法）指令。

Radix-4 的 SIMD 特化（DFT_VecR4<float>，第 651 行）更进一步：它把 4 个蝶形运算打包成一组 SIMD 操作，用 _mm_addsub_ps（SSE3 的交替加减指令）同时完成正负号处理，吞吐量比标量版本高 3-4 倍。

4.7 RealDFT：实数信号的 N/2 优化

图像数据是实数（灰度值 0-255），而通用 DFT 处理的是复数输入。如果把实数序列当复数处理（虚部全为 0），就浪费了一半的计算量。OpenCV 用 RealDFT<T>() 函数（第 1211 行）实现了一个巧妙的优化：

核心思路：把长度为 N 的实数序列 x[n] 重新解释为长度 N/2 的复数序列——x[0]+j·x[1], x[2]+j·x[3], …——先做一个 N/2 点的复数 DFT，然后用一趟 O(N) 的后处理还原出完整的 N 点实数 DFT 结果。

// dxt.cpp 第1274-1333行 RealDFT 核心逻辑（N为偶数时）
T scale2 = scale*(T)0.5;
int n2 = n >> 1;

// 修改 factors：将第一个因子除以2
c.factors[0] >>= 1;

// 把实数序列当成 N/2 点复数序列做 DFT
OcvDftOptions sub_c = c;
sub_c.n = n2;
DFT(sub_c, (Complex<T>*)src, (Complex<T>*)dst);

// 恢复 factors
c.factors[0] <<= 1;

// 后处理：从 N/2 点复数 DFT 结果还原 N 点实数 DFT
t = dst[0] – dst[1];
dst[0] = (dst[0] + dst[1])*scale;
dst[1] = t*scale;

const Complex<T> *wave = (const Complex<T>*)c.wave;

for( j = 2, wave++; j < n2; j += 2, wave++ )
{
// 分离偶数下标和奇数下标的频谱分量
h2_re = scale2*(dst[j+1] + t);
h2_im = scale2*(dst[n-j] – dst[j]);
h1_re = scale2*(dst[j] + dst[n-j]);
h1_im = scale2*(dst[j+1] – t);

// 用旋转因子校正奇数分量的相位
t = h2_re*wave->re – h2_im*wave->im;
h2_im = h2_re*wave->im + h2_im*wave->re;
h2_re = t;
t = dst[n-j-1];

// 组合成最终结果
dst[j-1] = h1_re + h2_re;
dst[n-j-1] = h1_re – h2_re;
dst[j] = h1_im + h2_im;
dst[n-j] = h2_im – h1_im;
}

这个优化让实数 DFT 的计算量几乎减半——N/2 点复数 DFT 的工作量约为 N 点的一半，加上一趟 O(N) 的后处理，总体比直接做 N 点复数 DFT 快近 2 倍。这也是为什么 OpenCV 要求输入是实数时不需要手动转成复数双通道矩阵的原因——它内部会自动走 RealDFT 这条更快的路径。

输出结果使用 CCS（Complex-Conjugate-Symmetrical）紧凑格式存储。实数信号的 DFT 结果具有共轭对称性：X[k] = X*[N-k]，所以 N 个复数输出中只有 N/2+1 个是独立的。CCS 格式把这些独立分量紧凑排列在一个单通道实数矩阵中：

Re(X[0]), Re(X[1]), Im(X[1]), Re(X[2]), Im(X[2]), …, Re(X[N/2])

这种格式节省了一半的存储空间，但也意味着后续用 mulSpectrums() 做频域乘法时，需要正确处理 CCS 格式的特殊布局——这一点我们马上就会看到。

4.8 getOptimalDFTSize()：预计算的 2^p·3^q·5^r 查找表

前面提到，混合基 FFT 只对因子为 2、3、5 的 N 高效。图像的实际尺寸很可能包含 7、11、13 等大素数因子，导致 FFT 性能急剧下降。getOptimalDFTSize() 的作用就是找到一个不小于给定尺寸的、只包含 2、3、5 因子的最小整数。

OpenCV 的实现出人意料地简单——一个预计算的查找表加上二分查找：

// dxt.cpp 第4461行开始
static const int optimalDFTSizeTab[] = {
1, 2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32,
36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, …
// 总共约 460+ 个值，覆盖到 2,125,764,000
};

// dxt.cpp 第4644行
int cv::getOptimalDFTSize( int size0 )
{
int a = 0, b = sizeof(optimalDFTSizeTab)/sizeof(optimalDFTSizeTab[0]) – 1;
if( (unsigned)size0 >= (unsigned)optimalDFTSizeTab[b] )
return -1;

while( a < b )
{
int c = (a + b) >> 1;
if( size0 <= optimalDFTSizeTab[c] )
b = c;
else
a = c+1;
}

return optimalDFTSizeTab[b];
}

这个查找表 optimalDFTSizeTab[] 包含了从 1 到约 21 亿范围内所有形如 2^p·3^q·5^r 的数，按升序排列，总共约 460 个条目。二分查找的时间复杂度 O(log 460) ≈ 9 次比较，几乎是零开销。

实际使用时，典型的模式是：

// 计算卷积时的标准做法
Size dftSize;
dftSize.width = getOptimalDFTSize(A.cols + B.cols – 1);
dftSize.height = getOptimalDFTSize(A.rows + B.rows – 1);
Mat tempA(dftSize, A.type(), Scalar::all(0));

比如一个 640×480 的图像，getOptimalDFTSize(640) 返回 640（640=2^7·5，正好是 2 和 5 的乘积），而 getOptimalDFTSize(480) 返回 480（480=2^5·3·5）。但如果图像是 641×481，那就会被填充到 648×486（648=2^3·3^4，486=2·3^5），多出来的像素用零填充。

零填充虽然增加了数据量，但 FFT 对 648 点的处理速度可能比 641 点快数倍，因为 641 是一个素数——它在因子分解时会作为一个"通用基"被处理，内层循环是 O(641²) 的。而 648 = 8×81 = 8×3⁴，所有因子都是 2 和 3，全程使用高效的 radix-4 和 radix-3 蝴蝶运算。

五、卷积定理：空间域卷积 = 频域乘法

5.1 卷积定理的数学本质

卷积定理是频域滤波的理论基石，它说的是：

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

用公式表达就是：

如果 h * x = y（空间域卷积）
那么 H · X = Y（频域逐元素乘法）

其中 H、X、Y 分别是 h、x、y 的 DFT。反过来也成立——频域的乘法等价于空间域的卷积。

这个定理的工程意义非常直接：空间域卷积的复杂度是 O(M·N·K²)（K 是卷积核大小），而频域乘法的复杂度是 O(M·N·log(M·N))（包含正反 DFT 和频域乘法）。当卷积核比较大时（比如 K>30），频域方法比空间域快得多。需要注意的是，OpenCV 的 filter2D() 并不会自动切换到频域实现——它始终使用空间域卷积（配合 IPP/SIMD 加速）。如果需要频域卷积的性能优势，必须像下面的示例那样手动实现。

5.2 频域卷积的标准流程

OpenCV 的 core.hpp 中给出了一个经典的频域卷积示例（第 2194 行），它完美展示了卷积定理的工程落地：

// core.hpp 第2194-2232行 DFT-based convolution示例
void convolveDFT(InputArray A, InputArray B, OutputArray C)
{
C.create(abs(A.rows – B.rows)+1, abs(A.cols – B.cols)+1, A.type());
Size dftSize;
// 1. 计算最优DFT尺寸（包含零填充）
dftSize.width = getOptimalDFTSize(A.cols + B.cols – 1);
dftSize.height = getOptimalDFTSize(A.rows + B.rows – 1);

// 2. 零填充
Mat tempA(dftSize, A.type(), Scalar::all(0));
Mat tempB(dftSize, B.type(), Scalar::all(0));
Mat roiA(tempA, Rect(0,0,A.cols,A.rows));
A.copyTo(roiA);
Mat roiB(tempB, Rect(0,0,B.cols,B.rows));
B.copyTo(roiB);

// 3. 正向DFT（利用nonzeroRows加速）
dft(tempA, tempA, 0, A.rows);
dft(tempB, tempB, 0, B.rows);

// 4. 频域乘法（卷积定理的核心）
mulSpectrums(tempA, tempB, tempA);

// 5. 逆向DFT
dft(tempA, tempA, DFT_INVERSE + DFT_SCALE, C.rows);

// 6. 裁剪结果
tempA(Rect(0, 0, C.cols, C.rows)).copyTo(C);
}

这段代码中 nonzeroRows 参数值得注意：dft(tempA, tempA, 0, A.rows) 中的 A.rows 告诉 DFT 函数"只有前 A.rows 行有非零数据"，这样 DFT 可以跳过后面全零行的计算，节省约 30-50% 的时间。这也是 OpenCV 为什么在 OcvDftImpl 中专门维护 nonzero_rows 字段的原因。

六、mulSpectrums() 源码：频域乘法的工程实现

mulSpectrums() 是 cv::dft() 的"搭档函数"，它实现了频域逐元素乘法（或共轭乘法）。表面上看这个函数应该很简单——就是逐元素做复数乘法嘛——但实际上它要处理两种不同的频谱格式，代码量并不小。

// dxt.cpp 第3719-3774行
void cv::mulSpectrums( InputArray _srcA, InputArray _srcB,
OutputArray _dst, int flags, bool conjB )
{
// OpenCL 加速分派
CV_OCL_RUN(_dst.isUMat() && …, ocl_mulSpectrums(…))

Mat srcA = _srcA.getMat(), srcB = _srcB.getMat();
int depth = srcA.depth(), cn = srcA.channels(), type = srcA.type();

CV_Assert( type == CV_32FC1 || type == CV_32FC2 ||
type == CV_64FC1 || type == CV_64FC2 );

_dst.create( srcA.rows, srcA.cols, type );
Mat dst = _dst.getMat();

// 支持 in-place：A 可以同时作为输出，B 不行（需clone）
if (dst.data == srcB.data)
srcB = srcB.clone();

bool is_1d = (flags & DFT_ROWS) || (rows == 1) || …;
bool isCN1 = cn == 1; // 是否是CCS格式（单通道实数）

// 根据类型分派到模板实现
if (depth == CV_32F)
{
if (!conjB)
mulSpectrums_Impl<float, false>(…);
else
mulSpectrums_Impl<float, true>(…);
}
// … double 类型类似 …
}

这里有两个关键的工程细节：

第一，CCS 格式的特殊处理。当输入是单通道实数矩阵（即 CCS 紧凑格式）时，第一列和最后一列（如果 N 是偶数）存的是纯实数，需要用 mulSpectrums_processCols() 单独处理。中间列存的是 Re、Im 交替排列的复数对，用 mulSpectrums_processRows() 处理。

第二，conjB 模板参数控制卷积 vs 相关。conjB=false 时做标准乘法 A·B，对应卷积操作；conjB=true 时做共轭乘法 A·B*，对应互相关操作。OpenCV 用模板参数而不是运行时 if 分支来区分这两种模式，编译器会为两种情况各生成一份特化代码，避免了循环内的分支预测开销。

频域乘法的核心逻辑极其简洁（dxt.cpp 第 3649-3658 行的宏展开）：

// 复数乘法核心
double a_re = dataA[j], a_im = dataA[j + 1];
double b_re = dataB[j], b_im = dataB[j + 1];
if (conjB) b_im = -b_im; // 编译时决定，不产生运行时分支
double c_re = a_re * b_re – a_im * b_im;
double c_im = a_re * b_im + a_im * b_re;
dataC[j] = (T)c_re;
dataC[j + 1] = (T)c_im;

注意乘法过程中用的是 double 精度——即使输入是 float，中间结果也用 double 计算再转回 float。这是为了避免复数乘法中的精度损失，尤其是当两个浮点数的乘积差值很小时（即 a_re*b_re 和 a_im*b_im 量级接近时），float 精度不够会产生显著误差。

七、频域滤波器设计：从低通到陷波

有了 DFT 和卷积定理作为工具，我们终于可以进入频域滤波器设计了。频域滤波的核心思路很简单：构造一个与频谱同尺寸的滤波器矩阵 H(u,v)，和图像的频谱 F(u,v) 逐元素相乘，再做逆 DFT 变回空间域。不同的 H(u,v) 设计，就对应不同的滤波效果。

7.1 低通滤波器：保留低频，去除高频

低通滤波器让低频分量通过、阻挡高频分量，效果类似于空间域的高斯模糊——图像变得平滑，边缘被弱化。

理想低通滤波器是最简单的设计：以频谱中心为圆心，半径 D₀ 以内的频率全部通过（H=1），以外全部阻止（H=0）。

// 理想低通滤波器
void idealLowPassFilter(Mat& filter, int rows, int cols, float D0) {
filter = Mat::zeros(rows, cols, CV_32F);
int cx = cols / 2, cy = rows / 2;
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
float d = sqrt((float)((i – cy) * (i – cy) + (j – cx) * (j – cx)));
if (d <= D0)
filter.at<float>(i, j) = 1.0f;
}
}
}

但理想低通有一个严重的问题——振铃效应（Ringing Effect）。因为 H(u,v) 在截止频率处有一个阶跃不连续点，对应到空间域上就是一个 sinc 函数，它的旁瓣会在图像边缘附近产生振荡伪影。

巴特沃兹低通滤波器通过引入平滑过渡来解决振铃问题：

// 巴特沃兹低通滤波器（n阶）
void butterworthLowPassFilter(Mat& filter, int rows, int cols, float D0, int n) {
filter = Mat(rows, cols, CV_32F);
int cx = cols / 2, cy = rows / 2;
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
float d = sqrt((float)((i – cy) * (i – cy) + (j – cx) * (j – cx)));
filter.at<float>(i, j) = 1.0f / (1.0f + pow(d / D0, 2.0f * n));
}
}
}

巴特沃兹滤波器的阶数 n 控制过渡带的陡峭程度：n 越大越接近理想滤波器（但振铃也会越明显），n=1 或 n=2 通常是实践中的最佳折中。

高斯低通滤波器的过渡更加平滑：

// 高斯低通滤波器
void gaussianLowPassFilter(Mat& filter, int rows, int cols, float D0) {
filter = Mat(rows, cols, CV_32F);
int cx = cols / 2, cy = rows / 2;
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
float d = sqrt((float)((i – cy) * (i – cy) + (j – cx) * (j – cx)));
filter.at<float>(i, j) = exp(-(d * d) / (2.0f * D0 * D0));
}
}
}

高斯滤波器完全没有振铃效应——因为高斯函数的傅里叶变换还是高斯函数，没有旁瓣。代价是过渡带相对较宽，不如巴特沃兹那么锐利地截断。

7.2 高通滤波器：提取边缘与细节

高通滤波器是低通的"互补"：H_hp(u,v) = 1 – H_lp(u,v)。它阻挡低频、保留高频，效果类似于空间域的拉普拉斯锐化——突出边缘和纹理细节。

// 从低通滤波器构造高通滤波器
Mat highPassFilter = Mat::ones(rows, cols, CV_32F) – lowPassFilter;

7.3 带通与带阻滤波器

带通滤波器只让某个频率范围内的信号通过，带阻滤波器则相反——阻止某个频率范围。带阻滤波器可以用两个不同半径的低通和高通组合得到，在消除特定频率范围的噪声时特别有用。

7.4 陷波滤波器：精准消除特定频率

陷波滤波器（Notch Filter）是频域去周期噪声的核心武器。它在频谱图上的特定位置"挖坑"——把那些对应周期噪声的频率分量压制为零（或接近零），而不影响其他频率。

由于频谱具有共轭对称性，陷波滤波器的零点必须成对出现：如果在 (u₀, v₀) 处挖一个坑，那么在 (-u₀, -v₀) 处也必须挖一个对称的坑，否则滤波结果会出现复数值。

巴特沃兹陷波滤波器的设计公式是：

H(u,v) = 1 / (1 + (D0² / (D1(u,v) · D2(u,v)))^n)

其中 D1 和 D2 分别是到两个对称陷波点的距离。

用 OpenCV 实现一个巴特沃兹陷波滤波器：

// 巴特沃兹陷波滤波器
// notchCenters: 陷波点坐标（频谱中心为原点）
// D0: 陷波半径
// n: 滤波器阶数
void butterworthNotchFilter(Mat& filter, int rows, int cols,
const vector<Point>& notchCenters,
float D0, int n) {
filter = Mat::ones(rows, cols, CV_32F);
int cx = cols / 2, cy = rows / 2;

for (const auto& center : notchCenters) {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
// 到正向陷波点的距离
float d1 = sqrt(pow(i – cy – center.y, 2) +
pow(j – cx – center.x, 2));
// 到对称陷波点的距离
float d2 = sqrt(pow(i – cy + center.y, 2) +
pow(j – cx + center.x, 2));

// 避免除零
if (d1 < 1e-5) d1 = 1e-5;
if (d2 < 1e-5) d2 = 1e-5;

float h = 1.0f / (1.0f + pow(D0 * D0 / (d1 * d2),
(float)n));
filter.at<float>(i, j) *= h;
}
}
}
}

陷波半径 D0 的选择至关重要：太小了，噪声消除不干净；太大了，会误伤周围的有用频率分量。通常的做法是先可视化频谱图，定位噪声亮点的位置和大小，再设置一个刚好覆盖亮点的半径。

八、实战：频域去周期性条纹噪声

现在我们把前面所有知识串联起来，实现一个完整的频域去周期噪声流程。假设我们有一张被水平条纹干扰的图像——这种干扰在工业检测、扫描文档等场景中很常见。

8.1 完整实现代码

#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>

using namespace cv;
using namespace std;

// 将频谱的四个象限交换到正确位置（低频移到中心）
void fftShift(Mat& mag) {
int cx = mag.cols / 2;
int cy = mag.rows / 2;
Mat q0(mag, Rect(0, 0, cx, cy));
Mat q1(mag, Rect(cx, 0, cx, cy));
Mat q2(mag, Rect(0, cy, cx, cy));
Mat q3(mag, Rect(cx, cy, cx, cy));
Mat tmp;
q0.copyTo(tmp); q3.copyTo(q0); tmp.copyTo(q3);
q1.copyTo(tmp); q2.copyTo(q1); tmp.copyTo(q2);
}

// 计算频谱幅度（用于可视化）
Mat computeMagnitudeSpectrum(const Mat& complexImg) {
Mat planes[2];
split(complexImg, planes);
Mat mag;
magnitude(planes[0], planes[1], mag);
mag += Scalar::all(1);
log(mag, mag);
normalize(mag, mag, 0, 255, NORM_MINMAX);
mag.convertTo(mag, CV_8U);
return mag;
}

// 巴特沃兹陷波滤波器
Mat createNotchFilter(int rows, int cols,
const vector<Point>& notchCenters,
float D0, int order) {
Mat filter = Mat::ones(rows, cols, CV_32F);
int cx = cols / 2, cy = rows / 2;

for (const auto& center : notchCenters) {
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
float d1 = sqrt(pow(i – cy – center.y, 2) +
pow(j – cx – center.x, 2));
float d2 = sqrt(pow(i – cy + center.y, 2) +
pow(j – cx + center.x, 2));
if (d1 < 1e-5f) d1 = 1e-5f;
if (d2 < 1e-5f) d2 = 1e-5f;
float h = 1.0f / (1.0f + pow(D0 * D0 / (d1 * d2),
(float)order));
filter.at<float>(i, j) *= h;
}
}
}
return filter;
}

int main() {
// 1. 读取灰度图像
Mat img = imread("striped_image.png", IMREAD_GRAYSCALE);
if (img.empty()) {
cerr << "无法读取图像" << endl;
return -1;
}

// 2. 扩展到最优DFT尺寸（关键步骤！）
int m = getOptimalDFTSize(img.rows);
int n = getOptimalDFTSize(img.cols);
Mat padded;
copyMakeBorder(img, padded, 0, m – img.rows, 0, n – img.cols,
BORDER_CONSTANT, Scalar::all(0));

// 3. 转为浮点数，构建双通道复数矩阵
Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
Mat complexImg;
merge(planes, 2, complexImg);

// 4. 正向DFT（使用DFT_COMPLEX_OUTPUT获取双通道复数输出）
dft(complexImg, complexImg);

// 5. 将低频分量移到频谱中心（方便分析和设计滤波器）
fftShift(complexImg);

// 6. 可视化频谱，定位噪声亮点
Mat spectrumBefore = computeMagnitudeSpectrum(complexImg);
imwrite("spectrum_before.png", spectrumBefore);

// 7. 定义陷波点（根据频谱分析结果设置）
// 假设分析发现水平条纹噪声在频谱中心上下方各有一组亮点
vector<Point> notchCenters;
// 这些坐标需要根据实际频谱分析确定
notchCenters.push_back(Point(0, 30)); // 上方亮点
notchCenters.push_back(Point(0, 60)); // 更远的谐波
notchCenters.push_back(Point(0, 90)); // 第三次谐波

// 8. 创建陷波滤波器
Mat notchFilter = createNotchFilter(padded.rows, padded.cols,
notchCenters, 10.0f, 2);

// 9. 将实数滤波器扩展为双通道（实部=滤波器，虚部=0）
Mat filterPlanes[] = {notchFilter, Mat::zeros(notchFilter.size(), CV_32F)};
Mat filter2ch;
merge(filterPlanes, 2, filter2ch);

// 10. 频域乘法（应用滤波器）
mulSpectrums(complexImg, filter2ch, complexImg, 0);

// 11. 可视化滤波后的频谱
Mat spectrumAfter = computeMagnitudeSpectrum(complexImg);
imwrite("spectrum_after.png", spectrumAfter);

// 12. 频谱中心移回原位
fftShift(complexImg);

// 13. 逆DFT
Mat result;
idft(complexImg, result, DFT_SCALE | DFT_REAL_OUTPUT);

// 14. 裁剪回原始尺寸，转为8位输出
result = result(Rect(0, 0, img.cols, img.rows));
result.convertTo(result, CV_8U);

imwrite("result.png", result);
cout << "处理完成！" << endl;
return 0;
}

8.2 关键步骤详解

步骤 2：getOptimalDFTSize 零填充。这一步直接调用我们前面分析过的 getOptimalDFTSize()，确保 DFT 在高效的 2^p·3^q·5^r 尺寸上执行。如果跳过这一步，对于某些尺寸（比如素数大小的图像），DFT 可能慢几十倍。

步骤 4：dft() 调用。这里直接传入双通道复数矩阵，所以 dft() 内部走的是 FwdComplexToComplex 模式，不会经过 RealDFT 路径。如果想利用 RealDFT 的 N/2 优化，可以传入单通道实数矩阵并加上 DFT_COMPLEX_OUTPUT 标志：dft(floatImg, complexOut, DFT_COMPLEX_OUTPUT)。

步骤 5：fftShift。OpenCV 的 dft() 输出的频谱中，直流分量（DC，即零频率）在左上角，而不是中心。fftShift() 把四个象限对调，让低频移到中心——这样频谱图"从中心向外频率递增"，方便直观分析和设计滤波器。

步骤 6：频谱分析。对频谱取幅度、做对数变换（因为频谱的动态范围通常很大，对数压缩后才能看到细节），再归一化到 0-255 显示。这一步产生的频谱图是设计陷波滤波器的依据——那些远离中心的孤立亮点就是周期噪声的"指纹"。

步骤 7：定位噪声亮点。水平条纹在频谱上表现为垂直方向上的亮点（因为空间域的水平条纹对应频域垂直轴上的频率分量）。亮点的距中心远近代表条纹的空间频率——距离越远，条纹越密集。

步骤 8：陷波滤波器参数。D0=10 表示陷波半径为 10 个像素，order=2 是二阶巴特沃兹。这两个参数需要根据实际亮点的大小和扩散程度调整。如果噪声亮点很尖锐（集中在 1-2 个像素），D0=5 就够了；如果亮点有一定扩散，可能需要 D0=15-20。

步骤 10：频域乘法。这里用 mulSpectrums() 做复数乘法。滤波器矩阵构造为双通道复数格式（实部=滤波器系数，虚部=0），mulSpectrums() 会按照复数乘法规则 (a+bj)·(h+0j) = ah + bhj 正确缩放频谱的实部和虚部。注意不能用 cv::multiply() 替代——multiply() 对双通道 Mat 做的是逐通道实数乘法而非复数乘法，虽然对纯实数滤波器碰巧结果一致，但语义不正确，遇到复数滤波器时会出错。

步骤 13：逆DFT。DFT_SCALE 标志让 idft() 自动除以 M×N（正变换和逆变换之间默认不做缩放，需要手动指定一方缩放）。DFT_REAL_OUTPUT 标志告诉 idft() 直接输出实数单通道结果。

8.3 调参建议

在实际应用中，频域去周期噪声的效果高度依赖参数选择。几个关键的调参经验：

陷波点定位：不要手动估计亮点坐标，而是用代码自动检测。最简单的方法是：先对频谱幅度图做阈值处理（排除中心的低频区域），然后用 findContours 或 findNonZero 定位高亮区域的质心。

陷波半径 D0：先从小半径开始（D0=5），逐步增大直到噪声完全消除。如果增大到 D0=30 还有残留，说明陷波点的坐标定位不准，需要重新检查频谱图。

阶数 n：n=2 是大多数场景的最佳平衡。n=1 过渡太平滑，可能不够力度；n=4 或更高则接近理想滤波，边界处可能出现伪影。

多次谐波：真实的周期噪声通常不是纯正弦波，它会在频谱上产生基频和多次谐波。需要同时在基频和谐波位置设置陷波点，否则只消基频不消谐波，效果会打折扣。

九、总结与工程建议

回顾全文，我们从 DFT 的数学定义出发，走过了一条完整的技术链路：

知识点

源码位置

核心要点

从工程角度，几条建议值得记住：

第一，永远使用 getOptimalDFTSize()。一个 641×481 的图像，DFT 可能比 640×480 慢几十倍。多花几微秒做零填充，换来几毫秒甚至几十毫秒的 DFT 加速，收益巨大。

第二，善用 nonzeroRows 参数。做频域卷积时，零填充后的图像大部分行都是零值，传入 nonzeroRows 可以跳过这些行的无用计算，实测能节省 30-50% 的时间。

第三，实数图像用单通道输入。传入单通道 Mat 让 OpenCV 走 RealDFT 路径，比自己手动构造双通道复数矩阵再调 dft() 快近 2 倍——前提是你不需要移频到中心后做可视化和手动滤波。如果需要移频操作，就用 DFT_COMPLEX_OUTPUT 获得双通道输出。

第四，频域滤波 vs 空间域滤波的选择。对于小核（3×3、5×5、7×7），空间域滤波（filter2D()、GaussianBlur() 等）更快；核大小超过约 30×30 时，频域方法开始有优势；当核超过 100×100，频域方法的速度优势会越来越明显。而对于周期噪声这种只有频域才能精准处理的问题，频域滤波是唯一选择。

OpenCV 的 dft() 只有一个函数调用那么简单，但背后是几十年 FFT 算法研究和工程优化的结晶——从 Cooley-Tukey 的分治思想到混合基蝶形运算，从位反转查找表到 SSE3 SIMD 指令，从 RealDFT 的 N/2 技巧到 CCS 紧凑格式，每一层优化都在平衡着性能、精度和通用性。理解这些底层细节，才能在实际项目中做出正确的工程决策，而不是把 dft() 当成一个黑盒子来调用。