傅里叶滤波去噪 - 频率域图像处理
傅里叶变换与图像的频率表示
傅里叶变换将图像从空间域转换到频率域,揭示图像中不同频率成分的分布。这为基于频率特性的图像处理提供了强大工具。
基本概念:
- 低频:图像中缓慢变化的区域(大面积均匀色块、渐变)。代表图像的整体结构和亮度分布
- 高频:图像中快速变化的区域(边缘、纹理、噪声)。代表细节和突变
- 频谱图:2D FFT 的幅度谱,中心为低频(DC 分量),向外频率递增
2D DFT(离散傅里叶变换):
- 将 MxN 的图像转换为 MxN 的复数矩阵
- 每个复数的幅度表示该频率的强度,相位表示位置信息
- 使用 FFT(快速傅里叶变换)算法高效计算,复杂度 O(N log N)
2D DFT 的定义式与计算量:傅里叶变换的英文为 Fourier Transform。MxN 图像 f(x,y) 的 2D DFT 定义为 F(u,v) = ΣΣ f(x,y) × exp(-j2π(ux/M + vy/N))。F(u,v) 是频率 (u,v) 处的复数值,具有幅度 |F(u,v)| 与相位 arg(F(u,v))。直接计算 DFT 为 O(N⁴),而 FFT (Fast Fourier Transform) 算法把它降到 O(N² log N):1,024 × 1,024 图像的 FFT 在现代 CPU 上只需数毫秒即可完成。
低通滤波器 - 去噪基础
低通滤波器保留低频成分(图像结构)、抑制高频成分(噪声和细节)。是最基本的频率域去噪方法。
常见低通滤波器:
- 理想低通:截止频率内完全通过,外部完全阻断。会产生振铃效应(Gibbs 现象)
- 巴特沃斯低通:平滑过渡,n 阶巴特沃斯在截止频率处衰减 3dB。阶数越高越接近理想滤波器
- 高斯低通:以高斯函数为传递函数,无振铃效应。等价于空间域的高斯模糊
频率域滤波流程:
- 对图像进行 FFT 得到频谱
- 将频谱与滤波器传递函数逐元素相乘
- 对结果进行逆 FFT 得到滤波后的图像
截止频率选择:太低会过度模糊(丢失细节),太高则去噪不充分。通常通过实验或基于噪声频率特性确定。
两种低通滤波器的传递函数:低通滤波器的英文为 Low-Pass Filter。Butterworth 低通滤波器定义为 H(u,v) = 1 / (1 + (D(u,v)/D0)^(2n)),具有平滑的截止特性;高斯低通滤波器定义为 H(u,v) = exp(-D(u,v)² / (2D0²)),与空间域的高斯模糊 (cv2.GaussianBlur) 等价,两者存在 σ = M/(2πD0) 的关系。它的去噪品质与 Butterworth 相当,优点是完全无振铃 (ringing-free)。
截止频率 D0 的定量确定:常用做法是取包含功率谱 90-95% 能量的频率作为 D0,或依据噪声的频谱特性来决定 D0。
高通和带通滤波器
高通滤波器保留高频(边缘、细节)、抑制低频(背景)。带通滤波器保留特定频率范围。
高通滤波器:
- 传递函数 = 1 - 低通传递函数
- 效果:提取边缘和细节,去除均匀背景
- 应用:边缘增强、锐化(原图 + 高通结果 = 锐化图像)
带通滤波器:
- 仅保留特定频率范围内的成分
- 实现:高通(低截止)x 低通(高截止)
- 应用:提取特定尺度的纹理,分析周期性结构
带阻(陷波)滤波器:
- 抑制特定频率范围,保留其余
- 用于去除已知频率的周期性噪声(如扫描线、摩尔纹)
由低通导出的各类滤波器:高通滤波器定义为 H_HP(u,v) = 1 - H_LP(u,v) (H_LP 为低通滤波器)。用于锐化时的合成形式为 enhanced = original + k × highpass (k 为锐化强度)。带通滤波器定义为 H_BP(u,v) = H_LP(D_H) - H_LP(D_L),只通过 D_L 到 D_H 的频带;带阻滤波器定义为 H_BR(u,v) = 1 - H_BP(u,v),用于去除特定频带。
陷波滤波器 - 周期性噪声去除
陷波滤波器(Notch Filter)针对频谱中特定位置的噪声峰值进行抑制,是去除周期性噪声的最有效方法。
周期性噪声的来源:
- 电磁干扰(如 50/60Hz 电源干扰)
- 扫描设备的机械振动
- 摩尔纹(两个周期性图案的干涉)
- 传感器的固定模式噪声
陷波滤波流程:
- 对图像进行 FFT,观察频谱
- 识别噪声对应的频率峰值(通常表现为频谱中的亮点)
- 在这些位置放置陷波(将对应频率的幅度设为 0 或大幅衰减)
- 逆 FFT 恢复去噪图像
陷波形状:可以是圆形(抑制特定频率)、带状(抑制特定方向的频率)或自定义形状。使用高斯衰减而非硬截断可减少伪影。
高斯型陷波的表达式:高斯型陷波 H(u,v) = 1 - exp(-D1²/(2σ²)) × exp(-D2²/(2σ²)) 可实现平滑的去除 (D1、D2 为到一对陷波中心的距离)。
维纳滤波与逆滤波
维纳滤波在频率域中同时处理去噪和去模糊,是最优线性滤波器(最小均方误差意义下)。
逆滤波:
- 直接除以退化函数(PSF 的傅里叶变换)恢复原图
- 问题:在 PSF 频谱接近零的频率处,噪声被极度放大
- 实际中几乎不可用,除非噪声极低
维纳滤波:
- 在逆滤波基础上加入噪声功率谱的正则化
- 在 PSF 频谱小的频率处自动降低恢复增益,避免噪声放大
- 需要估计信噪比(SNR)或噪声功率谱
- 是去模糊和去噪的最佳线性折中
参数估计:噪声功率谱可从图像的平坦区域估计。信号功率谱可用自然图像的统计模型(1/f 衰减)近似。实践中常用单一参数 K(噪声信号功率比)简化。
逆滤波的破绽与维纳滤波的定义:维纳滤波的英文为 Wiener Filter。退化图像 G(u,v) = H(u,v)F(u,v) + N(u,v) 复原原图 F 的最简做法是逆滤波 F_hat = G/H;但在 H(u,v) 接近 0 的频率处噪声项 N/H 被放大,结果随即崩坏。维纳滤波定义为 W(u,v) = H*(u,v) / (|H(u,v)|² + S_n(u,v)/S_f(u,v)),其中 H* 是 H 的复共轭,S_n 为噪声功率谱、S_f 为信号功率谱。
正则化参数与运动模糊:实际中噪声与信号的功率谱往往未知,因此把 S_n/S_f 用常数 K 近似。相机手持抖动造成的运动模糊,在频率域表现为特定方向上的 sinc 函数。
实现指南 - Python 中的傅里叶滤波
使用 NumPy 和 OpenCV 在 Python 中实现频率域图像滤波。
基本流程:
import numpy as np, cv2- FFT:
f = np.fft.fft2(img); fshift = np.fft.fftshift(f) - 创建滤波器(与图像同尺寸的掩码)
- 应用滤波:
filtered = fshift * mask - 逆 FFT:
result = np.abs(np.fft.ifft2(np.fft.ifftshift(filtered)))
高斯低通滤波器创建:
- 创建与图像同尺寸的网格坐标
- 计算每个点到中心的距离
- 应用高斯函数:
mask = np.exp(-(dist**2) / (2 * sigma**2))
注意事项:
- FFT 前对图像进行零填充(扩展到 2 的幂次)可加速计算并避免循环卷积
fftshift将零频移到中心,便于设计对称滤波器- 滤波后取实部或幅度,丢弃微小的虚部(数值误差)
- OpenCV 的
cv2.dft比 NumPy 的np.fft.fft2在大图像上更快
NumPy 的最小实现:FFT 用 F = np.fft.fft2(image) 计算,再用 F_shift = np.fft.fftshift(F) 把低频移到中心。滤波后用 F_ishift = np.fft.ifftshift(F_filtered) 还原,再用 result = np.abs(np.fft.ifft2(F_ishift)) 逆变换回空间域。滤波器可这样生成并应用:rows, cols = image.shape; crow, ccol = rows//2, cols//2; D = np.sqrt((np.arange(rows)[:,None]-crow)**2 + (np.arange(cols)[None,:]-ccol)**2); H = np.exp(-D**2/(2*D0**2)),然后 F_filtered = F_shift * H。
频谱可视化与加速:用 magnitude = 20*np.log(np.abs(F_shift)+1) 可视化频谱并定位峰值坐标 (u0, v0)。用 optimal_size = cv2.getOptimalDFTSize(n) 可取得最优尺寸;实数图像的 FFT 还可利用对称性用 np.fft.rfft2() 把计算量减半。GPU 加速可用 CuPy 的 cupyx.scipy.fft.fft2(),相对 CPU 可期待 10-50 倍提速。