尧图网站开发 
简介面向图像处理与光学信息研究者的MATLAB傅里叶变换与4F系统仿真资源适合学习数字图像频域分析、4F光学滤波实验模拟的本科生、研究生及工程师。资源共3个文件包含可直接运行的m源代码脚本、用于测试的示例图像以及fig格式结果图压缩包整体82.54MB能帮助快速复现从图像读取、fft2二维傅里叶变换、频域掩模滤波到ifft2反变换重建的完整流程。m脚本提供分步实现代码配合测试图像与结果图使用者可直观对比原始图像、频谱及滤波后图像的差异深入理解空间域与频域转换的物理意义。已有751人学习下载尤其适合作为光学信息处理、数字图像处理课程的课外实践参考也可为相关课题研究提供基础仿真框架有效缩短从理论到代码实现的探索时间。1. 4F系统与傅里叶变换为什么仿真里频谱位置总差半格调过光学平台的人都知道4F系统的对焦是一件磨人的事两个焦距相同的透镜相距 2f物面、频谱面、像面都要卡在焦平面上差一毫米频谱就糊了。而在 MATLAB 里做 4F 仿真问题却反了过来——fouriertransform 的数学公式清清楚楚但 fft2 出来的频谱图要么 DC 分量在角落要么滤波器的半径对不上物理尺寸要么像面左右颠倒。多数教材给的公式是连续的直接照搬到矩阵索引上坐标系和归一化全乱套。这篇文章要做的就是把这层窗户纸捅破把离散傅里叶变换、透镜焦平面的物理坐标和 4F 系统的频谱滤波串成一套可直接改参数的 MATLAB 仿真模板。适合做光信息处理、计算光学成像、全息与光学密码的工程师也适合需要合成“物函数—滤波结果”训练数据的算法工程师。2. 用MATLAB把FFT变成4F系统频谱面坐标映射与最小可跑仿真2.1 透镜后焦面为什么是傅里叶变换从夫琅禾费衍射到fft2单色平面波垂直照射透射率为 t(x, y) 的物体物体紧贴透镜前焦面时后焦面上的复振幅分布由夫琅禾费衍射给出U_f(x_f, y_f) e^(jkf) / (jλf) · ∫∫ t(x, y) · exp[-j2π(x·x_f y·y_f) / (λf)] dx dy积分号内就是二维傅里叶变换的核所以后焦面上的复振幅正比于 t(x,y) 的傅里叶变换。这里 λ 是波长f 是焦距x_f、y_f 是焦面上的物理坐标。这个公式里还挂着 e^(jkf) 和二次相位因子但二次相位因子只影响相位、不影响焦面强度且进入 4F 系统的第二个透镜后会被补偿掉因此在 4F 仿真框架里可以放心省略只保留傅里叶变换关系。把连续积分换成 MATLAB 的 fft2关键不在变换本身而在坐标换标空间频率 u x_f / (λf)v y_f / (λf)。也就是说透镜是把“空间频率轴”原样映射成了“焦面物理坐标轴”缩放系数是 λf。这一步没有做对滤波器设计得再漂亮也落不到正确的位置上。2.2 空间频率、频谱面坐标与数组索引的对应关系离散傅里叶变换的数组索引与空间频率的对应是 4F 仿真里最容易翻车的环节。fft2 输出的 (1,1) 位置对应零频但物理上我们希望零频在矩阵中心所以需要 ifftshift 再 fft2最后再 fftshift 把频谱搬回中心。这个组合里物体坐标的定义必须与 ifftshift 后的数组对齐否则偶数尺寸下会差半个像素。频率间隔由采样参数决定du 1 / (N·dx)其中 N 是单边采样数dx 是物面像素间距。矩阵里相邻元素对应的空间频率差是 du换算到频谱面上的物理间隔就是 λf·du。下表是四个量之间的换算关系物理量符号MATLAB 生成方式物面坐标x(-N/2 : N/2-1) * dx空间频率fx(-N/2 : N/2-1) / (N*dx)频谱面坐标xflambda * f * fx角频率滤波用rhosqrt(fx.^2 fy.^2)注意 fx 已经是 fftshift 坐标系的频率轴这样它与 fftshift(fft2(...)) 的结果天然对齐。如果你的物函数来自实测数据比如从 CSV 读入的波前采样矩阵只需要按同样的 dx 和 N 重新声明坐标网格后续流程完全一致。2.3 一个12行的夫琅禾费衍射最小仿真下面这段代码实现了从矩形孔到焦面频谱的完整过程是后续 4F 仿真的地基% 基本采样参数 N 512; % 采样数取 2 的幂 dx 5e-6; % 物面像素间距单位 m x (-N/2 : N/2-1) * dx; % 物面坐标向量 [X, Y] meshgrid(x, x); % 第一维对应 Y第二维对应 X % 光学参数 lambda 632.8e-9; % He-Ne 激光波长 f 100e-3; % 透镜焦距 100 mm % 物函数边长 1 mm 的矩形孔 t zeros(N, N); t(abs(X) 0.5e-3 abs(Y) 0.5e-3) 1; % 透镜后焦面复振幅 Uf fftshift(fft2(ifftshift(t))); % 频谱面物理坐标 du 1 / (N * dx); % 空间频率间隔 fxv (-N/2 : N/2-1) * du; % 空间频率单位 cycle/m xfv lambda * f * fxv; % 频谱面坐标单位 m % 显示频谱强度 imagesc(xfv*1e3, xfv*1e3, abs(Uf)); axis image; xlabel(x_f (mm)); ylabel(y_f (mm));这段代码里最值得注意的三个点一是 ifftshift(t) 把物体的几何中心搬回数组原点 (1,1)fft2 完成变换后 fftshift 再把零频搬回中心物体坐标 x 必须与这个流程配套定义二是 meshgrid(x, x) 生成的 X 沿矩阵第二维变化、Y 沿第一维变化与 fft2 先沿行后沿列的变换顺序一致三是显示时用频谱面物理坐标 xfv 而不是频率这样与实验测得的焦面光斑尺寸可以直接对比。把矩形孔换成圆孔、光栅或任意实值图像这套坐标体系都不需要改动。3. 4F系统空间滤波MATLAB实战低通、高通、暗场与相衬3.1 把4F系统封装成函数物函数进像强度出4F 系统在仿真里就是一个三步流水线物函数做傅里叶变换得到频谱频谱乘上滤波器再做逆傅里叶变换拿到像。我把这三步封装成一个函数后续所有滤波实验共用同一套代码function [I, xfv, Uf] sim4f(t, dx, lambda, f, H) % 4F 系统相干成像仿真理想透镜忽略孔径与像差 % t : 物面复振幅NxN 矩阵 % dx : 物面像素间距 (m) % lambda : 波长 (m) % f : 透镜焦距 (m) % H : 频谱面滤波器NxN与 t 同尺寸零频在中心 % I : 像面强度分布 % xfv : 频谱面坐标向量 (m)用于标注 N size(t, 1); Uf fftshift(fft2(ifftshift(t))); Uf_filtered Uf .* H; Uimg ifft2(ifftshift(Uf_filtered)); Uimg ifftshift(Uimg); I abs(Uimg).^2; xfv ((-N/2 : N/2-1) / (N * dx)) * lambda * f; end像面坐标与物面坐标相同但 4F 系统成的是倒像ifft2 对应的坐标翻转正好模拟了这一物理事实显示时用 fliplr 就能恢复正常朝向。函数里 Uimg 只取强度是因为 CCD 和 CMOS 记录的都是强度做相干成像仿真时强度才是与实验可比对的量如果后续要级联下一级光学系统则应该输出复振幅 Uimg。滤波器 H 必须与 fftshift 后的 Uf 同坐标系即零频在矩阵中心、频率轴按 fxv 定义。3.2 低通滤波用平滑截止代替硬截断低通滤波对应频谱面上放一个光阑拦掉高频分量。最直觉的实现是理想低通半径内为 1、半径外为 0但它会在像面产生明显的吉布斯振铃。我一般用 tanh 做平滑过渡既保留低通的物理含义又避免振铃污染定性结论N 512; dx 8e-6; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); [fx, fy] meshgrid((-N/2 : N/2-1) / (N*dx), (-N/2 : N/2-1) / (N*dx)); rho sqrt(fx.^2 fy.^2); fc 400; % 截止频率cycle/m w 100; % 过渡带宽cycle/m H_lp 0.5 * (1 - tanh((rho - fc) / w)); % 平滑低通 t im2double(imread(cameraman.tif)); t imresize(t, [N N]); [I, xfv] sim4f(t, dx, 632.8e-9, 100e-3, H_lp); imshow(I);fc 400 cycle/m 意味着物面上周期大于 2.5 mm 的结构能通过更细的结构被抑制。w 控制过渡带宽度w 越小越接近理想低通振铃越明显w 取 fc 的 1/4 到 1/3 是一个经验安全区间。注意这里网格生成用了 meshgrid频率轴 fx 沿第二维扩展与滤波器矩阵 H_lp 的行列角色一致。3.3 高通与边缘增强去掉DC分量后图像还剩什么高通滤波是低通的反向操作H_hp 1 - H_lp。它滤掉零频和低频背景像面上只保留边缘和轮廓视觉上呈现为黑色背景上的亮线。直接看 1 - H_lp 的输出会显得对比度偏低因为边缘以外的低频也被保留了一部分更强硬的边缘增强做法是高频提升滤波器H_edge 1 - exp(-(rho / fc).^2); % fc 为高频通过的尺度 [I_edge, ~] sim4f(t, dx, 632.8e-9, 100e-3, H_edge); imshow(I_edge);exp(-(rho/fc)^2) 在低频处接近 1H_edge 接近 0因此背景被压暗而在高频处 H_edge 趋近 1边缘细节保留。这种滤波器对应光学上的振幅挡板加中心光阑的组合是 4F 边缘增强实验的标准配置。调 fc 时观察像面fc 偏大则只剩最锐利的拐角fc 偏小则边缘有光晕通常从物函数最小特征尺寸的倒数起步调试。3.4 暗场与泽尼克相衬用一块挡板和一片相移片看相位物体纯相位物体比如透明生物切片透射率 t exp(j·φ(x,y))直接成像强度是均匀的什么都看不见。4F 系统处理这类物体的经典招数是在频谱面做手脚。暗场法是在中心放一个不透明小圆盘挡掉零级H_dark ones(N, N); H_dark(rho 60) 0; % 半径 60 cycle/m 的挡板 [I_dark, ~] sim4f(exp(1i*0.5*double(t0.3)), dx, 632.8e-9, 100e-3, H_dark);泽尼克相衬法更进一步把零级光的相位推迟 π/2让相位梯度转化为强度衬度H_pc ones(N, N); idx rho 60; H_pc(idx) exp(1i * pi/2); % 零级相移 [I_pc, ~] sim4f(exp(1i*0.5*double(t0.3)), dx, 632.8e-9, 100e-3, H_pc);暗场输出中相位变化剧烈的边界变亮缓慢变化区域仍是黑的相衬输出则把相位分布近似线性映射为灰度看起来像“浮雕”效果。两类滤波器都不改变振幅靠的是改变频谱各成分之间的相对相位所以滤波器 H 必须是复数矩阵用 exp(1i*pi/2) 而不是乘一个实数系数。下表汇总了本节四种滤波器的传递函数与适用场景滤波器传递函数要点像面效果参数调整方向理想低通rho fc 为 1其余 0平滑振铃明显fc 控制细节保留量平滑低通0.5(1-tanh((rho-fc)/w))平滑振铃弱w 控制过渡带高通 H_hp1 - H_lp边缘亮线配合 fc 抑制低频高频提升1-exp(-(rho/fc)^2)边缘增强、背景保留fc 取特征频率暗场中心圆形挡板置 0相位梯度可见挡板半径大于零级尺寸泽尼克相衬零级区乘 exp(jπ/2)相位映射为强度相移量为 π/2 最佳4. 4F仿真参数边界与常见坑采样率、零填充、振铃与维度顺序4.1 采样间隔、频谱面半径与混叠的定量关系物面像素间距 dx 决定了可表示的最高空间频率 f_max 1/(2dx)也就是离散傅里叶变换的奈奎斯特极限。换算到频谱面这个频率对应一个物理半径rho_max lambda · f / (2 · dx)用 2.3 节的参数 λ 632.8 nm、f 100 mm、dx 5 μm 代入rho_max ≈ 6.3 mm。意思是频谱面上半径 6.3 mm 之外的地方物理上不存在真实频率成分即使滤波器在那里开了通带通过的也只是混叠造成的伪频率。这个式子很有用把 dx 减半频谱面可用的物理半径就翻倍而把 N 翻倍只改变频谱采样的疏密不改变 rho_max。设计滤波器时先算 rho_max再决定 fc 的量级能避免“滤波器写对了但位置在物理上无意义”的尴尬。判断混叠是否严重可以看频谱强度正常频谱在 rho_max 附近强度应衰减到接近本底如果频谱在矩阵边缘仍然很强说明物面采样不足此时要减小 dx而不是靠增大 N 硬撑。4.2 零填充频谱插值还是信息增加零填充是在物面矩阵外围补零把尺寸从 N 扩展到 M。常见认知误区是“补零能提高分辨率”实际上补零只做了一件事让频谱的采样间隔从 1/(N·dx) 变细到 1/(M·dx)频谱曲线被插值得更平滑峰值定位更准但两个靠得很近的真实频率分量并不会因此被分开。补零的另一好处是消除了物体边缘与数组边界之间的周期性卷绕效应因为补零后物体的镜像副本离得更远。实现时注意坐标同步扩展M 变为原 N 的两倍后新坐标向量要用新的 M 重新生成dx 不变否则滤波器与频谱会对不齐。我一般把补零写成 padarray(t, [N/2 N/2], 0)之后所有涉及 x、fx 的向量全部按 M 重建。4.3 滤波器边缘用余弦过渡消除吉布斯振铃硬截断滤波器在截止频率处传递函数从 1 跳变到 0对应空域里出现 sinc 形状的旁瓣这就是振铃的来源。想要保留低通的“截止”语义、又不让振铃掩盖像面细节可以在截止频率附近做一段余弦过渡r rho / fc; % 归一化径向频率 H_cos ones(size(rho)); band r 0.8 r 1; H_cos(band) 0.5 * (1 cos(pi * (r(band) - 0.8) / 0.2)); H_cos(r 1) 0;过渡带从 0.8·fc 开始、到 fc 结束频带宽度固定为 0.2·fc。过渡带越宽振铃越弱但截止特性越不锐利。判断振铃是否可接受看像面中高对比度边缘附近有没有明暗交替的平行条纹有的话就把过渡带放宽到 0.3·fc 再跑一次。这个余弦渐变滤波器与 3.2 节的 tanh 低通本质上是同一种思路实际选哪个取决于你想不想保留一个显式的截止频率参数。4.4 meshgrid与fft2的维度顺序最容易错位的一处fft2 沿第一维行方向变换 y 分量沿第二维列方向变换 x 分量。meshgrid(x, x) 生成的 X 矩阵每行相同即 X 沿列方向变化对应第二维Y 矩阵每列相同即 Y 沿行方向变化对应第一维。所以 meshgrid(x, x) 与 fft2 的维度角色完全一致这是推荐用法。如果换成 ndgrid(x, x)得到的第一维是 x、第二维是 y与 fft2 的角色正好相反物函数和滤波器都会发生转置错位。判别方法很简单画 meshgrid 生成的 X看它沿哪个方向变化再画 ndgrid 生成的 X两者互为转置。同一段代码里物函数、滤波器、频谱三者必须统一用同一种网格生成函数混用 meshgrid 和 ndgrid 是 4F 仿真“像面横竖颠倒”最常见的原因。另外滤波器 H 如果是通过 meshgrid(fxv, fxv) 生成并用 rho 做半径判断的H 的类型应当是 double 或复数 double不能是 logical否则与 Uf 相乘时维度虽然合法、数值却会被隐式转换相位滤波直接失效。5. 用解析解验证4F仿真艾里斑位置与涡旋光生成5.1 艾里斑第一暗环的位置比对4F 仿真的正确性可以用圆孔的夫琅禾费衍射解析解来验证半径为 R 的圆孔频谱面强度是艾里斑第一暗环半径的解析值为 r_Airy 1.22·λ·f / D其中 D 2R 是圆孔直径。仿真里找到强度剖面第一极小值的位置与解析值对比偏差应在 5% 以内R 1e-3; t double(sqrt(X.^2 Y.^2) R); Uf fftshift(fft2(ifftshift(t))); profile abs(Uf(N/21, :)); % 过中心的水平剖面 [~, imin] min(profile(N/21 : end)); % 中心右侧第一极小 r_sim xfv(imin N/2); % 仿真暗环半径 r_theory 1.22 * lambda * f / (2*R);如果 r_sim 比理论值大且剖面边缘翘起多半是物面采样不足、频谱混叠如果 profile 在中心以外没有明显的零点检查 dx 是否小到能够分辨圆孔边缘的陡变。把 R 改成 0.5 mm 和 2 mm 各跑一遍理论值和仿真值同比例变化基本可以确认整个坐标链路没有系统性错误。5.2 光栅级次、能量守恒与涡旋光正弦光栅是另一个廉价而可靠的验证对象t 0.5·(1 cos(2π·x/Λ))频谱面应在 fx ±1/Λ 处出现一对对称的亮峰。用 max 定位峰值所在列换算成频率后与 1/Λ 对比误差小于一个频率采样间隔 du 就说明频率轴标定正确。滤波前后还可以做能量核对sum(abs(Uf_filtered).^2) / sum(abs(Uf).^2) 应等于滤波器透过率平方在频谱面上的占空比低通滤波器改小 fc 时这个比值单调下降任何不单调的情况都指向滤波器矩阵与频谱错位。最后给一个能立刻看出 4F 仿真价值的进阶手法在频谱面放螺旋相位板传递函数 H exp(j·m·θ)θ 是频域极角m 是拓扑荷theta angle(fx 1i*fy); H_vortex exp(1i * 1 .* theta); % m1 的涡旋滤波 [I_vortex, ~] sim4f(t_gaussian, dx, lambda, f, H_vortex);像面会得到中心强度为零的环形光斑环径随 m 增大而增大。这个结果可以直接与轨道角动量光学的理论预期对照也是用 MATLAB 验证“频谱面相位操作如何转化为像面振幅分布”最直观的例子。把 5.1 的艾里斑、光栅级次和这里的涡旋环放在一起4F 仿真的物面采样、频率轴、滤波器坐标系三个环节就都有了独立于代码本身的定量背书。本文还有配套的精品资源点击获取 
简介面向图像处理与光学信息研究者的MATLAB傅里叶变换与4F系统仿真资源适合学习数字图像频域分析、4F光学滤波实验模拟的本科生、研究生及工程师。资源共3个文件包含可直接运行的m源代码脚本、用于测试的示例图像以及fig格式结果图压缩包整体82.54MB能帮助快速复现从图像读取、fft2二维傅里叶变换、频域掩模滤波到ifft2反变换重建的完整流程。m脚本提供分步实现代码配合测试图像与结果图使用者可直观对比原始图像、频谱及滤波后图像的差异深入理解空间域与频域转换的物理意义。已有751人学习下载尤其适合作为光学信息处理、数字图像处理课程的课外实践参考也可为相关课题研究提供基础仿真框架有效缩短从理论到代码实现的探索时间。1. 4F系统与傅里叶变换为什么仿真里频谱位置总差半格调过光学平台的人都知道4F系统的对焦是一件磨人的事两个焦距相同的透镜相距 2f物面、频谱面、像面都要卡在焦平面上差一毫米频谱就糊了。而在 MATLAB 里做 4F 仿真问题却反了过来——fouriertransform 的数学公式清清楚楚但 fft2 出来的频谱图要么 DC 分量在角落要么滤波器的半径对不上物理尺寸要么像面左右颠倒。多数教材给的公式是连续的直接照搬到矩阵索引上坐标系和归一化全乱套。这篇文章要做的就是把这层窗户纸捅破把离散傅里叶变换、透镜焦平面的物理坐标和 4F 系统的频谱滤波串成一套可直接改参数的 MATLAB 仿真模板。适合做光信息处理、计算光学成像、全息与光学密码的工程师也适合需要合成“物函数—滤波结果”训练数据的算法工程师。2. 用MATLAB把FFT变成4F系统频谱面坐标映射与最小可跑仿真2.1 透镜后焦面为什么是傅里叶变换从夫琅禾费衍射到fft2单色平面波垂直照射透射率为 t(x, y) 的物体物体紧贴透镜前焦面时后焦面上的复振幅分布由夫琅禾费衍射给出U_f(x_f, y_f) e^(jkf) / (jλf) · ∫∫ t(x, y) · exp[-j2π(x·x_f y·y_f) / (λf)] dx dy积分号内就是二维傅里叶变换的核所以后焦面上的复振幅正比于 t(x,y) 的傅里叶变换。这里 λ 是波长f 是焦距x_f、y_f 是焦面上的物理坐标。这个公式里还挂着 e^(jkf) 和二次相位因子但二次相位因子只影响相位、不影响焦面强度且进入 4F 系统的第二个透镜后会被补偿掉因此在 4F 仿真框架里可以放心省略只保留傅里叶变换关系。把连续积分换成 MATLAB 的 fft2关键不在变换本身而在坐标换标空间频率 u x_f / (λf)v y_f / (λf)。也就是说透镜是把“空间频率轴”原样映射成了“焦面物理坐标轴”缩放系数是 λf。这一步没有做对滤波器设计得再漂亮也落不到正确的位置上。2.2 空间频率、频谱面坐标与数组索引的对应关系离散傅里叶变换的数组索引与空间频率的对应是 4F 仿真里最容易翻车的环节。fft2 输出的 (1,1) 位置对应零频但物理上我们希望零频在矩阵中心所以需要 ifftshift 再 fft2最后再 fftshift 把频谱搬回中心。这个组合里物体坐标的定义必须与 ifftshift 后的数组对齐否则偶数尺寸下会差半个像素。频率间隔由采样参数决定du 1 / (N·dx)其中 N 是单边采样数dx 是物面像素间距。矩阵里相邻元素对应的空间频率差是 du换算到频谱面上的物理间隔就是 λf·du。下表是四个量之间的换算关系物理量符号MATLAB 生成方式物面坐标x(-N/2 : N/2-1) * dx空间频率fx(-N/2 : N/2-1) / (N*dx)频谱面坐标xflambda * f * fx角频率滤波用rhosqrt(fx.^2 fy.^2)注意 fx 已经是 fftshift 坐标系的频率轴这样它与 fftshift(fft2(...)) 的结果天然对齐。如果你的物函数来自实测数据比如从 CSV 读入的波前采样矩阵只需要按同样的 dx 和 N 重新声明坐标网格后续流程完全一致。2.3 一个12行的夫琅禾费衍射最小仿真下面这段代码实现了从矩形孔到焦面频谱的完整过程是后续 4F 仿真的地基% 基本采样参数 N 512; % 采样数取 2 的幂 dx 5e-6; % 物面像素间距单位 m x (-N/2 : N/2-1) * dx; % 物面坐标向量 [X, Y] meshgrid(x, x); % 第一维对应 Y第二维对应 X % 光学参数 lambda 632.8e-9; % He-Ne 激光波长 f 100e-3; % 透镜焦距 100 mm % 物函数边长 1 mm 的矩形孔 t zeros(N, N); t(abs(X) 0.5e-3 abs(Y) 0.5e-3) 1; % 透镜后焦面复振幅 Uf fftshift(fft2(ifftshift(t))); % 频谱面物理坐标 du 1 / (N * dx); % 空间频率间隔 fxv (-N/2 : N/2-1) * du; % 空间频率单位 cycle/m xfv lambda * f * fxv; % 频谱面坐标单位 m % 显示频谱强度 imagesc(xfv*1e3, xfv*1e3, abs(Uf)); axis image; xlabel(x_f (mm)); ylabel(y_f (mm));这段代码里最值得注意的三个点一是 ifftshift(t) 把物体的几何中心搬回数组原点 (1,1)fft2 完成变换后 fftshift 再把零频搬回中心物体坐标 x 必须与这个流程配套定义二是 meshgrid(x, x) 生成的 X 沿矩阵第二维变化、Y 沿第一维变化与 fft2 先沿行后沿列的变换顺序一致三是显示时用频谱面物理坐标 xfv 而不是频率这样与实验测得的焦面光斑尺寸可以直接对比。把矩形孔换成圆孔、光栅或任意实值图像这套坐标体系都不需要改动。3. 4F系统空间滤波MATLAB实战低通、高通、暗场与相衬3.1 把4F系统封装成函数物函数进像强度出4F 系统在仿真里就是一个三步流水线物函数做傅里叶变换得到频谱频谱乘上滤波器再做逆傅里叶变换拿到像。我把这三步封装成一个函数后续所有滤波实验共用同一套代码function [I, xfv, Uf] sim4f(t, dx, lambda, f, H) % 4F 系统相干成像仿真理想透镜忽略孔径与像差 % t : 物面复振幅NxN 矩阵 % dx : 物面像素间距 (m) % lambda : 波长 (m) % f : 透镜焦距 (m) % H : 频谱面滤波器NxN与 t 同尺寸零频在中心 % I : 像面强度分布 % xfv : 频谱面坐标向量 (m)用于标注 N size(t, 1); Uf fftshift(fft2(ifftshift(t))); Uf_filtered Uf .* H; Uimg ifft2(ifftshift(Uf_filtered)); Uimg ifftshift(Uimg); I abs(Uimg).^2; xfv ((-N/2 : N/2-1) / (N * dx)) * lambda * f; end像面坐标与物面坐标相同但 4F 系统成的是倒像ifft2 对应的坐标翻转正好模拟了这一物理事实显示时用 fliplr 就能恢复正常朝向。函数里 Uimg 只取强度是因为 CCD 和 CMOS 记录的都是强度做相干成像仿真时强度才是与实验可比对的量如果后续要级联下一级光学系统则应该输出复振幅 Uimg。滤波器 H 必须与 fftshift 后的 Uf 同坐标系即零频在矩阵中心、频率轴按 fxv 定义。3.2 低通滤波用平滑截止代替硬截断低通滤波对应频谱面上放一个光阑拦掉高频分量。最直觉的实现是理想低通半径内为 1、半径外为 0但它会在像面产生明显的吉布斯振铃。我一般用 tanh 做平滑过渡既保留低通的物理含义又避免振铃污染定性结论N 512; dx 8e-6; x (-N/2 : N/2-1) * dx; [X, Y] meshgrid(x, x); [fx, fy] meshgrid((-N/2 : N/2-1) / (N*dx), (-N/2 : N/2-1) / (N*dx)); rho sqrt(fx.^2 fy.^2); fc 400; % 截止频率cycle/m w 100; % 过渡带宽cycle/m H_lp 0.5 * (1 - tanh((rho - fc) / w)); % 平滑低通 t im2double(imread(cameraman.tif)); t imresize(t, [N N]); [I, xfv] sim4f(t, dx, 632.8e-9, 100e-3, H_lp); imshow(I);fc 400 cycle/m 意味着物面上周期大于 2.5 mm 的结构能通过更细的结构被抑制。w 控制过渡带宽度w 越小越接近理想低通振铃越明显w 取 fc 的 1/4 到 1/3 是一个经验安全区间。注意这里网格生成用了 meshgrid频率轴 fx 沿第二维扩展与滤波器矩阵 H_lp 的行列角色一致。3.3 高通与边缘增强去掉DC分量后图像还剩什么高通滤波是低通的反向操作H_hp 1 - H_lp。它滤掉零频和低频背景像面上只保留边缘和轮廓视觉上呈现为黑色背景上的亮线。直接看 1 - H_lp 的输出会显得对比度偏低因为边缘以外的低频也被保留了一部分更强硬的边缘增强做法是高频提升滤波器H_edge 1 - exp(-(rho / fc).^2); % fc 为高频通过的尺度 [I_edge, ~] sim4f(t, dx, 632.8e-9, 100e-3, H_edge); imshow(I_edge);exp(-(rho/fc)^2) 在低频处接近 1H_edge 接近 0因此背景被压暗而在高频处 H_edge 趋近 1边缘细节保留。这种滤波器对应光学上的振幅挡板加中心光阑的组合是 4F 边缘增强实验的标准配置。调 fc 时观察像面fc 偏大则只剩最锐利的拐角fc 偏小则边缘有光晕通常从物函数最小特征尺寸的倒数起步调试。3.4 暗场与泽尼克相衬用一块挡板和一片相移片看相位物体纯相位物体比如透明生物切片透射率 t exp(j·φ(x,y))直接成像强度是均匀的什么都看不见。4F 系统处理这类物体的经典招数是在频谱面做手脚。暗场法是在中心放一个不透明小圆盘挡掉零级H_dark ones(N, N); H_dark(rho 60) 0; % 半径 60 cycle/m 的挡板 [I_dark, ~] sim4f(exp(1i*0.5*double(t0.3)), dx, 632.8e-9, 100e-3, H_dark);泽尼克相衬法更进一步把零级光的相位推迟 π/2让相位梯度转化为强度衬度H_pc ones(N, N); idx rho 60; H_pc(idx) exp(1i * pi/2); % 零级相移 [I_pc, ~] sim4f(exp(1i*0.5*double(t0.3)), dx, 632.8e-9, 100e-3, H_pc);暗场输出中相位变化剧烈的边界变亮缓慢变化区域仍是黑的相衬输出则把相位分布近似线性映射为灰度看起来像“浮雕”效果。两类滤波器都不改变振幅靠的是改变频谱各成分之间的相对相位所以滤波器 H 必须是复数矩阵用 exp(1i*pi/2) 而不是乘一个实数系数。下表汇总了本节四种滤波器的传递函数与适用场景滤波器传递函数要点像面效果参数调整方向理想低通rho fc 为 1其余 0平滑振铃明显fc 控制细节保留量平滑低通0.5(1-tanh((rho-fc)/w))平滑振铃弱w 控制过渡带高通 H_hp1 - H_lp边缘亮线配合 fc 抑制低频高频提升1-exp(-(rho/fc)^2)边缘增强、背景保留fc 取特征频率暗场中心圆形挡板置 0相位梯度可见挡板半径大于零级尺寸泽尼克相衬零级区乘 exp(jπ/2)相位映射为强度相移量为 π/2 最佳4. 4F仿真参数边界与常见坑采样率、零填充、振铃与维度顺序4.1 采样间隔、频谱面半径与混叠的定量关系物面像素间距 dx 决定了可表示的最高空间频率 f_max 1/(2dx)也就是离散傅里叶变换的奈奎斯特极限。换算到频谱面这个频率对应一个物理半径rho_max lambda · f / (2 · dx)用 2.3 节的参数 λ 632.8 nm、f 100 mm、dx 5 μm 代入rho_max ≈ 6.3 mm。意思是频谱面上半径 6.3 mm 之外的地方物理上不存在真实频率成分即使滤波器在那里开了通带通过的也只是混叠造成的伪频率。这个式子很有用把 dx 减半频谱面可用的物理半径就翻倍而把 N 翻倍只改变频谱采样的疏密不改变 rho_max。设计滤波器时先算 rho_max再决定 fc 的量级能避免“滤波器写对了但位置在物理上无意义”的尴尬。判断混叠是否严重可以看频谱强度正常频谱在 rho_max 附近强度应衰减到接近本底如果频谱在矩阵边缘仍然很强说明物面采样不足此时要减小 dx而不是靠增大 N 硬撑。4.2 零填充频谱插值还是信息增加零填充是在物面矩阵外围补零把尺寸从 N 扩展到 M。常见认知误区是“补零能提高分辨率”实际上补零只做了一件事让频谱的采样间隔从 1/(N·dx) 变细到 1/(M·dx)频谱曲线被插值得更平滑峰值定位更准但两个靠得很近的真实频率分量并不会因此被分开。补零的另一好处是消除了物体边缘与数组边界之间的周期性卷绕效应因为补零后物体的镜像副本离得更远。实现时注意坐标同步扩展M 变为原 N 的两倍后新坐标向量要用新的 M 重新生成dx 不变否则滤波器与频谱会对不齐。我一般把补零写成 padarray(t, [N/2 N/2], 0)之后所有涉及 x、fx 的向量全部按 M 重建。4.3 滤波器边缘用余弦过渡消除吉布斯振铃硬截断滤波器在截止频率处传递函数从 1 跳变到 0对应空域里出现 sinc 形状的旁瓣这就是振铃的来源。想要保留低通的“截止”语义、又不让振铃掩盖像面细节可以在截止频率附近做一段余弦过渡r rho / fc; % 归一化径向频率 H_cos ones(size(rho)); band r 0.8 r 1; H_cos(band) 0.5 * (1 cos(pi * (r(band) - 0.8) / 0.2)); H_cos(r 1) 0;过渡带从 0.8·fc 开始、到 fc 结束频带宽度固定为 0.2·fc。过渡带越宽振铃越弱但截止特性越不锐利。判断振铃是否可接受看像面中高对比度边缘附近有没有明暗交替的平行条纹有的话就把过渡带放宽到 0.3·fc 再跑一次。这个余弦渐变滤波器与 3.2 节的 tanh 低通本质上是同一种思路实际选哪个取决于你想不想保留一个显式的截止频率参数。4.4 meshgrid与fft2的维度顺序最容易错位的一处fft2 沿第一维行方向变换 y 分量沿第二维列方向变换 x 分量。meshgrid(x, x) 生成的 X 矩阵每行相同即 X 沿列方向变化对应第二维Y 矩阵每列相同即 Y 沿行方向变化对应第一维。所以 meshgrid(x, x) 与 fft2 的维度角色完全一致这是推荐用法。如果换成 ndgrid(x, x)得到的第一维是 x、第二维是 y与 fft2 的角色正好相反物函数和滤波器都会发生转置错位。判别方法很简单画 meshgrid 生成的 X看它沿哪个方向变化再画 ndgrid 生成的 X两者互为转置。同一段代码里物函数、滤波器、频谱三者必须统一用同一种网格生成函数混用 meshgrid 和 ndgrid 是 4F 仿真“像面横竖颠倒”最常见的原因。另外滤波器 H 如果是通过 meshgrid(fxv, fxv) 生成并用 rho 做半径判断的H 的类型应当是 double 或复数 double不能是 logical否则与 Uf 相乘时维度虽然合法、数值却会被隐式转换相位滤波直接失效。5. 用解析解验证4F仿真艾里斑位置与涡旋光生成5.1 艾里斑第一暗环的位置比对4F 仿真的正确性可以用圆孔的夫琅禾费衍射解析解来验证半径为 R 的圆孔频谱面强度是艾里斑第一暗环半径的解析值为 r_Airy 1.22·λ·f / D其中 D 2R 是圆孔直径。仿真里找到强度剖面第一极小值的位置与解析值对比偏差应在 5% 以内R 1e-3; t double(sqrt(X.^2 Y.^2) R); Uf fftshift(fft2(ifftshift(t))); profile abs(Uf(N/21, :)); % 过中心的水平剖面 [~, imin] min(profile(N/21 : end)); % 中心右侧第一极小 r_sim xfv(imin N/2); % 仿真暗环半径 r_theory 1.22 * lambda * f / (2*R);如果 r_sim 比理论值大且剖面边缘翘起多半是物面采样不足、频谱混叠如果 profile 在中心以外没有明显的零点检查 dx 是否小到能够分辨圆孔边缘的陡变。把 R 改成 0.5 mm 和 2 mm 各跑一遍理论值和仿真值同比例变化基本可以确认整个坐标链路没有系统性错误。5.2 光栅级次、能量守恒与涡旋光正弦光栅是另一个廉价而可靠的验证对象t 0.5·(1 cos(2π·x/Λ))频谱面应在 fx ±1/Λ 处出现一对对称的亮峰。用 max 定位峰值所在列换算成频率后与 1/Λ 对比误差小于一个频率采样间隔 du 就说明频率轴标定正确。滤波前后还可以做能量核对sum(abs(Uf_filtered).^2) / sum(abs(Uf).^2) 应等于滤波器透过率平方在频谱面上的占空比低通滤波器改小 fc 时这个比值单调下降任何不单调的情况都指向滤波器矩阵与频谱错位。最后给一个能立刻看出 4F 仿真价值的进阶手法在频谱面放螺旋相位板传递函数 H exp(j·m·θ)θ 是频域极角m 是拓扑荷theta angle(fx 1i*fy); H_vortex exp(1i * 1 .* theta); % m1 的涡旋滤波 [I_vortex, ~] sim4f(t_gaussian, dx, lambda, f, H_vortex);像面会得到中心强度为零的环形光斑环径随 m 增大而增大。这个结果可以直接与轨道角动量光学的理论预期对照也是用 MATLAB 验证“频谱面相位操作如何转化为像面振幅分布”最直观的例子。把 5.1 的艾里斑、光栅级次和这里的涡旋环放在一起4F 仿真的物面采样、频率轴、滤波器坐标系三个环节就都有了独立于代码本身的定量背书。本文还有配套的精品资源点击获取