简介面向物理光学、光电信息科学与工程等专业的学习者提供一套基于Matlab的多光束干涉透射光条纹仿真工具包。资源围绕不同反射率下透射光强度分布的变化展开对应物理光学教材中图4.3的典型曲线适合课程学习、实验预习、教师课堂演示、课程设计或自学复现。压缩包共4个文件仅77KB包含完整Matlab仿真脚本、生成的条纹图像、Word分析文档及txt说明脚本内置参数设置、干涉计算和绘图流程用户可按需替换波长、入射角或反射率参数快速比较不同条件下条纹宽度、对比度及半高全宽的变化。已有1596人浏览学习配套图、码、文齐全可直接运行并对照教材理解法布里-珀罗干涉仪或多光束叠加原理也可作为进一步扩展多层膜干涉仿真的起点。1. 光学仿真实操多光束干涉透射光强度分布曲线为什么反射率越高峰越尖第一次在物理光学教材的图4.3上看到多光束干涉透射光强曲线时多数人都会有个直观疑问同样是干涉为什么双光束的条纹是平滑余弦而多光束在反射率接近1时变成了细锐的“针”这个问题用MATLAB仿真最直接——写一个几十行的脚本把不同反射率下的透射光强随相位差的变化画出来就能亲眼看到条纹从宽到窄、从平缓到陡峭的整个过程。这篇笔记定位很明确面向光学工程、光电专业学生和需要复现教材图形的工程师给出完整可复现的MATLAB代码、参数选择逻辑和几个容易翻车的细节。做完之后你不但能画出和书上图4.3一致的曲线还能顺手验证精细度与自由光谱范围的关系。2. 多光束干涉公式的物理底细从平行板内反射到透射光强的Airy分布2.1 Fabry-Perot腔内多束光叠加的等比级数模型多光束干涉的经典场景是两片平行且部分反射的镜面构成的法布里-珀罗腔。一束光射入腔内后会在两个镜面之间反复反射每次透射出一部分光最终在透射方向上形成多光束叠加。关键点在于相邻两束透射光的相位差是恒定的由光程差决定等于δ 4π n d cosθ / λ其中n是腔内介质折射率d是腔长θ是腔内光线与镜面法线的夹角λ是波长。这个相位差不随反射次数变化所以所有透射光的复振幅构成一个公比为r^2 e^{iδ}的等比级数其中r是镜面振幅反射系数。级数求和后取模平方就得到标准透射强度公式I_t / I_0 1 / (1 F sin^2(δ/2))其中F 4R / (1 - R)^2R是强度反射率F被称为精细度系数coefficient of finesse。在MATLAB里我们不需要推导级数只要直接把这个公式写成数组运算即可。实际仿真中I_0通常取1纵轴变成归一化透射率。注意这里有个容易混淆的细节sin^2(δ/2)的最大值是1所以透射光强的最小值是1/(1F)而最大值恒等于1。这意味着无论反射率多高透射峰值永远都在1不会因为反射增强而变小。曾有人误以为R大了透射光应该更弱实际上在理想无吸收腔中共振波长处透射率仍是100%只是非共振波长被强烈反射曲线变“挑剔”了。2.2 反射率R如何改变透射光强分布把F表达式代入观察不同R时的曲线形状。当R 0.04时F 4×0.04 / 0.96^2 ≈ 0.17sin^2项的权重很小透射率在1和1/(10.17)≈0.85之间波动表现为类似双光束干涉的宽阔余弦条纹。当R 0.9时F 4×0.9 / 0.1^2 360此时只有当sin^2(δ/2)接近0时分母才接近1透射率才能达到峰值稍微偏离δ 2mπ透射率迅速跌落。这就是反射率对透射条纹的“锐化”作用。从数值仿真角度看R越接近1F值越大曲线峰附近的变化越剧烈。这带来一个直接后果画图时如果相位差采样点数不够峰值处的数值可能没有被采到画出来的曲线峰值会明显小于1。更严重的会出现“欠采样锯齿”让一个本该光滑的尖峰变成折线。因此在脚本设计中需要根据最大R值估算所需的最小采样点数这一点会在第4章详细展开。2.3 公式数值化MATLAB中表示相位差的正交与量纲问题MATLAB的sin函数默认输入单位是弧度。物理光学书里横轴常用相位差δ弧度有时会用δ/π作为单位让干涉级次恰好落在整数位置。仿真时建议先统一成“δ/π”显示这样曲线峰值位于δ/π 0, 2, 4, …这些偶数点上便于核查。代码中linspace(0, 8*pi, 2001)生成从0到8π的相位差数组2001个点让每个π区间约250个点对R0.9的尖峰已经够用但R0.98时就需要更多点。另一个量纲陷阱公式里用的是sin^2(δ/2)不是sin(δ)^2也不是sind。如果用角度制必须先把δ转成度或者直接用rad比较稳妥的写法是构造delta_rad然后调用sin(delta_rad/2)再点平方。代码里点运算很重要因为delta是数组正弦后的数组再平方必须用.^避免矩阵乘法报错。实际调试中遇到“错误使用 ^ 输入必须为标量”这类提示十有八九是把.^写成了^。3. 用MATLAB画透射光强曲线一个可以直接跑的脚本与逐行拆解3.1 完整脚本不同反射率下的透射光强分布这里给出一个完整可运行的MATLAB脚本它对应物理光学书上图4.3的典型画法横轴为δ/π纵轴为归一化透射强度多条曲线对应不同反射率。代码中特意保留注释和可调参数。% 多光束干涉透射光强分布仿真 clear; clc; close all; % 反射率数组覆盖低反射、中反射到高反射 R_list [0.04, 0.1, 0.5, 0.9, 0.98]; % 相位差范围0 ~ 8π每个π区间分配251个点 delta linspace(0, 8*pi, 2001); figure(Color,w, Position,[100 100 800 500]); hold on; grid on; for R R_list % 精细度系数 F 4*R / (1 - R)^2; % 透射光强归一化公式 I_t 1 ./ (1 F * sin(delta/2).^2); plot(delta/pi, I_t, LineWidth, 1.5, ... DisplayName, sprintf(R%.2f, R)); end xlabel(相位差 \delta / \pi); ylabel(透射光强 I_t / I_0); ylim([0 1.08]); legend(Location,northeast); title(不同反射率下多光束干涉透射光强分布);运行后你会得到5条曲线。R0.04的曲线几乎是余弦形状峰值宽而平缓R0.9和R0.98的曲线则在δ/π 0、2、4、6、8处出现极窄的尖峰且底座几乎贴近零。这个图形趋势和教材图4.3完全一致。3.2 逐行逻辑解释从F到I_t的点运算代码里最关键的是F 4R / (1 - R)^2。注意这里R是标量(1 - R)^2也是标量所以用普通除法没问题。但如果你将来把R做成数组就一定要改成F 4R ./ (1 - R).^2。工程上我的习惯是能点运算的地方全用点运算避免修改参数后触发维度错配。比如sin(delta/2).^2这里delta是1×2001的行向量sin(delta/2)也是同尺寸数组.^2对每个元素平方得到包含2001个值的数组。1 ./ (...)同样用点除把标量1除以每个元素。绘图部分delta/pi把横轴转成π的倍数。这样x轴刻度默认在0、1、2、3处而峰值出现在偶数刻度看起来更清晰。使用DisplayName配合legend可以自动生成图例比手动写R0.04之类的文本更省事。这里强调一点不要直接用legend({R0.04,R0.1...})这种硬编码方式否则调整R_list后需要同步修改图例容易漏改。3.3 运行结果与教材图4.3的对照要点运行脚本后重点观察两点。第一所有曲线的透射峰值都等于1无论反射率多高峰值位置都在δ/2 mπ处即δ 2mπ。第二R越大曲线的半高全宽越窄。物理上半宽度Δδ以相位差计近似为4/√F而Finesse定义为2π/Δδ因此精细度约为π√F/2。拿R0.9来算√F19精细度约为π×19/2≈29.8这意味着一个自由光谱范围内能容纳约30个可分辨的峰对应到图上就是峰值非常锋利。对照教材时如果发现书上某些曲线峰值似乎低于1那往往是印刷时纵轴归一化方式不同或者镜面存在吸收损耗我们的仿真假设无损所以峰值恒为1。4. 参数深入反射率、相位扫描范围和采样点数如何决定曲线质量4.1 反射率R从0.04到0.98条纹锐化的临界点反射率选取直接决定图形面貌。低反射区R ≤ 0.1的条纹接近余弦适合展示多光束向双光束的过渡中反射区R ≈ 0.5条纹已有明显变锐趋势高反射区R ≥ 0.9则需要处理数值上的尖峰采样问题。实际做参数扫描时我一般会取对数间隔比如0.04、0.1、0.3、0.6、0.9、0.98这样曲线在图上均匀分布不会因为R0.9和0.98靠得太近而挤压图例。如果你要复现教材某张具体图先看书上给的R值再决定是否加入额外中间值。R接近1时F值急剧增大。以R0.98为例F 4×0.98/0.02^2 9800这已经是R0.9的27倍。此时曲线的半高宽Δδ ≈ 4/√9800 ≈ 0.04弧度而单个π区间长度是π≈3.14尖峰只占1/78的范围。如果采样点每个π区间只给100点平均每隔0.0314弧度采一个点运气不好峰值点会被漏掉画出来就是一个略低于1的平顶峰。经验法则一个π区间至少设置200点对R0.98建议500点。脚本里2001点覆盖8π每个π区间250点基本够用。4.2 相位差扫描范围为什么推荐2π的整数倍仿真扫描区间选择0到2π理论上可以完整展示一个级次但物理光学教材通常会把多个级次并列好让读者理解“等间隔峰排列”。因此我推荐扫描至少4π或8π让横轴出现2到4个峰。这样也能顺带验证相邻峰的间距是否严格等于2π即自由光谱范围。注意如果用0到8π横轴长度是4个自由光谱范围R0.98的曲线会出现4个同形状的峰中间隔着一大片接近零的区域。如果你的目的是放大观察单个峰可以直接把区间缩小比如delta linspace(3pi, 4pi, 1000)只画一个峰这样能更清楚地看到FWHM的细节。有一个常见误区有人把δ直接从头到尾线性增加却忘了δ正比于频率或波长变化所以真正的实验扫描通常不是相位差线性而是波长线性。不过在理论仿真中直接以δ为自变量是最标准的不需要考虑非线性问题。如果你要模拟实际光谱仪扫描就需要把δ映射到λ例如δ 4πL/λ这时横轴会变成波长曲线形状在反射率很高时会出现非对称展宽——那是另一层话题了。4.3 采样点数与峰值插值避免“尖峰抹平”的采样量经验欠采样是这类仿真中最常见的“玄学”问题。你以为代码没错但画出来的R0.98曲线峰值只有0.85怎么调反射率都没用最后发现是delta点数太少。具体原因在于峰值在δ2mπ处是一个孤立的极窄尖点如果采样点中恰好没有δ等于2π的整数倍点那么附近的最大采样值就会低于1。解决办法有三个方向一是增加网格点数把linspace的第三个参数提到5000甚至10000二是显式把峰值点插入数组中三是用fplot匿名函数绘图让MATLAB自适应取样。我个人推荐第三种尤其在高反射率时。用fplot((d) 1./(19800sin(d/2).^2), [0 8pi])MATLAB会在峰附近自动加密取样画出的曲线不会缺角。但fplot返回的默认横轴是弧度你要单独设置xticklabels才能显示为δ/π。如果你坚持用linspace就记住这个经验点数至少为最大反射率对应的FWHM宽度的十倍所覆盖的π区间数量的20倍。简单估算对R0.98FWHM约0.04弧度扫描8π约25弧度你需要至少25/0.04 × 10 6250个点。所以2001点确实不够想画细致就得加到8000以上。这也是为什么我在第3章的代码里只写到2001点——它适合R≤0.9的展示不适合R0.98的精细分析。5. 避坑与常见问题让仿真曲线稳定可复现的五个细节5.1 现象透射峰值始终小于1最大值只有0.8左右这个问题在R0.98时尤其明显。原因是linspace采样点没有恰好落在δ2mπ上而尖峰宽度极窄采样点落在峰两侧斜坡上导致近似最大值低于真实峰值。解决方法是把delta序列中手动插入2mπ点例如在生成delta后用[0:2pi:8pi]补进数组再排序去重。也可以用fplot自适应采样。排查时先用max(I_t)看返回值如果接近1且相差小于0.001就能接受如果差很多就检查点数。5.2 现象R接近0.999时曲线变成一条直线或者峰值爆炸当R非常接近1时F 4R/(1-R)^2在R0.999时约4×0.999/1e-6 ≈ 4×10^6。理论上峰值还是1但数值计算中在远离峰值处sin^2(δ/2)如果因浮点误差出现非零极小值分母会变得极大I_t变成极小值接近0。而峰值处sin^2为0分母为1所以曲线会在极窄处跳变。如果网格点数不够可能完全没采到峰值整条曲线看上去就是一条接近于0的直线。解决方法是不要用R0.999这种极端值做展示除非你确实要演示极限精细度如果必须用就把δ范围缩小到单个峰附近且用高精度计算避免浮点截断。5.3 现象横轴标注的是弧度和书上角度/相位差不一样物理光学教材里图4.3有时横轴用“相位差δ”但单位可能是π或2π有时候直接用入射角。MATLAB默认的plot横轴只是坐标向量并不带单位。解决方法是自己控制横轴语义如果你画的是delta/pi那xlabel就要写成“δ/π”如果你直接画delta弧度就写“δ (rad)”。有些人忘记给横轴除以π导致峰间距在数值上显示为6.28而不是2看上去很奇怪。归根结底是你决定横轴是什么只要标注和绘图一致就行。5.4 现象用sind(角度)画图峰值位置和相位差对不上MATLAB里sind期望输入度数而δ是弧度。如果你用了sind(delta/2)那么只有当delta/2是整数时sin值为0峰值位置出现在delta/2180k度时也就是δ360k度≈6.28k弧度。这本身没问题但y轴公式中的sin^2却变成了以度为单位导致半高宽数值完全不同。更麻烦的是如果你混用sin和sind同一张图两条曲线宽度不可比。解决方法是全脚本统一用sin且把delta明确为弧度。偶尔有人把2π误写为360而MATLAB的pi就是π不存在角度制所以没必要用sind。5.5 现象导出的图放大后锯齿严重文章里看不清峰值这不算公式错误而是出图规范问题。锯齿有两个来源采样点不足以及LineWidth太细。采样点问题按4.3节解决线宽建议设置到2且用“Alias”关掉实际上在MATLAB figure里锯齿主要来自屏幕显示导出矢量图例如pdf或eps后锯齿会消失。我的做法是print(gcf, -dpdf, multi_beam.pdf)再用Inkscape转成png时选择高分辨率。图表符合期刊投稿要求的关键点包括图例字体不小于8pt曲线颜色按反射率从浅到深排列峰值纵坐标标到1.0坐标轴范围比数据区略大而不是紧紧贴边。6. 进阶从强度曲线提取FWHM与Finesse用解析值交叉验证仿真当你已经能画出漂亮的曲线后下一步就是量化验证测量的半高全宽是否和理论公式一致。这一步能有效排除坐标轴单位、公式系数等隐性问题。理论值很好算对透射光强公式峰值一半对应的条件为I_t 1/2也就是1 F sin^2(δ/2) 2解得sin^2(δ/2) 1/F。假设δ_0 2mπ是峰值位置那么半高半宽Δδ/2满足sin((δ_0 Δδ/2)/2) sin(Δδ/4) ≈ Δδ/4 1/√F小角近似所以全宽Δδ ≈ 4/√F。由此Finesse 2π/Δδ ≈ π√F/2。用这个解析值去对比仿真提取值如果偏差超过5%说明仿真中可能有欠采样或坐标轴换算错误。提取FWHM我习惯用MATLAB自带findpeaks然后半高交叉法。下面给一个单峰提取的例子% 选择单条高反射率曲线 R 0.9; F 4*R/(1-R)^2; delta linspace(0, 4*pi, 8000); I_t 1 ./ (1 F*sin(delta/2).^2); % 只取第二个峰附近区域避免峰值过多干扰 ind_center find(abs(delta - 2*pi) 0.01); center_val max(I_t(ind_center)); half center_val / 2; % 找左边界和右边界 left_idx find(delta 0 delta 2*pi I_t half, 1, last); right_idx find(delta 2*pi delta 4*pi I_t half, 1, first); fwhm delta(right_idx) - delta(left_idx); finesse_exp 2*pi / fwhm; finesse_theory pi*sqrt(F)/2;解释一下这段代码的思路先把8000个点的曲线算出来然后锁定δ2π处的峰。通过center_val找到峰值理论上应为1取半高阈值half。接着在左右两侧分别找到第一个低于half的点这两个点对应的横坐标之差就是FWHM。这里使用了find加方向限制避免了多个峰之间的干扰。对于R0.9理论Finesse≈29.8仿真提取值通常落在29到31之间误差主要来自线性插值而非物理偏差。如果你想把精度提高可以在left_idx和right_idx之间用interp1做线性插值精确求解I_thalf对应的δ值。工程上插值后的偏差可以压到0.1%以下。这个方法有一个隐藏价值当你不确定某个复杂光学腔的实际精细度时可以先按理想公式仿真再拿实测曲线验算。我曾经在搭建可调谐滤波器时用这个思路定位到镜片反射率低于设计值——因为提取出的Finesse只有预期的一半反解等效R约为0.82和厂家标称值差距不小最后换货解决了问题。从那以后每次仿真干涉曲线我都会强制跑一遍FWHM交叉验证哪怕只是写个教学Demo也不例外。这样做的习惯帮我避开了无数次“图看着对、数据却对不上”的尴尬。希望这篇笔记能帮你在光学仿真这条路上少走几个弯路。本文还有配套的精品资源点击获取