☰
离散傅里叶变换(DFT)在图像处理中的应用与频域滤波实战
2026/10/15 11:33:05 网站建设 项目流程

做数字图像处理的人,无论做什么方向,最后都会碰到同一个绕不开的概念——DFT,也就是离散傅里叶变换。我的数字图像处理笔记写到第四篇,前几篇基本都在空域打转:灰度变换、直方图、平滑锐化,这些操作直观,但很多问题在空域里很难一眼看穿。这一次我们把图像变换到频域,从另一个维度重新理解图像。这篇笔记我打算把DFT的概念、原理、实现和应用完整地串一遍,最后附带我踩过的几个坑,适合正在学数字图像处理的学生,也适合工作中需要做图像滤波或去噪的朋友。读完你可以独立完成图像的频谱分析、频域滤波和周期性噪声去除。

1. DFT是什么:从空间域到频域的思维转换

1.1 为什么图像处理绕不开DFT

DFT的全称是Discrete Fourier Transform,中文译为离散傅里叶变换。它是连续傅里叶变换在离散数据上的对应版本。图像本身是离散的像素点,所以图像处理中我们用DFT而不是连续傅里叶变换。

傅里叶变换的核心思想是:任何一个信号都可以分解为一系列不同频率的正弦波的叠加。图像可以视为一种二维信号,只不过自变量是空间坐标x和y,而不是时间t。所以图像里的“频率”是空间频率,表示灰度在空间上变化的快慢。灰度变化缓慢的平坦区域对应低频;灰度剧烈变化的边缘、细节、纹理对应高频。

这里有个类比可以帮你建立直觉:把图像想象成一块布料,低频相当于布料的底色和整体图案,高频相当于布面上细小的织纹和毛刺。DFT就是一台光谱仪,能够告诉你这块布料里包含了哪些“空间频率成分”以及每种成分有多强。一旦拿到这张频谱,你就可以有针对性地过滤掉不想要的成分,这就是频域滤波的基础。

为什么说图像处理绕不开DFT?因为在空域里,一个高斯滤波核只能作用于局部像素,无法很好地区分“周期性噪声”和“真实边缘”;而在频域里,周期噪声和真实边缘有着截然不同的分布。所以很多空域里需要反复试参数才能勉强处理的问题,换到频域只需看一眼频谱图就能定位。这是经验的直观感受,后面我会用例子演示。

1.2 数学定义与离散化的必要性

先看一维DFT的定义。设离散信号x[0], x[1], ..., x[N-1]长度为N,它的DFT定义为:

X[k] = Σ_{n=0}^{N-1} x[n] * exp(-j * 2π * k * n / N),k = 0, 1, ..., N-1

其中j是虚数单位,exp(-j2πk*n/N)可以理解为旋转因子。每一个X[k]都是复数,代表该信号在第k个频率分量上的幅度和相位。

二维DFT是它在一维基础上的直接推广。设图像f(x, y)尺寸为M×N,则:

F(u, v) = Σ_{x=0}^{M-1} Σ_{y=0}^{N-1} f(x, y) * exp(-j * 2π * (u*x/M + v*y/N))

u和v分别是水平方向和垂直方向的空间频率。F(u, v)同样是一个复数,包含了幅度和相位两部分信息。

为什么要专门强调“离散”?因为在模拟世界中,光照通过透镜成像,理论上可以用连续函数描述;但计算机能处理的只有有限个采样点。连续傅里叶变换需要积分,在计算机上无法精确计算。离散化等价于在时域和频域同时抽样,DFT把有限长度信号变换为有限长度频谱,这才是可计算的。换句话说,DFT是连续傅里叶变换在数字世界里的“可执行版本”。

这里有一个很多人困惑的点:DFT输出的长度等于输入信号长度。图像是M×N个点,DFT结果也是M×N个点。因为变换是可逆的,没有信息丢失,这也是我们后面能够用逆变换重建图像的前提。

1.3 幅度谱与相位谱到底在说什么

由于DFT输出是复数,我们可以把它拆成两种表示方式:实部虚部,或者幅度相位。工程上更常用的是幅度谱和相位谱。幅度谱表示每个频率成分的强度,计算公式为|F(u,v)|;相位谱表示每个频率成分的相位偏移,计算公式为arctan(虚部/实部)。

直观地看,幅度谱告诉我们图像中哪些频率占主导地位。自然图像的能量大多集中在低频区,所以频谱图上中心区域亮、四周暗。而相位谱包含的信息常常被忽略,但实际上它比幅度谱更能描述图像的结构特征。

做过实验就会印象深刻:把一张图的幅度谱保留、相位谱替换成另一张随机图,逆变换后图像几乎面目全非,只剩下一些模糊的色块;反过来,把相位谱保留、幅度谱全部置为1,逆变换后虽然丢失了原始的明暗对比,但图像的主要边缘和轮廓仍然清晰可辨。这让我在工程里养成了一个习惯:凡是要修改频域的系数,一定要先想清楚改的是幅度还是相位。大多数滤波器只改动幅度,尽量不碰相位,就是为了避免破坏图像的结构关系。

2. DFT在图像处理中的实际用途

2.1 频域滤波:一键实现低通与高通

频域滤波的思路很简单:先对图像做DFT,得到频谱F(u,v);再构造一个频域滤波器H(u,v),用H点乘F;最后对结果做逆DFT,得到滤波后的图像。整个流程可以写成一个公式:g = IDFT(F * H)。

这里的H是跟图像同尺寸的矩阵,取值通常在0到1之间。如果H在低频区域为1、高频区域接近0,那就是低通滤波器,效果相当于空域平滑,可以去除噪声但也会让边缘变模糊。如果H在低频区域接近0、高频区域为1,那就是高通滤波器,效果相当于锐化,可以提取边缘但也会放大噪声。

很多教科书会把低通、高通、带通、带阻等每种滤波器都画一遍频谱图,实际用起来其实没那么复杂:你先看频谱中哪些位置需要保留、哪些需要削减,然后构造对应的掩码矩阵即可。在空域里需要设计卷积核、考虑边界、调整尺寸,在频域里这些都简化为一次矩阵乘法,这是频域处理最大的工程优势。

来看一个高斯低通滤波的例子。高斯低通滤波器的传递函数可以用下面这段代码构造,其中cutoff是截止频率,控制保留低频的程度:

def gaussian_lowpass(shape, cutoff): rows, cols = shape crow, ccol = rows // 2, cols // 2 y = np.arange(rows).reshape(-1, 1) x = np.arange(cols).reshape(1, -1) d = np.sqrt((y - crow)**2 + (x - ccol)**2) return np.exp(-0.5 * (d / cutoff)**2)

cutoff设得越小,保留的低频越少,图像越模糊;cutoff设得越大,滤波效果越接近原图。我习惯把cutoff设为图像短边尺寸的5%到10%作为一个起点,再根据实际效果微调。如果你发现滤波后的图像出现振铃现象,多半是截止频率设置过低或滤波器边缘太陡,后面第5节会详细讨论。

2.2 去周期性噪声的利器

周期噪声是DFT最拿手的场景。实际项目中,扫描仪的影线、相机传感器上的规律坏点、显示器拍摄时的摩尔纹,都会在图像中形成周期性纹理。这种噪声的可怕之处在于:它们在空域中与真实内容混在一起,普通的均值滤波去不掉,中值滤波也只能缓解,而且很容易把图像细节都抹掉。

但在频域里,周期性噪声的表现非常独特:它会集中在频谱上少数几个离散的亮点(也称为频点)上。因为一个在空间域呈正弦变化的条纹,其傅里叶变换就是一对位于特定频率处的冲激。也就是说,只要在频谱图上找到这些亮点,把它们周围一个小区域置零,就能把对应频率的噪声去掉,而几乎不损伤其他频段的内容。

这种方法叫陷波滤波(notch filtering)。实际工程里要比教科书麻烦一点:噪声频点往往不是严格的一个像素点,而是呈一个小十字或小圆斑的形状,所以通常需要把亮点附近一个小半径的区域都置零,而不是只去掉一个点,否则会出现残留的条纹。这个细节我建议你在自己动手试的时候一定留意,我在第四部分会给出完整代码。

2.3 图像压缩与特征表达的频域视角

除了滤波,DFT在图像压缩中也有很强的关联。经典的JPEG压缩虽然不是直接用DFT,而是用DCT(离散余弦变换),但DCT可以看作DFT的实部形式,核心思路完全一致:把图像变换到频域,能量集中到少数低频系数上,然后对高频系数进行量化甚至置零来压缩数据。

从频域的视角也能解释很多空域现象。比如图像缩放后的模糊感,本质是重采样时混叠了高频分量;图像锐化,本质上是在增强高频分量的权重;图像插值,本质上是在频域里补充缺失的采样。理解DFT之后,再回头看空域的很多操作,会有一种“原来如此”的感觉。

此外,DFT的移位性、卷积定理也常被用于模板匹配和快速计算。卷积定理说:空域的卷积等于频域的乘法。这个定理的实际价值就是“用FFT加速大核卷积”。当卷积核特别大时,空域卷积的计算量随核尺寸平方增长,而频域方案几乎不受核尺寸影响。这个特性是后面第三节里性能对比的理论基础。

3. 从公式到代码:DFT的实现完全拆解

3.1 FFT:让DFT从理论走向工程

直接按照上一节的定义来实现DFT,计算量是O(N^2)的。如果你要处理一张512×512的图像,直接做二维DFT需要几百亿次复数乘加运算,任何主流PC跑起来都会非常吃力。这也是为什么早期数字图像处理停留在理论阶段,直到快速傅里叶变换(FFT)出现才真正进入工程。

FFT的核心思想是分治。它利用旋转因子的周期性、共轭对称性,把一个N点DFT不断拆分为两个N/2点DFT,从而把计算复杂度从O(N^2)降到O(N log N)。这个差距是数量级的:同样的512×512图像,直接二维DFT和FFT的计算量相差大约是上万倍,从不可能到毫秒级完成。今天所有numpy、OpenCV里的傅里叶变换实现,底层都是FFT,只是封装得非常好用,你不需要自己写蝴蝶运算,但要理解它解决了什么问题。

在使用FFT接口时,我建议你养成两个习惯:一是在内存足够的情况下优先选择64位的浮点类型,不要用8位整型,因为傅里叶结果大量介于小数之间,整型计算会带来明显误差;二是不要重复计算,如果同一张图像要做多次滤波,把DFT结果缓存起来,只更新滤波掩码再逆变换,效率会快很多。

3.2 用numpy快速上手二维DFT与频谱可视化

python里做DFT,最常用的就是numpy的fft模块。下面是一段完整的频谱分析代码,注释已经写在里面:

import cv2 import numpy as np import matplotlib.pyplot as plt img = cv2.imread('input.jpg', cv2.IMREAD_GRAYSCALE) print(img.shape) # 二维DFT,默认变换轴上后两个维度 f = np.fft.fft2(img) # 平移频谱,把低频分量移到中心 fshift = np.fft.fftshift(f) # 幅度谱并做log压缩 magnitude = np.abs(fshift) magnitude_log = np.log(1 + magnitude) # 可视化 plt.figure(figsize=(10, 4)) plt.subplot(1, 2, 1), plt.imshow(img, cmap='gray'), plt.title('Original') plt.subplot(1, 2, 2), plt.imshow(magnitude_log, cmap='gray'), plt.title('Spectrum(log)') plt.show()

这里有两个地方新手几乎都会踩坑。

第一个坑是fftshift漏掉。numpy的fft2输出结果里,零频分量(DC分量)位于数组的左上角(0,0)位置,而不是中心。如果不做fftshift,你会看到频谱的四个角很亮,中间是黑的,很多人第一次看到会觉得“怎么我的频谱是反着的”。fftshift做的事情就是把数组重新排列:让低频出现在中心,高频围绕在四周,这和教科书里画的频谱图才是一致的。

第二个坑是幅度谱的动态范围。图像的傅里叶幅度谱动态范围巨大,最大的DC值可能达到几百万甚至上亿,而大部分高频分量只有几百甚至几十。如果不加处理直接显示,整张图除了中心一个白点外全是黑的。所以要取log(1 + magnitude),用对数压缩把动态范围缩小,才能看到频谱图里的丰富细节。注意这里加1是为了避免log0,同时保证幅度为0时输出也是0。

3.3 逆变换与无损重建的细节

把频域的图变回空间域,用ifft2即可:

# 如果之前做过fftshift,逆变换前先ifftshift f_ishift = np.fft.ifftshift(fshift) # 逆DFT img_recon = np.fft.ifft2(f_ishift) # 取实部 img_recon = np.real(img_recon) # 确认重建误差 err = np.mean(np.abs(img - img_recon)) print(err)

逆变换后数组是复数,但虚部理论上应该为0。实际计算中会残留10^-12量级的虚部,这是浮点运算的累积误差,正常现象,直接取real即可。如果你发现虚部很大,要检查一下是不是某一步把实部和虚部弄混了。

关于归一化,numpy的ifft2已经内置了1/(M*N)缩放因子,所以fft2之后再ifft2可以无损恢复原图。但如果你自己写DFT实现,就要格外注意缩放因子放在正变换还是逆变换侧。常见约定有三种,不管哪种,只要保证正逆变换的乘积关系即可。我的建议是直接用库封装好的函数,少去手动实现,因为边界条件和性能优化都很容易出问题。

还需要注意:cv2.dft和np.fft.fft2的输出布局略有不同。cv2.dft输出的是多通道数组,通道0为实部、通道1为虚部,类型为float32;而np.fft.fft2输出复数数组。在切换API时不要直接拿来就乘,先检查shape和dtype。有一次我把cv2.dft的输出直接拿去np.fft.fftshift,结果因为多通道维度的问题导致所有系数错位,排查了半天,后来才意识到是两个库的表示方式不同。

4. 实战:用DFT干掉周期性条纹噪声

4.1 先看频谱,再定方案

假设你拿到了这样一张工业品表面图片:本底是细密的编织纹理,但上面叠了一层肉眼可见的横向条纹。这层条纹的间距很规律,属于周期性干扰。想用中值滤波或者高斯滤波去掉它,要么留下残影,要么把编织纹理也破坏了。这时候就是DFT上场的时候。

第一步不是直接写滤波器,而是先对图像做频谱分析。用第三节里的log幅度谱可视化,你会看到在频谱中心垂直方向(对应空间域的水平条纹)上出现一对对称的亮点。只要找到这对亮点的坐标,噪声频率就定位了。为什么要看频谱再定方案?因为不做分析就去调滤波参数,相当于蒙着眼睛修电路,大概率会把别的频段也误伤。而频谱图会把噪声位置清清楚楚地标出来,后续设计陷波器就有了依据。

要强调的是,亮点的位置不是固定的。如果条纹在空间域里很细、间隔很小,亮点就会离频谱中心更远;如果条纹很粗、间隔很大,亮点离中心就更近。也就是说,噪声频率和空间周期是倒数关系,距离中心越远对应频率越高。这个关系在动手时非常实用:你可以通过频谱图直接估算条纹的周期,而不必在空间域拿尺子量。

4.2 陷波滤波的完整实现

拿到噪声亮点坐标后,我通常是这么处理的:构造一个全为1的掩码图像,把要抑制的频点周围一小块区域置为0,然后频谱乘以掩码,再做逆变换。下面的代码演示了完整流程:

h, w = img.shape crow, ccol = h // 2, w // 2 # 假设通过频谱分析发现两个噪声亮点坐标如下 noise_points = [(crow - 52, ccol), (crow + 52, ccol)] # 构造掩码,默认全1(保留所有频率) mask = np.ones((h, w), dtype=np.float32) # 把噪声亮点周围半径5的区域置0 radius = 5 for cy, cx in noise_points: cv2.circle(mask, (cx, cy), radius, 0, -1) # 原频谱(在第三节中已经拿到fshift) fshift_filtered = fshift * mask # 逆变换回空间域 f_ishift = np.fft.ifftshift(fshift_filtered) img_clean = np.fft.ifft2(f_ishift) img_clean = np.real(img_clean)

这里有一个我在实际项目里踩过的坑:cv2.circle的圆心参数是(x, y),而numpy数组索引是(row, col),两者顺序相反。我第一次写成cv2.circle(mask, (cy, cx), ...)时,噪声点被标到了镜像位置,滤波器没覆盖真正的噪点,条纹还在。所以上面的代码我特意写成(cx, cy),提醒你注意坐标顺序。

半径的选择也有讲究。如果半径设得太小,没有完全覆盖噪声能量扩散的区域,残余条纹还会若隐若现;设得太大,会把附近真实频率也削掉,导致图像细节损失。我的经验法则是:先用半径3试,看频谱中亮点周围是否还有残留,如果还有就逐步加大到5、7,直到空间域图像里看不到明显条纹为止。如果你观察到噪声频点呈十字扩散(很多图像的噪声都是横竖方向的),可以考虑用矩形掩码而不是圆形,覆盖效果更准。

陷波滤波之后,图像会有一定程度的轻微模糊或细节损失,这是正常现象。为了尽量保住细节,可以用一个加权方式让掩码边缘平滑过渡,而不是直接0和1硬切换。过渡带越宽,振铃效应越弱,图像越自然。我常用高斯型陷波器替代硬圆掩码,效果会更好。不过绝大多数项目里硬圆掩码已足够用,先把功能跑通再优化不迟。

4.3 性能对比:什么时候该用频域

前面提到卷积定理,这里给一组直观的数字让大家感受性能差异。假设图像是512×512,空域卷积核是15×15。空域卷积需要对每个像素做225次乘加,总计算量大约是512×512×225,约5900万次。如果换成FFT:先fft2一次、频域乘法一次、ifft2一次,总计大约几百万次量级,实际耗时往往优于空域大核卷积。

但当卷积核只有3×3或5×5时,空域卷积只需要9或25次乘加,加上边界处理也很轻量,这时用FFT反而额外引入频域变换的开销,不仅没有提速,还平白增加内存占用。所以我的经验是:卷积核大于15×15时值得考虑频域方案;小于等于7×7时直接用空域;中间地带可以两个都试一下,以实际运行时间为准。如果滤波器需要频繁调整,比如做交互式参数调节,频域方案还有个优势:每次只需要修改H并再做一次逆变换,而空域方案每次都要重新卷积,频域缓存会带来额外收益。

还有一点容易被忽略的是内存。FFT变换要求处理复数数组,内存占用是原始灰度图的好几倍。在一些嵌入式或移动端环境里,512×512的灰度图实际需要约4MB以上的浮点存储,这还不算中间结果。如果内存紧张,可以考虑使用OpenCV的cv2.dft并设置DFT_REAL_OUTPUT或者对图像分块处理,但分块不能随意切了分别变换再拼,因为DFT的全局性不允许把图像切片独立变换再拼接,这会引入严重的边界伪影。需要分块时,常用重叠窗函数和拼接加权平均的方法,这是另一个话题,这里先提个醒。

5. 常见问题与避坑指南

5.1 频谱移位与频率坐标的理解

fftshift不只是为了好看,它直接影响你对频率坐标的解读。未移位时,数组索引(0,0)对应零频,也就是直流分量;索引(k, l)对应的实际频率是k/(MΔx)和l/(NΔy)。做完fftshift后,零频位于中心,索引到实际频率的换算也会变化:中心就是0频率,离中心越远频率越高。理解这个换算关系,你才能正确判断频谱中某个亮点对应的空间结构到底是多少周期每像素。

举个例子:在512×512的图像中,shfit后中心坐标是(256, 256)。如果某个噪声亮点在(256, 204),距离中心52个像素,意味着噪声在垂直方向上的周期大约是512/52≈9.8个像素。这个推算对后续设计滤波器很有用,因为它帮你验证当前处理的频点是否与你在空域目测到的条纹周期一致。

另外要注意,频率坐标的轴方向正负其实不重要,因为幅度谱是共轭对称的,正负频率出现成对亮点。通常我们只需要处理其中一两个亮点,但要记得:如果噪声是垂直条纹(水平方向周期变化),频谱中心水平方向对称出现两个亮点;如果噪声是水平条纹,亮点在垂直方向。不要搞反了。

5.2 幅度谱动态范围:为什么显示出来全黑

几乎每次演示DFT代码,都会有人问“为什么我的频谱图全黑”或者“为什么只有一个白点”。原因在3.2里已经提过:幅度谱的数值范围跨越好几个数量级,而显示设备(或者plt.imshow默认)会把数据线性映射到0~255或0~1的灰度范围。DC分量如果值是一亿,其他分量只有几百,其他分量相对的亮度值就只有0.000几,肉眼当然看不到。

解决的办法就是对幅度取对数:log(1+|F(u,v)|)。取对数之后,一个从1到1e8的区间会被压缩到大约0到18.4的区间,动态范围小了,低幅度的细节才在显示上可见。注意log压缩会影响你对频谱真实强弱对比的感知,它只是可视化手段,不是数据处理手段。如果你要做滤波,操作对象应该还是原始幅度谱或复数谱,不要用log后的频谱去乘滤波器,否则滤波器的衰减关系会被对数扭曲,得到的结果完全不对。

5.3 补零、边界振铃与DC分量

最后整理几个常见问题:

现象原因处理建议
频谱图角落很亮没做fftshift做ifftshift
频谱图全黑动态范围过大用log(1+幅度)显示
滤波后图像出现黑白波纹滤波器过渡带太陡使用平滑过渡带或窗函数
逆变换后虚部很大某一步操作破坏了共轭对称性检查是否只修改了实部或虚部
重建图像整体亮度偏差动了DC分量处理时保留F(0,0)或相应中心点

关于补零,我多说两句。很多人为了让频谱更光滑,在fft2之前把图像用0扩充成2倍大小。补零确实会增加频率采样点数,让频谱的显示更平滑,但它并不会提高频率分辨率,也就是说无法区分原本就靠得很近的两个频率。它的本质是在频域做插值,并没有增加新的信息。补零还有一个代价,就是频谱会出现明显的插值波纹,看起来像加了窗泄漏,如果后续要做精确的频率位置测量,可能会带来误导。所以我认为:除了FFT尺寸有要求(比如某些硬件加速要求2的幂次)或者做卷积需要避免周期卷绕,一般图像分析不必盲目补零。

边界振铃是频域滤波的老朋友。图像在DFT看来是一个周期性重复的无限信号,边界处左边缘和右边缘不连续,相当于在信号里注入了一个阶跃,频谱中就会出现振铃。处理办法有三种:一是滤波前用反射或复制的方式扩边,再裁掉;二是对图像施加一个窗函数,让边界平滑地衰减到0;三是滤波器设计时保留少量高频,不要一刀切到0。我自己的经验是,工业检测场景里用反射扩边配合平滑过渡的滤波器,振铃最小。

DC分量是频谱正中心那个点,它的值等于图像所有像素灰度之和。你把DC分量调大,整张图像变亮;调小,图像变暗。如果你想做的是亮度归一化或者对比度调整,可以直接操作空域像素,没必要绕到频域来;但如果你在频域处理边缘或纹理,一定注意不要无意间改动DC点,否则重建后图像的整体亮度会跟原图对不上。我调试时习惯在处理前把DC分量单独记下来,最后对比一下,排查问题时很有用。

最后分享一个我自己的习惯。现在凡是接到图像处理的活,我第一件事不是急着写滤波代码,而是先对样张做一次DFT频谱分析,把原始图像、幅度谱、相位谱三张图并排打印出来看一眼。这一步的成本只有几十毫秒,但能帮你快速判断问题属于哪种类型:是背景条纹、传感器噪声还是高频干扰。踩过几次坑之后,我发现多数看似复杂的空域困扰,在频谱上一眼就能定位根因。希望这篇笔记能让你在处理图像时也多一个频域的视角,少走一些弯路。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询