从fft2到自动截止频率：二维傅里叶变换与频域滤波MATLAB实战

做图像频域滤波这件事我自己最开始翻车翻得很彻底拿一张照片直接fft2再imshow(abs(F))屏幕上是一片刺眼的白色中间一个小亮点周围什么都看不清。后来才知道问题不在算法而在我对频谱的排布和动态范围完全没有概念。这篇东西就是记录我从只会调用fft2到能自己写出一套二维傅里叶变换滤波器源代码的完整过程包含可直接运行的MATLAB代码覆盖二维傅里叶变换、频域低通/高通/带通滤波、以及如何让滤波器自动确定截止频率这几个核心点适合正在做图像处理课设、需要对二维信号做频域分析、或者想把MATLAB这堆函数真正串起来的读者参考。1. 为什么直接fft2看不出门道先建立频域直觉1.1 图像在频域里到底长什么样图像本质上是二维离散信号f(x,y)傅里叶变换把它拆成一堆不同方向、不同频率的正弦平面波。高频部分对应像素值变化剧烈的位置也就是边缘、纹理、噪声低频部分对应亮度缓慢变化的平坦区域比如天空、背景、大块的色块。但初次接触的人很容易把二维频谱的坐标和图像坐标搞混。fft2返回的矩阵F(u,v)原点在左上角F(1,1)是直流分量也就是整幅图像的平均亮度。只有用fftshift把四个象限对调之后低频才会跑到正中央高频跑到四周。这一步看起来无关紧要实际上一漏掉它后面所有半径、截止频率、中心坐标全都会算错。1.2 动态范围问题是白屏的元凶另一个新手必踩的坑就是频谱直接abs显示全是白的。原因很简单直流分量和低频成分的幅度往往是高频的几十上百倍直接imshow会把显示范围拉伸到最大值低频之外的地方全部压成黑色。解决方法是做对数压缩F fft2(double(img)); Fc fftshift(F); S log(1 abs(Fc)); imshow(S, []);log(1 x)既能压缩大数值又能避免log(0)的问题。这一步做不做直接决定你能否从频谱里看出任何东西。1.3 从频谱能读出什么把频谱可视化之后你会看到中心亮斑是直流分量对应图像平均亮度。经过中心的亮线或十字通常来自图像中的水平/垂直边缘或传感器噪声。偏离中心的对称亮点对应周期性纹理比如布纹、栅栏、印刷网纹。低通滤波的直观目标就是保留中心附近削弱外围高通滤波相反。有了这个直觉后面的滤波器就只是一个怎么画半径的问题。2. 从矩阵到物理频率频谱布局与频率网格的生成2.1 用通用公式构造二维频率坐标所有频域滤波器的第一步都是生成一幅和图像同样大小的网格图每个像素位置记录它到频谱中心的距离。MATLAB里没有现成的fftfreq但手动构造非常简单。假设图像有M行、N列[M, N] size(img); cx floor(M/2) 1; cy floor(N/2) 1; [cols, rows] meshgrid(1:N, 1:M); D sqrt((cols - cy).^2 (rows - cx).^2);这里cols是每个像素的列坐标rows是行坐标。D就是每个像素到中心点(cx, cy)的欧氏距离。所有理想低通、巴特沃斯、高斯滤波器本质上都是关于D的函数。注意如果图像尺寸是偶数中心点其实落在四个像素的交界上简单取floor(M/2)1就可以实际使用误差可以忽略。2.2 具体频率值怎么换算如果你需要把像素半径换算成物理频率比如单位是周期/像素可以这样理解经过fftshift后从中心向右移动一个像素相当于水平方向增加1/N周期/像素的频率。对于采样率fs的连续信号这个频率就是fs/N。大多数图像处理的截止频率讨论都直接以像素半径为准不需要额外换算。2.3 频谱中心化和ifftshift的成对使用fftshift和ifftshift一定要成对出现。滤波流程必须按照正变换→移位→滤波→反移位→逆变换的顺序执行F fft2(double(img)); Fc fftshift(F); Gc H .* Fc; % H是滤波器尺寸和Fc一致 G ifftshift(Gc); g real(ifft2(G));漏掉ifftshift会导致输出图像在四个角出现强烈的块状错位看起来就像图像被切碎重排了。这个错误非常典型也特别容易排查。3. 滤波器设计从理想矩形到平滑过渡的数学模型3.1 理想低通为什么会产生振铃理想低通滤波器是最直观的写法H_ideal double(D D0);它的频率响应是矩形半径D0以内全保留以外全丢弃。但空域里的脉冲响应是sinc函数带有正负振荡高频被硬生生截断后会引发吉布斯现象图像边缘附近会出现一系列明暗交替的纹路这就是振铃。我做实验时第一次用理想低通去除噪声噪点确实没了但人像轮廓周围出现了一圈水波纹看起来比噪声还难受。所以实际项目中理想滤波器主要用于教学演示不建议直接用于图像平滑。3.2 巴特沃斯和高斯滤波器平滑过渡带来的收益巴特沃斯滤波器有一个可以调节的阶数n阶数越高过渡带越窄越接近理想滤波器H_butter 1 ./ (1 (D ./ D0).^(2*n));当n1时过渡带很宽振铃几乎不可见但也会混入一些中频成分当n增大到 4、5锐利度和振铃同时增加。折中建议从n2起步。高斯滤波器没有振铃问题频域是平滑的高斯形状空域也是高斯形状不会出现负值振荡H_gauss exp(-(D.^2) ./ (2 * D0^2));它的代价是过渡带偏宽选择D0偏小时可能会压低没想删除的中频信息。图省事且追求稳妥的时候高斯是我的默认选择。3.3 高通、带通与带阻的统一定义高通滤波器就是1减低通H_hp 1 - H_lp;带通和带阻可以用两个截止半径相减得到% 带通保留 D1 到 D2 之间的频带 H_bp (D D1) (D D2); % 带阻剔除 D1 到 D2 之间的频带 H_bs 1 - H_bp;巴特沃斯和高斯版本只需要把逻辑表达式换成组合运算。高通滤波实际做下来会得到一个中心暗、四周亮的滤波器图和低通的视觉感受完全相反。4. 卷积定理与边界效应为什么频域滤波会卷出黑边4.1 频域相乘等价于循环卷积教科书会告诉你空间域的卷积对应频域相乘。但这段关系更严谨的表述是离散傅里叶变换下的频域相乘对应的是循环卷积不是普通线性卷积。循环卷积意味着图像左右边界、上下边界是首尾相连的左边缘和右边缘会互相影响。所以直接用和图像同尺寸的滤波器做频域乘法边界处的伪影不可避免尤其在使用高通或窄带带通滤波器时图像的四周常常出现深色或亮色条带看起来非常不自然。4.2 先用镜像扩展再滤波是缓解边界伪影最实用的办法如果想减少边界效应比较通用的做法是在滤波前对图像做镜像扩展滤波完成后裁掉扩展部分pad floor(min(M, N) * 0.1); img_pad padarray(img, [pad pad], symmetric); Fp fft2(double(img_pad)); Fcp fftshift(Fp); % 生成和 extend 后尺寸相同的滤波器 Hp Gcp Hp .* Fcp; gp ifftshift(Gcp); g_pad real(ifft2(gp)); g g_pad(pad1:padM, pad1:padN);镜像扩展的英文选项是symmetric它比replicate的边界特性更好。replicate是复制边缘像素会引入一个阶跃symmetric让边界处一阶连续频谱衰减更快伪影更少。不过需要特别说明扩展后生成滤波器时D网格也必须用扩展后的尺寸重新生成否则矩阵尺寸不匹配MATLAB会直接报维度错误。4.3 MATLAB数组复数、溢出与显示范围频域滤波之后拿到的是复数矩阵取实部直接显示是正确的但前提是你确认虚部足够小。多数情况下由于浮点误差虚部数量级在1e-10左右直接real(ifft2(G))没问题。如果发现虚部很大通常说明H的频域共轭对称性被破坏了比如用不太规整的逻辑运算或手动赋值导致不对称。另外imshow默认认为输入是[0,255]范围。如果你的处理结果在[0,1]附近直接用imshow(g)会变成全黑。正确做法是imshow(mat2gray(g));或者显式指定范围imshow(g, [min(g(:)), max(g(:))]);这是我认为图像处理里最容易让人怀疑自己代码出错的地方其实只是显示范围的问题。5. 让滤波器自动起来基于频谱能量的截止频率选择5.1 为什么需要自动确定截止频率很多滤波器教程把D0写死成一个常数比如50、80。但实际面临不同分辨率的图像、不同的噪声强度、不同的应用场景时手调D0极其痛苦。你在这个图像上调好的值换一张图就可能要么残留噪声要么把文字细节抹掉。如果能让滤波器根据频谱本身自动判断多少频率范围该保留体验会好很多。一个非常务实的方法是从低频往高频逐步累积频谱能量当累积能量达到总能量的某个阈值时把这个半径作为自动确定的截止频率D0。5.2 累积能量谱的MATLAB实现[M, N] size(img); F fft2(double(img)); Fc fftshift(F); S abs(Fc).^2; % 功率谱 [cols, rows] meshgrid(1:N, 1:M); cx floor(M/2) 1; cy floor(N/2) 1; D sqrt((cols - cy).^2 (rows - cx).^2); % 把所有像素按半径从小到大排序同时带上对应的功率 [D_sorted, idx] sort(D(:)); S_sorted S(idx); cum_energy cumsum(S_sorted); cum_energy cum_energy / cum_energy(end); % 归一化到 0~1 energy_ratio 0.99; cut_idx find(cum_energy energy_ratio, 1, first); D0_auto D_sorted(cut_idx); fprintf(自动确定的截止频率 D0 %.2f 像素\n, D0_auto);然后直接把D0_auto丢进任意滤波器公式即可。这里最关键的是用sort把半径排序再用累积和找到对应分位点整个计算量极小对512×512的图像也是瞬间完成。5.3 阈值的经验取值与注意事项阈值选多少取决于你的目的应用场景建议阈值说明去除明显噪声保留主体结构0.950.98低频能量占比高0.95通常已经能去掉大部分高频噪声基本无损保真0.990.999适合保留细节但可能残留部分噪声高通分析边缘1-0.95高通截止频率可以从低通自动值的倒推中取得需要注意0.99在合成测试图上往往会给出比较大的半径因为周期条纹的能量集中在少数频点上。面对自然图像时0.95可能更接近人工手调的观感。这个方法是自动化的起点不是终点具体阈值要根据视觉结果微调。6. 完整示例自动低通去噪与高通锐化一条龙6.1 造一张带噪声的测试图像为了不依赖外部图片我直接用正弦条纹加高斯噪声合成一张测试图M 256; N 256; [x, y] meshgrid(1:N, 1:M); clean 128 60*sin(2*pi*8*x/N) 30*cos(2*pi*14*y/M); noisy clean 12 * randn(M, N); figure; subplot(1,2,1); imshow(clean, [0 255]); title(原始条纹); subplot(1,2,2); imshow(noisy, [0 255]); title(加入高斯噪声);这里噪声标准差取12肉眼能看出颗粒感但主体条纹依然清晰比较接近真实拍摄图像加噪后的状态。6.2 自动低通滤波代码F fft2(noisy); Fc fftshift(F); S abs(Fc).^2; [cols, rows] meshgrid(1:N, 1:M); cx floor(M/2) 1; cy floor(N/2) 1; D sqrt((cols - cy).^2 (rows - cx).^2); [D_sorted, idx] sort(D(:)); cum_energy cumsum(S(idx)) / sum(S(:)); cut_idx find(cum_energy 0.95, 1, first); D0 D_sorted(cut_idx); H_lp exp(-(D.^2) ./ (2 * D0^2)); % 高斯低通 G ifftshift(H_lp .* Fc); denoised real(ifft2(G)); figure; imshow(denoised, [min(noisy(:)), max(noisy(:))]); title(自动低通去噪结果);这里我把自动阈值设为0.95D0会偏小一些噪声压制更明显。如果你发现条纹边缘有点糊把阈值调到0.98或0.99再试。6.3 高通锐化和边缘提取高通滤波器同样可以通过自动低通反推H_hp 1 - H_lp; G_hp ifftshift(H_hp .* Fc); hp_result real(ifft2(G_hp)); figure; imshow(mat2gray(hp_result)); title(高通分量);高通结果通常会出现负值直接imshow会得到一个很奇怪的黑白分布。用mat2gray之后正负值被拉伸到0到1边缘信息会以灰底亮线的方式呈现。如果不想要灰底也可以取绝对值再归一化。这个步骤看你要做边缘检测还是锐化锐化通常是把高通分量按比例加回原图sharp noisy 0.8 * hp_result; imshow(sharp, [min(noisy(:)), max(noisy(:))]);系数0.8是我常用的起点调得太大边缘会出现白边。6.4 用PSNR和SSIM做量化对比只靠眼睛判断不够客观MATLAB里可以算峰值信噪比psnr_noisy psnr(uint8(noisy), uint8(clean)); psnr_denoised psnr(uint8(denoised), uint8(clean)); fprintf(加噪 PSNR %.2f dB, 滤波后 PSNR %.2f dB\n, psnr_noisy, psnr_denoised);一般来说高斯白噪声在12个灰度级标准差时PSNR大概在26~28 dB自动低通后能提到32 dB以上就算效果不错。我实测下来0.95阈值对这张条纹图能把PSNR从26.5提升到33.1左右条纹边缘虽然没有原图锐利但视觉上已经非常接近。这个完整流程就是标题里说的二维傅里叶变换滤波器自动处理的最小可运行版本。7. 排错记录我在这个代码上踩过的几个坑7.1 忘了fftshift导致滤波器把图像切碎我的第一个坑是把fftshift忘在了脑后直接在F fft2(img)的基础上生成滤波器。结果滤波器的低频区域在矩阵左上角高频在右下角相乘后输出看起来像被平均分成了四块每块内容都不完整。教训如果你用D距离公式设计滤波器必须先fftshift频谱再用ifftshift把滤波结果恢复。也可以用不经过fftshift的方式设计滤波器但需要把滤波器四个象限重新排列麻烦得多。统一采用移位→滤波→反移位的流程最不容易出错。7.2 滤波器尺寸和图像尺寸不一致导致维度错误业余玩家很容易在扩展图像后用旧尺寸生成滤波器img_pad padarray(img, [pad pad], symmetric); [M2, N2] size(img_pad); % 这里必须用 M2 N2 重新生成 cols/rows/D一次两次维度错误MATLAB会直接提示Matrix dimensions must agree。这个是好事比静默出错强。但如果用了旧版H新版本MATLAB有时候会隐式扩展或广播反而掩盖了问题结果图像就变成混叠状态。7.3 逆变换后虚部偏大输出有重影有一次我为了做带阻滤波器直接对H做了一些非对称处理比如H 1 - ((D D1) (D D2))结果逆变换后虚部居然和实部一个量级。一开始我没取实部直接imshow(g)画面出现了明显的重影。后来意识到H破坏了共轭对称性。解决办法检查H是否关于中心点对称。标准的低通、高通、带通、带阻按上述距离公式生成都是对称的不要在某一边手动修改单个区域否则就会破坏对称性。7.4 显示全黑不等于处理失败这个坑最容易让人心态崩溃。滤波完图像明明有内容但窗口里全黑。原因多半是输出数据范围不是[0, 255]imshow(denoised, []); % 最简单 imshow(mat2gray(denoised));我喜欢imshow(denoised, [])这种写法让 MATLAB 自动找最小值和最大值映射到显示范围。用在分析阶段很方便但如果你想保存结果需要用im2uint8(mat2gray(denoised))把数据转回uint8。7.5 边界伪影被误认为滤波效果差用高通滤波去噪用稍大一点的D0图像四边会出现明显的亮框。很多新手以为是滤波器设计错误其实这是循环卷积的边界效应。先用镜像扩展滤波后再裁剪回来能显著缓解。我在实际项目里会把边界伪影是否明显作为选择滤波器的重要参考。如果一副图像本身就带黑边比如摄影作品边界伪影会被期刊或显示器放大那就必须做扩展处理。7.6 调参时先看频谱再决定滤波策略最后一个经验也是我觉得最值钱的一条遇到任何图像问题先把频谱画出来。如果频谱中心外缘有一圈均匀分布的高频能量那就是高斯噪声如果频谱里出现一对集中的亮斑那就是周期性纹理噪声用带阻滤波比低通滤波有效得多如果频谱能量沿某一方向拉长那说明图像有方向性纹理这时候应该用扇形或方向性滤波器而不是圆形滤波器。仅靠高斯低通高通这两个按钮解决不了所有二维信号问题。先把频谱看成图像诊断报告再决定用哪种滤波器效率会高很多。我保存了这套源代码之后最常用的部分反而是 5.2 节那段累积能量计算后来遇到不同分辨率、不同内容的图像我都是先跑一遍它拿到D0_auto再决定手写滤波器的参数。对刚开始接触二维傅里叶变换的朋友我的建议是不要只在MATLAB里用封装好的imgaussfilt一定亲手写一遍fft2Hifft2的完整链路跑通之后你对频域滤波的理解会有一次实质上的提升。