MATLAB仿真法布里-珀罗干涉仪：多光束干涉与Airy函数全解析

做光学仿真的朋友很多都有这种感觉干涉公式背得滚瓜烂熟杨氏双缝、薄膜等倾等厚用MATLAB也能画出漂亮的条纹可是一碰到“多光束干涉”书上那几行Airy函数就让人犯怵。前面几篇我们把双光束干涉、等厚干涉、衍射这些基础仿真轮着做了一遍这篇终于轮到多光束干涉了。选谁当主角法布里-珀罗干涉仪它是多光束干涉最纯粹、最经典、也最实用的一个模型。激光器选模靠它精密光谱测量靠它光通信里的可调谐滤波器还是靠它。用一个MATLAB仿真把它彻底吃透比你背十遍公式都管用。这篇到底搞什么我们把法布里-珀罗干涉仪的透射光谱用MATLAB完整仿真一遍包括一维的透射率-波长曲线、精细度随反射率的变化规律以及二维的等倾干涉圆环条纹。代码全部给出参数从物理意义出发来定不搞“调参调到图好看”那套。适合正在学光学课程、做光电方向毕业设计或者单纯想把干涉理论变成生动图像的朋友。读的过程中你会反复看到一些坑——采样密度不够导致峰值消失、单位搞混导致相位爆炸、二维成像动态范围不对导致条纹看不见——这些全是实际仿真里天天遇到的事。1. 从双光束到多光束法布里-珀罗干涉仪到底做了什么1.1 双光束干涉的“盲区”杨氏双缝也好迈克尔逊干涉仪也好我们习惯于一个图像光走了两条路两束光之间有个相位差δ然后强度是I I₀cos²(δ/2)。这个想法本身没错但它暗含一个前提——参与干涉的有效光束只有两束。实际情况往往不是这样。两片平行的镜子放在一起光进去之后会在两个镜面之间来回反射每一次反射都会从第二个镜面“漏”出一部分光。你放在后面的探测器接收到的是一个无穷序列第一束透射光、第二束、第三束……每一束都比前一束在腔内多走了两个来回。这个序列能不能直接用强度叠加不能。因为它们来自同一个光源彼此相干必须先让它们的电场复数振幅相加再取模平方才得到真正的干涉强度。我第一次做这个仿真时就踩了这个坑直接用几何级数去累加光强结果高反射率下出来的曲线是一条平滑的底噪完全没有尖峰。后来才反应过来多光束干涉的“多”不是多路强度相加而是多路振幅相加本质上是把“多路”压缩成一个复振幅几何级数求和。1.2 多光束干涉的核心物理先振幅、后强度写法很简单。设两片镜子相同振幅反射系数为r振幅透射系数为t强度反射率R |r|²强度透射率T |t|²腔内往返一次的相位差是δ。那么从第二面镜子透射出来的第m束光相对第一束光的电场振幅可以写成t²·(r²e^{iδ})ᵐ。把所有项加起来E_t E₀ t² [1 r²e^{iδ} (r²e^{iδ})² ...] E₀ t² / (1 - r²e^{iδ})这里用到了无穷等比级数求和公比是r²e^{iδ}模长小于1级数必然收敛。取模平方之后整理成强度形式就是I_t I₀ T² / [ (1 - R)² 4R sin²(δ/2) ]如果镜子没有吸收损耗T 1 - R那么I_t I₀ (1 - R)² / [ (1 - R)² 4R sin²(δ/2) ]这就是Airy函数的标准形式。注意分母里那个sin²(δ/2)是它决定了透射率随相位差做周期性变化也是它决定了峰有多尖。1.3 法布里-珀罗干涉仪两片镜子加一个腔为什么值得单独写法布里-珀罗干涉仪就是一个典型的多光束系统两片高反射镜之间夹着一层介质通常就是空气间距为d介质折射率为n入射角为θ。光在里面来回反射每一圈光程差是2nd·cosθ所以往返相位差δ (4π/λ) · n·d·cosθ注意这里有个4π而不是2π——因为光程差是2nd·cosθ相位差是(2π/λ)乘以光程差所以前面是4π。这个系数很多人写代码时会弄错后面我会专门说。从应用角度讲法布里-珀罗几乎是所有高精度光学测量设备的核心。激光谐振腔就是在两个高反镜之间放增益介质纵模选择靠的就是法布里-珀罗的选频特性天文光谱仪里的法布里-珀罗标准具能分辨极窄的谱线。这块内容放一篇单独的仿真文章非常值因为它的物理图像清晰、公式简洁同时又有一堆值得抠的数值实现细节。2. Airy函数与三个必须搞懂的参数2.1 透射光强公式一次到位别再被符号绕晕很多教材里会把Airy函数写成下面这种形式I_t I₀ / [ 1 F·sin²(δ/2) ]其中F 4R/(1-R)²叫做“精细度系数”。这个形式和上一节那个公式是等价的只是把分母里公共的(1-R)²提出来了。用这个形式有个好处讨论问题的焦点就变成了一个很干净的东西分母什么时候取得最小值sin²(δ/2) 0也就是δ 2mπ此时I_t I₀透射率100%。分母什么时候最大sin²(δ/2) 1此时透射率降到最低。这个图像非常直观相位差处在共振点所有透射光同相叠加能量全透射偏离共振点反射光把能量“憋”在腔内透射率断崖式下跌。反射率R越高分母中sin²项的“杠杆效应”越强透射峰就越细。我习惯在仿真代码里同时算出透射率和反射率两个加起来做一个自检无吸收损耗时I_t I_r I₀恒成立。一旦发现两条曲线加起来不为1说明公式抄错了。反射率的表达式可以从能量守恒直接写出来I_r I₀ · 4R·sin²(δ/2) / [ (1 - R)² 4R·sin²(δ/2) ]2.2 精细度峰有多尖分辨本领就有多强“精细度”这个概念中文教材里经常和“精细度系数”F混在一起其实它们是两码事。F是刚才那个4R/(1-R)²而精细度finesse通常用N表示定义为自由光谱范围除以峰的半高全宽在数学上可以推出N π√R / (1 - R)精细度衡量的是在一个自由光谱范围内能放下多少个半高全宽。也就是说N越大峰越尖仪器分辨相邻谱线的能力越强。我把不同R对应的数值列一下大家感受一下量级反射率R精细度N0.04普通玻璃面0.650.54.440.929.80.9561.20.98155.50.99312.6普通玻璃表面反射率只有4%精细度不到1它根本看不出“多光束”的锋利效果只有把R做到0.9以上才能称得上高精细度法布里-珀罗。这个表格在仿真里非常有用——你看到自己画的峰太胖不用瞎猜先算一下R对应的理论N是多少就知道是参数没设对还是采样不够细。2.3 自由光谱范围两个相邻级次之间的“安全距离”自由光谱范围FSR是指相邻两个透射极大之间的频率或波长间隔。这个参数决定了仪器一次能观察多宽的谱而不发生级次重叠。频率域的FSR表达式非常优美Δν_FSR c / (2nd·cosθ)注意它和波长无关只和腔长、折射率、入射角有关。换算到波长域用近似关系Δλ ≈ λ²/(2nd·cosθ)。举个例子。腔长d 0.55mm空气n 1波长550nm正入射那么FSR ≈ (550×10⁻⁹)² / (2×1×5.5×10⁻⁴) 2.75×10⁻¹⁰m也就是0.275nm。这个数字意味着如果你想在550nm附近用这台仪器测谱线两条谱线间隔超过0.275nm就不会被同一个级次搞混。再结合分辨本领ℜ m·Nm是级次m 2nd/λ 2000R 0.9时N ≈ 29.8所以ℜ ≈ 59600可以分辨的波长极限约是550nm/59600 ≈ 0.0092nm。做仿真之前先把这些手算估值写在草稿纸上后面看结果就知道自己有没有写错。3. MATLAB实操一维透射光谱扫描与参数曲线绘制3.1 先画相位轴最直接的验证第一个程序我建议从相位轴开始目的是验证Airy函数本身。相位差δ取0到4π横轴以π为单位纵轴是透射率。画多条曲线对比不同反射率的效果。% 法布里-珀罗干涉仪透射率随相位差变化 R [0.04, 0.5, 0.9, 0.98]; % 强度反射率 delta linspace(0, 4*pi, 8001); % 往返相位差单位 rad lineSpec {k-, b-, r-, m-}; legStr cell(1, 4); figure(Color, white, Position, [100 100 900 600]); hold on; for i 1:length(R) % Airy函数无吸收损耗 T (1 - R(i)).^2 ./ ((1 - R(i)).^2 4*R(i).*sin(delta/2).^2); plot(delta/pi, T, lineSpec{i}, LineWidth, 1.6); legStr{i} sprintf(R %.2f, R(i)); end hold off; xlabel(\delta / \pi (rad)); ylabel(透射率 I_t / I_0); ylim([0, 1]); legend(legStr, Location, northeast); grid on; set(gca, FontName, Times New Roman, FontSize, 12);运行之后应该看到R 0.04时曲线几乎是平缓的余弦状波纹峰值靠近1但峰很宽R 0.5时峰开始变细R 0.9时已经是明显的窄峰R 0.98时峰窄到像一根根针。同时低谷处的透射率随R增大迅速趋向0。这个图里藏着两个细节。第一为什么R 0.04的透射率低谷不是0因为单次反射率太低腔内往返一次的能量损耗太小根本不足以形成有效的破坏性干涉。第二所有曲线在δ 2mπ处都严格等于1这是无损耗Airy函数的必然结果可以用来验证代码里的公式有没有抄错。3.2 换成波长轴贴近真实实验的结果相位轴虽好但实验里你只能扫波长或扫腔长。我更推荐把仿真做到波长轴这样可以直接和光谱仪数据对照。参数还是用前面那套d 0.55mmn 1λ₀ 550nm正入射。扫描范围取549nm到551nm一共2nm。% 法布里-珀罗干涉仪透射光谱波长扫描 lambda0 550e-9; % 中心波长 550 nm d 0.55e-3; % 腔长 0.55 mm n 1; % 空气折射率 R 0.9; % 反射率 theta 0; % 正入射 lambda linspace(549e-9, 551e-9, 10001); % 波长扫描范围 delta 4*pi*n*d*cos(theta) ./ lambda; % 往返相位差 T (1 - R).^2 ./ ((1 - R).^2 4*R.*sin(delta/2).^2); figure(Color, white, Position, [100 100 900 500]); plot((lambda - lambda0)*1e9, T, b-, LineWidth, 1.5); xlabel(波长偏移 \lambda - \lambda_0 (nm)); ylabel(透射率); title(sprintf(法布里-珀罗透射谱 d %.3f mm, R %.2f, d*1e3, R)); grid on;代码里用的是向量化操作lambda是一个10001个元素的数组delta和T也是同样长度的数组整个过程只需要一次内存分配没有for循环。这是MATLAB仿真里最基本的性能意识。运行结果里应该看到大约7个透射峰因为扫描范围2nm除以FSR 0.275nm约等于7.3。峰与峰之间的间隔并不是严格的等距——这个现象在窄范围内看不明显但如果把扫描范围拉宽到20nm就会看到波长间隔略有变化因为FSR本身就是波长的函数。3.3 从仿真结果反推精细度动手自检一次画完漂亮的峰别急着收工。我自己每次做完都会做一步验证从仿真曲线里量峰的位置和宽度反推精细度看看和理论值是否吻合。如果对不上十有八九是采样密度不够。用islocalmax找峰再用半高交叉点估FWHM% 在整段数据上找透射峰 TF islocalmax(T, MinProminence, 0.5); locs find(TF); if length(locs) 2 % 用前两个峰的位置估 FSR FSR_lambda abs(lambda(locs(2)) - lambda(locs(1))); end % 挑第一个峰量它的半高全宽 pk locs(1); % 左右两侧找 T 第一次降到 0.5 的位置 idxL find(T(1:pk-1) 0.5, 1, last); idxR pk find(T(pk1:end) 0.5, 1, first); fwhm_lambda abs(lambda(idxR) - lambda(idxL)); % 实验精细度 finesse_exp FSR_lambda / fwhm_lambda; finesse_theory pi*sqrt(R)/(1 - R); fprintf(FSR %.4f nm\n, FSR_lambda*1e9); fprintf(FWHM %.4f nm\n, fwhm_lambda*1e9); fprintf(精细度(实验) %.2f, 精细度(理论) %.2f\n, ... finesse_exp, finesse_theory);这里有个选择用0.5作为半高阈值是合理的因为无损耗情况下峰高正好是1。如果用的是有吸收损耗的模型峰顶不到1就得用峰高的一半做阈值代码要相应改成T(pk)*0.5。采样密度会对这个验证产生直接影响。理论上峰宽是FSR/NR 0.9时约0.0092nm而我用了10001个点扫2nm每点间距0.0002nm一个峰内大约46个点反推出来的精细度误差在百分之一量级。如果你偷懒只取500个点每点间距0.004nm一个峰里连3个点都没有半高交叉点就找不到FWHM会严重偏高。这类采样问题在二维仿真里会更致命。3.4 顺手加一条反射曲线能量守恒的天然检验Airy公式算完透射反射可以直接用1减去透射得到无吸收时。把两条曲线画在同一个图里你会发现一个很有趣的现象透射峰的位置正好对应反射谷而且是严格互补的。Rr 4*R.*sin(delta/2).^2 ./ ((1 - R).^2 4*R.*sin(delta/2).^2); figure(Color, white, Position, [100 100 900 500]); plot((lambda - lambda0)*1e9, T, b-, LineWidth, 1.5); hold on; plot((lambda - lambda0)*1e9, Rr, r--, LineWidth, 1.2); hold off; xlabel(波长偏移 (nm)); ylabel(强度); legend({透射率, 反射率}, Location, northeast); grid on;这个图在汇报展示时会非常加分因为它直观展示了能量守恒无论反射率多高共振波长处所有能量都穿过腔体反射为零。这也就是法布里-珀罗能当“透射式梳状滤波器”的原因。4. 二维等倾圆环条纹把一维谱线变成真实实验图4.1 为什么是同心圆环等倾干涉的几何关系真实法布里-珀罗实验中如果用扩展光源照明并在干涉仪后面加一个会聚透镜在焦平面上看到的是明暗相间的同心圆环。原因不复杂对一个固定的腔长透射率只依赖入射角θ。凡是入射角相同的光相位差相同透射率就相同。由于系统具有旋转对称性同一入射角在焦平面上对应一个圆环。问题是θ和半径r怎么对应。假设透镜焦距为f光线以角度θ入射经透镜聚焦后离焦点的横向距离近似是r f·tanθ小角度下r ≈ f·θ。所以相位差可以写成δ(r) (4πnd/λ) · cos( arctan(r/f) ) ≈ (4πnd/λ) · (1 - r²/(2f²))这个表达式说明中心r 0时相位差最大向外r增加相位差单调减小。由于每变化2π透射率就经历一个周期所以从中心向外条纹是一圈圈往外排的环而且相邻环的间距会越来越小——因为r²项是二次的。注意这跟牛顿环不一样。牛顿环是等厚干涉条纹代表空气劈尖的等厚度线这里的圆环是等倾干涉每个环代表一个入射角。两者看起来都是同心圆但物理来源完全不同。4.2 二维仿真的实现细节二维仿真用meshgrid生成网格然后一次性算出每个像素点的透射率。为了减少计算量可以用极坐标的旋转对称性但直接用直角坐标网格更直观而且能原样展示圆形条纹。% 法布里-珀罗等倾干涉二维环形条纹 lambda0 550e-9; d 0.55e-3; n 1; R 0.95; f 0.3; % 会聚透镜焦距 0.3 m rmax 0.05; % 观察面半宽 5 cm Npix 801; x linspace(-rmax, rmax, Npix); [X, Y] meshgrid(x, x); r sqrt(X.^2 Y.^2); theta atan(r / f); % 入射角 delta 4*pi*n*d*cos(theta) / lambda0; T (1 - R).^2 ./ ((1 - R).^2 4*R.*sin(delta/2).^2); figure(Color, white, Position, [100 100 820 760]); imagesc(x*1e3, x*1e3, T); axis xy; axis tight; axis equal; colormap(gray); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(FP等倾圆环条纹 d%.3fmm, R%.2f, d*1e3, R));这段代码的核心是把角度θ和位置r通过atan(r/f)关联起来而不是用近似式。虽然atan比直接r/f稍慢一点但在这个尺寸的网格上完全可以忽略。用精确关系的好处是即使你为了看更多条纹而把rmax调大入射角接近大角度时结果依然正确。801×801个点double类型内存大约是5MBMATLAB完全扛得住。如果想出更高清的图把Npix加到1601内存会涨到20MB多数电脑也没问题。2048以上就要注意了数值计算里“看起来没多大”的网格翻倍一次内存就翻四倍。4.3 看看中心剖面怎么确认画的是对的二维图好看是好看但没法精确定量。我习惯把中心一行的数据单独画出来看成一维曲线来分析。这样既能看条纹间距又能检查采样是否平滑。figure(Color, white, Position, [100 100 900 400]); plot(x*1e3, T(ceil(Npix/2), :), b-, LineWidth, 1.2); xlabel(x (mm)); ylabel(透射率); title(中心行剖面); grid on;R 0.95时剖面图会看到一串极窄的尖峰峰值高度为1尖峰之间几乎贴零。如果你用的是R 0.2你会发现峰很宽、峰谷不落零整体对比度低。这个差别清清楚楚地告诉你低反射率的法布里-珀罗没有实用价值它和普通双光束干涉没什么本质区别。参数敏感性方面腔长d和反射率R的影响我之前说过这里补一个透镜焦距f的影响f越大同样的入射角对应的横向位置越靠外所以条纹整体变疏f越小条纹越密集。如果你用白色光源做宽带仿真还会看到不同波长的环半径不同形成彩色晕圈但那是另一个话题这里不展开。5. 做这个仿真我踩过的坑查错速查表5.1 采样密度不足透射峰“失踪”或变矮这是所有法布里-珀罗仿真里最高频的问题。高反射率下峰宽极窄比如R 0.98时FWHM在2π相位内只有约0.04rad如果你在相位轴上只取了100个点每个点间距0.06rad大概率一个峰都采不到曲线会变得“矮胖”甚至基本平躺。判断原则不复杂先算理论N π√R/(1-R)再看单周期2π相位内的采样点数至少要有5N个点建议10N以上。R 0.98时N ≈ 155单周期至少需要1500个点我在仿真里取了8001个点覆盖4π即单周期4000点绰绰有余。现象可能原因解决办法透射峰高度达不到1采样点太少没有采到真正的峰值加采样密度按10N规则重取峰看起来特别胖R设得太低或相位范围算错核对R核对δ公式中的4π因子低谷不为0R太低或镜面有吸收损耗低R无吸收时用能量守恒自检二维图几乎全黑或全白动态范围太大线性色标看不清用开方/对数色标或单独看剖面相位出现Inf或NaN波长数组出现0或负数检查波长范围和单位转换5.2 相位差计算4π因子和单位别搞混我最常看学生代码里犯的错是把δ写成(2*pi*d)./lambda少乘了一个2。光程差是2nd·cosθ相位差是(2π/λ)乘以光程差所以前面是4π。如果你定义的是“单程相位差”公式又不一样关键是写代码前明确自己用的定义然后在共振条件δ 2mπ下验算一遍。单位问题也很阴险。MATLAB的sin函数用弧度但很多从C语言转过来的朋友习惯给sin传角度一传就直接炸。波长用米还是纳米扫频用频率还是角频率我自己的做法是全程用国际单位只在画图时再把横轴转换成nm或THz显示这样最不容易出错。5.3 二维显示动态范围太大线性色标“骗”了你的眼睛R 0.98时峰值和背景的强度比可以达到上万倍。线性colormap下最亮的尖峰会把整个色标撑满背景看起来就是纯黑一片好像只有几条细线实际信息被动态范围吃掉了。这不是仿真算错了是显示问题。处理方法有三类用sqrt(T)或T.^(0.3)做幂律拉伸让弱信号也能看见用对数色标imagesc(log10(T eps))把动态范围压缩或者干脆不靠色标靠中心剖面图定量看。报告里如果只贴二维图建议明确标注“强度经过非线性变换显示”避免误导读者。5.4 性能优化向量化是MATLAB的基本功二维仿真天然适合向量化用meshgrid生成网格后整个T就是一个大矩阵运算没有任何循环。但如果你习惯写嵌套for循环for i 1:N for j 1:N theta atan(sqrt(X(i,j)^2 Y(i,j)^2) / f); T(i,j) (1-R)^2 / ((1-R)^2 4*R*sin(delta(i,j)/2)^2); end endN 801时这个循环要跑64万次MATLAB的脚本模式下会卡到你怀疑人生。向量化版本用一套矩阵运算就完成速度差两个数量级。这不是技巧问题是写MATLAB仿真最基本的习惯。你可以记住一个原则只要看到矩阵计算可以用点运算完成就不要用循环。5.5 如何判断结果“看起来对”最后给一个经验法则。做完仿真不要只盯着图觉得“挺好看”按下面三步自检手算FSR和理论精细度和仿真结果对比误差应在几个百分点以内检查透射峰的位置是否满足共振条件δ 2mπ换算成波长后与手算的级次m一致检查透射与反射曲线是否互补无吸收时这是Airy函数推导里的内禀约束。这三步过了仿真结果的置信度就很高了。后续做再复杂的扩展比如加损耗、加高斯光束、扫描入射角都是在这个可信地基上盖楼。6. 从理想仿真到真实系统几个值得尝试的扩展方向6.1 加入镜面吸收损耗峰值立刻“秃”给你看实际高反镜不可能无损多少有点吸收。设单程吸收率A那么T 1 - R - A。透射公式中的分子要换成T²分母仍然是(1-R)² 4R·sin²(δ/2)。于是共振峰的峰值透射率变成T_peak T² / (1-R)² (1 - R - A)² / (1 - R)²举个例子R 0.98A 0.01峰值透射率 (0.01/0.02)² 0.25。也就是说反射率做得越高一点点吸收就会让透射峰从1掉到25%。这个问题在真实激光腔设计中非常要命——高反镜的损耗直接决定激光器的输出功率和效率。仿真里加这项不需要改多少代码但对理解实际系统帮助巨大。6.2 入射角连续分布从理想单色平面波走向真实光源真实光源不是严格的平面波也不是严格的单色光。光源有一定角分布时不同入射角对应不同δ探测器接收到的是把这些Airy函数做加权平均。结果往往是峰被展宽、对比度下降。你可以写一个极简单的扩展用一个高斯分布生成θ的权重然后对T(θ)加权求和。这一步做完你会理解为什么法布里-珀罗仪器对光源准直度要求那么高。6.3 从Airy函数走向多层膜转移矩阵法法布里-珀罗干涉仪本质上就是“两个界面的多光束干涉”。如果镜子换成多层介质膜比如十几层高低折射率交替的镀膜就必须用转移矩阵法TMM来算。TMM可以看成Airy公式的推广每层介质是一个2×2矩阵整叠膜就是这些矩阵的乘积。光学镜片的增透膜、高反膜、截止滤光片全是从这套理论来的。如果你对这块感兴趣建议先把本文的法布里-珀罗仿真吃透因为TMM里所有关键概念——相位厚度、反射率、透射率、能量守恒——都能在法布里-珀罗模型里找到直观对应。等哪天有空我可以再把TMM的MATLAB仿真单独写一篇那又是一个能玩很久的方向。做这个仿真最大的收获其实不是学会画了几张漂亮的图而是彻底理解了“为什么激光器要用高反镜”“为什么法布里-珀罗能当滤波器”这些工程问题的底层逻辑。仿真这个东西就是这样公式摆在那儿你可能无感但一旦自己亲手在代码里把R从0.04改到0.98看着透射峰从胖墩墩变成一根根细针那种感受比背十遍概念都来得扎实。下次调试光学系统时你会下意识地先去估算FSR和精细度再去怀疑仪器——这种直觉就是仿真带给你的最大回报。