OpenCV 傅里叶变换实战指南：Numpy 与 cv.dft 频谱分析、频域滤波及 FFT 性能优化

本文基于 OpenCV 官方 Python 教程 py_fourier_transform.markdown 展开，系统讲解如何在 OpenCV 中计算图像的二维离散傅里叶变换（2D DFT）：既覆盖 np.fft.fft2 / np.fft.ifft2 的 Numpy 路线，也覆盖 cv.dft() / cv.idft() 的 OpenCV 路线，并深入源码解释 cv.getOptimalDFTSize() 背后的最优尺寸表与质因数分解机制，最终带你理解"Laplacian 为什么是高通滤波器"这一经典问题的频域视角。读完本文，你将能够独立完成图像频谱可视化、频域低通/高通滤波、以及 DFT 计算的性能调优。

理论背景：图像的频率域表示

傅里叶变换用于分析各种滤波器的频率特性。对图像而言，使用 二维离散傅里叶变换（2D DFT） 来获取频率域表示，实际计算时采用快速算法 FFT（Fast Fourier Transform）。其基本思想可以概括为：

• 对正弦信号 $x(t) = A \sin(2 \pi f t)$，其频谱在 $f$ 处出现一个尖峰；
• 离散化后频谱在 $[-\pi, \pi]$ 或 $[0, 2\pi]$（即 N 点 DFT 的 $[0, N]$）内呈周期性；
• 图像可以看作在两个方向上采样的信号，因此在 X、Y 两个方向各做一次傅里叶变换，就得到图像的频率表示。

一个更直观的理解是：幅度变化快的部分对应高频，变化慢的部分对应低频。在图像中，幅度的剧烈变化发生在边缘点和噪声处，因此边缘和噪声属于图像的高频成分；幅度变化平缓的区域则属于低频成分。这也是后文频域滤波操作的立论基础。

用 Numpy 求傅里叶变换

Numpy 提供了 FFT 包来实现这一功能。np.fft.fft2() 返回复数数组形式的频谱，其参数行为如下：

参数	说明
第一个参数	输入图像，须为灰度图
第二个参数（可选）	指定输出数组尺寸。若大于输入尺寸，输入先被零填充再计算 FFT；若小于输入尺寸，输入会被裁剪；不传则输出与输入同尺寸

拿到频谱后，零频率分量（直流分量 DC）位于左上角。若希望将其移到中心便于分析，需在两个方向上平移 $N/2$ 个位置，这正是 np.fft.fftshift() 的作用。之后即可计算幅度谱：

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
assert img is not None, "file could not be read, check with os.path.exists()"
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
magnitude_spectrum = 20*np.log(np.abs(fshift))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

可以看到图像中心有更亮的区域，说明低频成分占主导——图像的大部分数据集中在频谱的低频区域。

频域高通滤波与振铃效应

得到频谱后，就可以在频域做操作，例如高通滤波后再求逆 DFT 重建图像。具体做法：用一个 60x60 的矩形窗将低频区域置零，然后用 np.fft.ifftshift() 把直流分量移回左上角，最后用 np.fft.ifft2() 求逆 FFT。结果仍是复数，取其绝对值（实部）即可：

rows, cols = img.shape
crow, ccol = rows//2, cols//2
fshift[crow-30:crow+31, ccol-30:ccol-31] = 0
f_ishift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f_ishift)
img_back = np.real(img_back)

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(img_back, cmap = 'gray')
plt.title('Image after HPF'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back)
plt.title('Result in JET'), plt.xticks([]), plt.yticks([])

plt.show()

该结果印证了两点结论：

1. 高通滤波本质上是一种边缘检测操作（与图像梯度章节的结论一致）；
2. 图像的大部分数据确实存在于频谱的低频区域。

仔细观察 JET 伪彩色结果中的波纹状结构（原文档以红色箭头标注），这就是所谓的振铃效应（ringing effects）。它由掩膜所用的矩形窗引起：矩形窗在频域中对应 sinc 形状的响应，从而产生波纹。因此实际滤波中不宜使用矩形窗，更好的选择是高斯窗。

用 OpenCV 求傅里叶变换

OpenCV 提供了 cv.dft() 和 cv.idft() 两个函数。与 Numpy 返回复数数组不同，OpenCV 返回两通道数组：第一通道是实部，第二通道是虚部；输入图像需要先转换为 np.float32：

import numpy as np
import cv2 as cv
from matplotlib import pyplot as plt

img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
assert img is not None, "file could not be read, check with os.path.exists()"

dft = cv.dft(np.float32(img),flags = cv.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

magnitude_spectrum = 20*np.log(cv.magnitude(dft_shift[:,:,0],dft_shift[:,:,1]))

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(magnitude_spectrum, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

提示：也可以用 cv.cartToPolar() 一次性同时获得幅度和相位，替代 cv.magnitude() 的手动取通道方式。

频域低通滤波（LPF）

上一节在 Numpy 中做的是高通滤波（去低频），这里用 OpenCV 做低通滤波（去高频）——效果相当于对图像做模糊。做法是构造一个"低频为 1、高频为 0"的掩膜：

rows, cols = img.shape
crow, ccol = rows//2, cols//2

# create a mask first, center square is 1, remaining all zeros
mask = np.zeros((rows,cols,2),np.uint8)
mask[crow-30:crow+30, ccol-30:ccol+30] = 1

# apply mask and inverse DFT
fshift = dft_shift*mask
f_ishift = np.fft.ifftshift(fshift)
img_back = cv.idft(f_ishift)
img_back = cv.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_back, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.show()

关于掩膜的两个实现细节值得注意：

• 掩膜创建为 (rows, cols, 2) 的三通道数组，与 dft_shift 的两通道（实部/虚部）结构直接相乘；
• 重建时调用 cv.magnitude() 对 idft 结果的实/虚两通道求模，还原为单通道灰度图。

官方教程同时指出：OpenCV 的 cv.dft() 与 cv.idft() 比 Numpy 对应函数更快，但 Numpy 的 API 更友好易用，性能差异详见下一节。

DFT 性能优化：cv.getOptimalDFTSize()

DFT 的计算性能与数组尺寸密切相关：尺寸是 2 的幂时最快；尺寸可分解为 2、3、5 的乘积时也较高效。因此若在意性能，可以在求 DFT 之前通过零填充把数组调整为"最优尺寸"。OpenCV 中需要手动填充零，而 Numpy 只要指定新的 FFT 尺寸即可自动填充。

查找最优尺寸的函数是 cv.getOptimalDFTSize()，它同样适用于 np.fft.fft2()。教程用 IPython 的 %timeit 实测（输入 messi5.jpg，原始尺寸 342x548）：

In [15]: img = cv.imread('messi5.jpg', cv.IMREAD_GRAYSCALE)
In [16]: assert img is not None, "file could not be read, check with os.path.exists()"
In [17]: rows,cols = img.shape
In [18]: print("{} {}".format(rows,cols))
342 548

In [19]: nrows = cv.getOptimalDFTSize(rows)
In [20]: ncols = cv.getOptimalDFTSize(cols)
In [21]: print("{} {}".format(nrows,ncols))
360 576

尺寸 (342, 548) 被调整为 (360, 576)。接下来用两种等价方式补零：

# 方式一：新建全零大数组并拷贝
nimg = np.zeros((nrows,ncols))
nimg[:rows,:cols] = img

或者用 cv.copyMakeBorder()：

# 方式二：copyMakeBorder
right = ncols - cols
bottom = nrows - rows
bordertype = cv.BORDER_CONSTANT
nimg = cv.copyMakeBorder(img,0,bottom,0,right,bordertype, value = 0)

教程给出的实测数据（具体数值随硬件环境而异）：

操作	耗时
%timeit np.fft.fft2(img)	40.9 ms/loop
%timeit np.fft.fft2(img,[nrows,ncols])	10.4 ms/loop（约 4 倍加速）
%timeit cv.dft(np.float32(img),flags=cv.DFT_COMPLEX_OUTPUT)	13.5 ms/loop
%timeit cv.dft(np.float32(nimg),flags=cv.DFT_COMPLEX_OUTPUT)	3.11 ms/loop（约 4 倍加速）

两个结论与教程一致：最优尺寸优化可带来约 4 倍加速；OpenCV 函数比 Numpy 快约 3 倍。逆 FFT 的测试则留作练习。

源码级佐证：最优尺寸表与质因数分解

这一节的行为在 OpenCV 源码中可以完整印证。

cv::getOptimalDFTSize() 的实现位于 dxt.cpp：它并不现场计算，而是对一张预编译的 optimalDFTSizeTab 表（见 dxt.cpp）做二分查找，返回不小于 size0 的最小表项；该表中的每一项都是 2、3、5 的幂次乘积（如 342 → 360 = 2³×3²×5，548 → 576 = 2⁶×3²，与教程输出一致）。

而"为什么偏偏是 2、3、5"，答案在同文件的 DFTFactorize()（dxt.cpp）：该静态函数把变换尺寸 n 分解为因子序列，优先提取 2，然后从 3 开始只按奇数步长试探 3、5 等小质因子，配合内部的 Cooley-Tukey 类分治结构实现高效计算。API 文档也在 core.hpp 中明确说明：dft 支持任意尺寸，但只有可分解为小质数（当前实现为 2、3、5）乘积的数组才能被高效处理。DFT 相关测试用例见 test_dxt.cpp。

为什么 Laplacian 是高通滤波器？

这是一个经典的频域分析问题：为什么 Laplacian（以及 Sobel）是高通滤波器（HPF），而均值/高斯核是低通滤波器（LPF）？方法是直接对这些核做较大尺寸的 FFT 并观察幅度谱——核在哪个频率区域"能量大"（通过），哪个区域"能量小"（阻挡），它就在图像上起相应的滤波作用：

import cv2 as cv
import numpy as np
from matplotlib import pyplot as plt

# simple averaging filter without scaling parameter
mean_filter = np.ones((3,3))

# creating a gaussian filter
x = cv.getGaussianKernel(5,10)
gaussian = x*x.T

# different edge detecting filters
# scharr in x-direction
scharr = np.array([[-3, 0, 3],
                   [-10,0,10],
                   [-3, 0, 3]])
# sobel in x direction
sobel_x= np.array([[-1, 0, 1],
                   [-2, 0, 2],
                   [-1, 0, 1]])
# sobel in y direction
sobel_y= np.array([[-1,-2,-1],
                   [0, 0, 0],
                   [1, 2, 1]])
# laplacian
laplacian=np.array([[0, 1, 0],
                    [1,-4, 1],
                    [0, 1, 0]])

filters = [mean_filter, gaussian, laplacian, sobel_x, sobel_y, scharr]
filter_name = ['mean_filter', 'gaussian','laplacian', 'sobel_x', \
                'sobel_y', 'scharr_x']
fft_filters = [np.fft.fft2(x) for x in filters]
fft_shift = [np.fft.fftshift(y) for y in fft_filters]
mag_spectrum = [np.log(np.abs(z)+1) for z in fft_shift]

for i in range(6):
    plt.subplot(2,3,i+1),plt.imshow(mag_spectrum[i],cmap = 'gray')
    plt.title(filter_name[i]), plt.xticks([]), plt.yticks([])

plt.show()

从图中可以直接读出每个核阻挡和通过的频率区域：mean_filter 与 gaussian 的能量集中在频谱中心（低频通过，故为 LPF）；laplacian、sobel_x、sobel_y、scharr 的中心能量低而外围能量高（低频被阻断，故为 HPF）。这也与 Laplacian 核系数之和为 0（对常数/低频信号响应为零）的直观解释互为印证。

小结与参考

本文的核心知识点与对应实现依据如下：

主题	关键点	依据
Numpy 路线	np.fft.fft2 / fftshift / ifftshift / ifft2，零填充由 size 参数隐式完成	原文档 Numpy 章节
OpenCV 路线	cv.dft 需先转 np.float32，输出为实部/虚部两通道，DFT_SCALE 需显式指定才会缩放	core.hpp
性能优化	cv.getOptimalDFTSize() 查预编译表（二分查找）+ 手动 copyMakeBorder 补零	dxt.cpp、dxt.cpp
滤波选型	矩形窗产生 sinc 响应与振铃伪影，频域滤波宜用高斯窗	原文档 ringing effects 段落
HPF/LPF 判定	对核本身做 FFT，看其幅度谱能量分布	原文档 Laplacian 章节

延伸阅读可参考仓库内的 C++ 示例 dft.cpp（DFT 卷积）与 Python 示例 dft.py（傅里叶象限重排）、deconvolution.py（Wiener 反卷积），API 声明集中在 core.hpp，核心实现在 dxt.cpp。