简介以MATLAB为工具的傅里叶变换与4F系统仿真资源面向图像处理入门者和光学信息处理方向学生解决频域分析、滤波操作及光学系统模拟中的动手实践问题。压缩包共3个文件包含MATLAB脚本、示例图像和FIG图窗文件整体82.54MB脚本覆盖从读取图像、fft2正变换、频域掩模处理到ifft2逆变换的完整流程示例图片供测试效果FIG文件可直接查看滤波结果。该资源已有751人学习/下载。运行后可掌握二维傅里叶变换的核心函数调用、滤波器设计与频谱可视化方法并借助4F系统模拟理解光学图像处理的频域机制参照代码修改掩模即可实现高通、低通等不同滤波实验适合课程实验、期末设计或自学进阶。1. 4F 系统仿真为什么绕不开傅里叶变换4F 系统是光学信息处理里最经典的级联结构两个焦距相等的透镜相距 2f物面放在第一个透镜的前焦面像面落在第二个透镜的后焦面中间那个公共焦面就是频谱面。物面到像面的总距离正好是 4 个焦距所以叫 4F。它之所以重要是因为两个透镜在这里分别承担“傅里叶变换”和“逆傅里叶变换”的角色——你在频谱面上加一个挡片或相位板就等价于对图像的二维频谱做了一次滤波。真实光路上调 4F 系统很麻烦透镜要对准、焦面位置要用剪切板反复确认、频谱面上的针孔孔径稍微偏一点像面就完全变形。而在 MATLAB 里这件事可以用 fft2 和 ifft2 两句核心调用复现。你只需要把物光场铺成离散网格做一次正变换进频谱域再在频谱面乘一个滤波 mask最后逆变换回像面就能看到实验里基本一致的滤波结果。这篇文章就沿着这条路走一遍先解决网格与采样再讲透镜怎么“变成”傅里叶变换器然后给出完整的 4F 滤波仿真代码最后用角谱法做交叉验证。2. 用 MATLAB 搭 4F 仿真前的光场网格采样间隔与 fftshiftFFT 本身不知道你算的是光学还是信号它只处理一个复数矩阵。光学仿真里每个矩阵元素代表一个微小面元上的复振幅模长是幅度幅角是相位单位是 V/m 或直接归一化都行。要让 fft2 算出来的结果能对应到 4F 系统里真实的频谱面第一步是你得把坐标网格造对。2.1 光场网格与物理参数我一般先用一组固定参数搭骨架后面所有代码都复用这套变量。波长为 632.8 nmHe-Ne 激光焦距 f0.5 m网格取 1024×1024采样间隔 8 μm。网格总宽 8.192 mm。lambda 632.8e-9; % 波长 (m) f 0.5; % 透镜焦距 (m) N 1024; % 一维采样点数 pitch 8e-6; % 物面采样间隔 (m) L N * pitch; % 物面总宽度 (m) x (-N/2 : N/2-1) * pitch; [X, Y] meshgrid(x); k 2*pi / lambda; % 波数坐标从-L/2排到L/2之前原点在矩阵中心这符合光学里“光轴沿 z、横截面以光轴为原点”的约定。meshgrid生成二维坐标网格后面生成物函数、透镜相位都要在X,Y上求值。注意meshgrid的输出形状第一个输出X的每一行相同第二个输出Y的每一列相同。如果后面把X当横坐标、Y当纵坐标就不会搞反。采样间隔 8 μm 对可见光来说是偏粗的但仿真物是微米到毫米尺度的特征这个尺度足够而且减小 N 能明显加快 FFT。2.2 fftshift 的必要性直流分量不在中心MATLAB 的fft2输出遵循 DFT 的约定零频直流分量在矩阵的四个角上而不是在中心。光学里频谱面关注的是零频在中心、高频向外扩散的分布所以每一轮fft2之后都要做一次fftshiftUf fft2(U0); % 直接变换直流在四角 Uf_shifted fftshift(Uf); % 移中心直流在正中对应地逆变换之前要先ifftshift再ifft2。一个容易出错的地方是fftshift和ifftshift在 N 为奇数时行为不同前者把原索引 0 移到中心后者把中心元素移回索引 0二者互为逆操作但实现细节不同。用 1024 这类偶数点数时两者等价所以我建议网格尺寸固定取 2 的幂能少踩一个坑。2.3 空间频率坐标和物理频谱面坐标fft2输出的每个点对应一个空间频率(fx, fy)单位是 cycles/m周/米。频率坐标必须和网格匹配否则滤波 mask 不知道自己的带宽在哪里。fx (-N/2 : N/2-1) / L; % 频率分辨率 1/L [FX, FY] meshgrid(fx);fftshift之后的Uf_shifted矩阵中中心点对应(0,0)向右一格频率增加1/L。最高频率是1/(2*pitch)即 62500 cycles/m这就是该网格能表示的频谱范围边界。把这些参数的关系整理成一个表调参时对着查就行。物理量符号MATLAB 表达式本示例数值物面采样间隔dxpitch8 μm物面宽度LN * pitch8.192 mm可用最高频率f_max1 / (2*pitch)62500 cycles/m频率分辨率df1 / L122.07 cycles/m频谱面采样间隔dxilambda * f / L3.86e-5 m频谱面物理坐标单位 m由空间频率换算得到xi FX * lambda * f; eta FY * lambda * f;这里FX是空间频率乘以lambda*f得到的是 4F 系统里真实频谱面上的位置。后面做滤波 mask 时mask 半径若按物理尺寸写就必须用这套坐标。2.4 物函数从哪来合成图案优先4F 仿真最适合的输入是二元几何图案和灰度图像两者我都常用。合成图案的好处是不依赖外部文件任何机器上都能跑% 两个圆孔中心分别在 x ±1.5 mm r 0.5e-3; U0 double(sqrt((X - 1.5e-3).^2 Y.^2) r) ... double(sqrt((X 1.5e-3).^2 Y.^2) r);如果物函数来自实验测量数据比如 CCD 采集的振幅分布存成了 CSV用readmatrix读进来后确保它是一个二维矩阵再插值或裁剪到N×N转成 double 就能进入同一条流水线。灰度图也同理imread后取单通道并转 double 即可。这类输入的重点是确保矩阵尺寸与N一致否则频谱坐标会错位。3. 4F 系统的两次傅里叶变换透镜相位、频谱面与像面FFT 只是数值工具真正让 4F 系统成立的是透镜的相位变换作用。理解了透镜相位你才知道为什么频谱面能放滤波器也知道纯fft2仿真里丢掉了哪个全局相位。3.1 透镜是把光场变换到频域的器件薄透镜近似的复振幅透过率是phiLens exp(-1j * k / (2*f) * (X.^2 Y.^2));它的作用是对入射光场加一个与半径平方成正比的相位延迟。入射场U_in经过透镜后变成U_out U_in .* phiLens;这个二次相位因子不是随便写的。入射光场带上球面相位后继续向前传播在透镜后焦面上菲涅尔衍射积分里的二次相位项恰好被透镜相位抵消剩下的积分式正好是入射光场的二维傅里叶变换U_f(fx, fy) ∝ ∫∫ U_in(x, y) * exp(-j*2*pi*(fx*x fy*y)) dx dy其中fx xi/(lambda*f)fy eta/(lambda*f)。这就是“透镜做傅里叶变换”的来历。焦距相同的两个透镜级联第一个把物场变到频域第二个把频域场变回空域于是构成 4F 系统。3.2 频谱面上的球面相位与 4F 的相位补偿严格推导时透镜后焦面上的场不是纯傅里叶变换而是要乘一个球面相位exp(j*k/(2f)*(xi^2 eta^2))。这个相位的存在意味着频谱面上的复振幅带有弯曲波前直接再放一个透镜做逆变换时前面累积的相位和第二级透镜的相位会再次相消最后像面得到的是近似几何像倒像。4F 的精妙之处就在于两段对称的结构自动完成了二次相位的补偿。离散仿真里这个相位怎么处理我见过不少写法是直接fft2拿频谱、ifft2还原完全不加透镜相位这在“只看强度”的场景下完全够用因为频谱面上的球面相位只影响相位分布不影响abs(Uf).^2的功率谱。但如果你要在频谱面放相位型滤波器涡旋相位板、相位二分板等或者关心复振幅的干涉结果就应该把这一项补上Uf_phys fftshift(fft2(U0)) .* exp(1j * k/(2*f) * (xi.^2 eta.^2));xi,eta就是 2.3 节里换算出的物理坐标。这里乘上的相位模拟的是“真实频谱面上的波前弯曲”最终ifft2还原时会自然消掉因此像面强度结果基本不变但中间过程的相位结构更接近实验。写滤波器代码时我通常把这两条路径都保留显示时用Uf_phys滤波时直接用fftshift(fft2(U0))避免多一次复数乘法。3.3 频谱面的可视化频谱动态范围往往跨越好几个数量级线性显示只能看到中心亮点所以可视化用对数强度S log10(1 abs(Uf_phys).^2); imagesc(xi, eta, S); axis image; colormap(gray);xi和eta作为imagesc的坐标轴传入x 轴、y 轴的单位就是米。这一步建议每次滤波前后都看一眼它能确认 mask 中心是否与零频对齐也能直观看到截断半径对应的频谱范围。4. MATLAB 里的 4F 滤波链低通去噪、高通边缘与像面重建这一章给出完整可运行的 4F 滤波仿真。物函数用双圆孔目的是观察衍射条纹和滤波后的细节变化。实际替换成灰度图同样成立。4.1 完整仿真代码close all; clear; % ---------- 物理参数 ---------- lambda 632.8e-9; f 0.5; N 1024; pitch 8e-6; L N * pitch; x (-N/2 : N/2-1) * pitch; [X, Y] meshgrid(x); k 2*pi / lambda; % ---------- 物函数两个圆孔 ---------- r 0.5e-3; U0 double(sqrt((X - 1.5e-3).^2 Y.^2) r) ... double(sqrt((X 1.5e-3).^2 Y.^2) r); % ---------- 频谱面 ---------- Uf fftshift(fft2(U0)); fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx); % ---------- 低通滤波 ---------- rho_c 3162; % 截止频率单位 cycles/m Hlp double(sqrt(FX.^2 FY.^2) rho_c); U_lp ifft2(ifftshift(Uf .* Hlp)); % ---------- 高通滤波 ---------- Hhp 1 - Hlp; U_hp ifft2(ifftshift(Uf .* Hhp)); % ---------- 显示 ---------- figure(Color,w); subplot(1,3,1); imagesc(x*1e3, x*1e3, abs(U0)); axis image; colormap(gray); title(物面); subplot(1,3,2); imagesc(fx, fx, log10(1 abs(Uf).^2)); axis image; title(频谱对数强度); subplot(1,3,3); imagesc(x*1e3, x*1e3, abs(U_lp).^2); axis image; title(低通像面);代码的核心逻辑是fftshift(fft2(U0))把频谱中心移到坐标原点Hlp按空间频率半径做二值 mask乘到频谱上相当于把高于截止频率的分量清零ifftshift把频谱还原成 FFT 的原始排列最后ifft2回到空域。注意abs(U_lp).^2是强度实验里探测器只能看到这个量。4.2 截止频率怎么换算rho_c的单位是 cycles/m它和真实频谱面上针孔半径的换算关系是rho_physical lambda * f * rho_c比如rho_c 3162cycles/m 时对应物理半径632.8e-9 * 0.5 * 3162 ≈ 1e-3m也就是 1 mm 的针孔。这个换算在做实验预演时特别重要你在代码里设的是 1 mm光路里就应该准备直径 2 mm 的针孔。反过来你知道实验上只想放一个 0.5 mm 的小孔代码里的截止频率就是rho_c 0.5e-3 / (lambda * f);算出来约 1581 cycles/m。常见截止半径对应的效果参考下表。物理针孔半径 (mm)截止频率 (cycles/m)像面表现0.51581几乎只剩低频背景双孔细节消失1.03162边缘开始模糊中心结构可辨2.06324接近原图边缘锐利注意频谱面上能表示的最大物理半宽是0.5 * N * lambda*f / L本例约为 19.3 mm所以 2 mm 仍远小于频谱面尺寸属于窄带滤波。rho_c最大不能超过最大频率1/(2*pitch)即 62500否则 mask 失去滤波意义。4.3 高通滤波与边缘提取高通滤波把零频附近的低频分量去掉只保留高频边缘信息效果等价于图像处理里的边缘提取。代码如下U_hp ifft2(ifftshift(Uf .* (1 - Hlp))); imagesc(x*1e3, x*1e3, abs(U_hp).^2); axis image;高通输出里圆孔内部变暗而圆孔的边界会变成一条亮环。这是因为边界对应频谱中的高频分量被原样保留甚至增强。实验上高通滤波常用一个小黑屏挡住频谱中心而不是用大孔径光阑正是对应这里的1 - Hlpmask。一个值得说的是高通滤波去除直流后像面强度的平均值为零所以abs(U_hp).^2里会出现大量暗区这是正常现象不是数值误差。如果你想让边界环更亮可以把 mask 改成带增益的形式Hhp 1 a * Hlp通过调节a控制边缘增强强度。4.4 把物函数换成灰度图换输入只需替换U0那一句cam imread(cameraman.tif); % 灰度图256x256 cam imresize(cam, [N, N]); % 重采样到网格 U0 double(cam) / 255; % 归一化到 [0,1]imresize要凑到N×N否则频率坐标全部错位。摄像师图本身的频谱能量集中在中低频rho_c3162时像面会明显变平滑细节如支架、人脸轮廓被抹掉rho_c6324时基本恢复原图。这一步能直观看到低通在图像上的效果。由于仿真里没有加入噪声低通的去噪效果看不出来。想模拟实验场景可以在物面加高斯噪声U0 double(imread(...)) / 255 0.05 * randn(N, N);噪声的能量分布在整个频谱面上低通 mask 切掉高频后像面的颗粒感明显减少这就更进一步接近真实实验了。5. 用角谱传播交叉验证 4F 仿真三个必查项纯fft2的两步变换算得再顺也只能说明数值自洽不能证明它等价于真实 4F 光路。最稳的验证方式是用角谱传播把 4F 拆成五段传播到透镜、乘透镜相位、再传播到频谱面、乘第二个透镜相位、最后传播到像面。两种路径的结果如果误差在可忽略范围就说明简化模型没有用错。5.1 角谱交叉验证代码% 空间频率网格沿用 4.1 节 fx (-N/2 : N/2-1) / L; [FX, FY] meshgrid(fx); % 自由空间传播 z 米的角谱算子菲涅尔近似 prop (U, z) ifft2(ifftshift(fftshift(fft2(U)) .* ... exp(-1j * pi * lambda * z * (FX.^2 FY.^2)))); phiLens exp(-1j * k / (2*f) * (X.^2 Y.^2)); Uv prop(U0, f); % 物面到透镜1 Uv Uv .* phiLens; % 过透镜1 Uv prop(Uv, f); % 透镜1到频谱面 Uv Uv .* phiLens; % 过透镜2 Uv prop(Uv, f); % 透镜2到像面 err norm(abs(Uv).^2 - abs(U_lp).^2) / norm(abs(U_lp).^2);prop里先fftshift把频谱中心化乘角谱转移函数再ifftshift还原后逆变换。这里值要强调的是角谱路径里每段都是真实的物理传播而简化路径只有两次 FFT所以误差主要来自简化路径省略的球面相位项。对强度而言err通常在 1e-6 量级基本可以忽略。5.2 三个必查项第一个必查项是像面方向。真实 4F 系统成倒像而代码里fft2之后再ifft2得到的是正立像因为两次变换的方向约定抵消了。所以判断仿真是否正确不能看方向要看结构是否与原物一致方向问题属于坐标约定留到和实验对照时再处理。第二个必查项是全通滤波恢复。把 mask 设成全 1像面强度应该与原物强度完全一致误差同样是 1e-6 量级。如果你发现像面有明显条纹或光晕先检查fftshift/ifftshift是否成对使用这是最常见的错误源。第三个必查项是频谱 center 对齐。滤波 mask 的圆心必须对准频谱中心(N/21, N/21)有偏移时像面会出现明显的整体倾斜调制。精确做法是把 mask 构造在频域网格上而不是用imcrop之类的图像工具手动画圆。最后给一个实用技巧如果你要在频谱面模拟带方向性的滤波器比如让图案只在竖直方向模糊可以使用椭圆 maskH_ell double((FX.^2 / a^2 FY.^2 / b^2) 1);其中a,b分别是两个主轴方向的截止频率。这样一套代码就能覆盖 4F 系统的大部分频域滤波实验从各向同性低通到方向滤波参数只需要改一个 mask 表达式。本文还有配套的精品资源点击获取