简介本资源是一套面向电子信息、通信工程及生物医学工程等专业本科生的毕业设计实践材料聚焦光学相干断层扫描OCT成像原理的MATLAB仿真实现解决课程设计与毕设中缺乏可运行、可复现OCT系统建模案例的痛点。压缩包共135个文件含109个核心MATLAB脚本如光谱域OCT信号生成、干涉图重建、A/B-scan图像合成等关键模块、12张结果可视化PNG图、4个预置参数MAT数据、3份PDF格式的设计文档与研究报告以及配套说明文本与Markdown使用指南整体大小仅4.92MB轻量易部署。资源已获47人学习下载代码经实测稳定运行覆盖从光源建模、干涉信号仿真、频谱处理到图像重构的完整OCT链路附带清晰函数调用关系与参数注释特别适合初学者理解OCT物理机制也便于教师快速构建教学演示原型或学生开展二次开发拓展功能。 做OCT仿真毕设的人最开始几乎都会卡在同一个地方光学原理看了不少但打开MATLAB不知道从哪里写起。要么是搜了一堆讲OCT成像原理的PPT要么是找到了论文里的公式但公式和代码之间的路怎么走没人讲透。我当年做这个题目也一样光学教材翻了三天直觉上觉得该用干涉公式但一落到代码就蒙了。后来把整条链路拆成光源-干涉-重建三个模块思路才真正打开。这篇内容就把这条链路完整捋一遍频域OCT仿真在MATLAB里到底怎么建模、怎么实现、怎么验证。不管你是正在做相关毕业设计还是想用MATLAB实现光学仿真但找不到切入点都能按这套流程走通。我会把核心代码贴在正文里也会把那些跑出来结果不对时最容易被忽视的细节单独拎出来讲。至于文档结构和工程组织方式最后也提一下因为你拿到一个OCT仿真工程包的时候第一步就是看懂它的文件结构。1. 先搞清楚OCT仿真到底在“仿”什么1.1 物理过程拆解从迈克尔逊干涉说起OCT的中文名是光学相干层析成像本质上可以理解成光学版本的超声成像。超声成像靠的是声波在组织界面上的反射回波OCT靠的是光在样品内部不同折射率界面上的背向反射散射。光速太快没法像超声那样直接计时所以OCT用了一个非常巧妙的手段——低相干干涉。整个系统是一个迈克尔逊干涉仪宽谱光源发出的光经过分束器一路打到参考镜一路打到样品上。参考镜的反射光和样品内部的背向散射光在探测器位置相干叠加。只有当参考臂和样品臂的光程差小于光源的相干长度时两者才能发生明显的干涉。所以通过扫描参考镜的位置时域OCT或者分析干涉光谱的调制频率频域OCT就能还原出样品内部不同深度处的反射强度分布。在MATLAB仿真里我们不需要真的模拟光的传播过程也不需要做光线追迹只需要把这个干涉过程的数值模型写出来。核心就一个公式I(k) S(k) * |r_R Σ r_n e^(i·2k·d_n)|²展开之后会有三项参考镜的直流项、样品各层的自相干项、以及参考光与样品光之间的交叉干涉项。其中交叉干涉项才是真正包含深度信息的项它的形式是I_cross(k) 2·S(k)·r_R·Σ r_n·cos(2k·d_n)这里k是波数2π/λd_n是第n层与参考镜之间的光程差r_R是参考镜反射系数r_n是样品第n层的反射系数。这个cos项的振荡频率正比于深度d_n所以对I(k)做傅里叶变换就能把不同深度的反射信息在频率域也就是深度域分离出来。1.2 时域与频域为什么仿真首选频域方案OCT有两大类实现方式时域OCTTD-OCT和频域OCTFD-OCT。两者物理原理一样但信息获取方式完全不同。时域OCT需要机械移动参考镜在每个参考镜位置记录干涉信号的强度包络得到一个深度点。要在深度方向成像就得让参考镜快速往复扫描机械结构复杂成像速度慢。而且机械扫描的精度直接决定成像质量仿真的时候要模拟这个扫描过程代码里得嵌套大量循环效率低。频域OCT则不需要移动参考镜。它用光谱仪探测干涉光谱I(k)然后对I(k)做一次傅里叶逆变换就能得到整个深度方向的A-scan信号。也就是说一次测量就能获得所有深度信息速度快了几个数量级。灵敏度方面频域OCT也比时域OCT高很多这就是所谓的Fellgett优势。在MATLAB仿真里频域方案的优势更加明显没有循环扫描参考镜的过程只需要构造一个干涉光谱数组然后一次ifft完事。代码量少、容易调试、物理含义清晰。所以做毕设仿真除非题目明确要求实现时域OCT否则直接选频域方案。提示频域OCT在文献里有两种叫法光谱域OCTSD-OCT和扫频源OCTSS-OCT。前者用宽谱光源加光谱仪后者用扫频激光器加单点探测器。但在仿真层面两者数学模型几乎完全一样都是对I(k)做傅里叶变换。1.3 仿真模块划分光源、样品、干涉、重建在开始写代码之前先把整个仿真系统拆成五个模块。这一步非常重要因为毕设仿真最怕的就是把所有东西堆在一个脚本里后面调试起来根本无从下手。模块物理系统仿真建模内容MATLAB实现光源SLD超发光二极管 / 宽谱激光器高斯型光谱S(k)一个高斯函数分光与合束迈克尔逊干涉仪干涉光谱叠加向量运算样品被测组织 / 多层膜结构深度-反射系数分布数组d和r探测器光谱仪波数域采样k的等差数列重建信号处理傅里叶逆变换A-scan提取ifft这样拆完之后每个模块都是一个数组或者几个数组的运算代码逻辑一目了然。后面不管是换光源参数、改样品结构还是调采样点数都只需要改对应模块不用动其他代码。2. 频域OCT在MATLAB中的完整实现过程2.1 光源模型参数怎么定、代码怎么写光源模型是整个仿真中最简单的部分但参数的物理含义必须搞清楚。OCT中常用的光源是超辐射发光二极管SLD它的光谱近似为高斯形状。两个关键参数是中心波长λ_c和半高全宽Δλ。中心波长决定了成像的穿透深度和分辨率基准。生物组织OCT常用的波段是850nm和1300nm。850nm分辨率更好但组织散射强穿透浅1300nm穿透更深适合皮肤、血管等 imaging。毕设仿真里选850nm就好各方面参数都比较好算。半高全宽Δλ决定了轴向分辨率。轴向分辨率公式为δz (2·ln2 / π) · (λ_c² / Δλ) ≈ 0.44 · λ_c² / Δλ举个例子λ_c 850nmΔλ 50nm分辨率就是δz 0.44 × (850×10⁻⁹)² / (50×10⁻⁹) ≈ 6.36μm这个分辨率足以分辨细胞层面的结构。同样中心波长下带宽越宽轴向分辨率越高。但需要注意带宽变宽意味着光谱仪的采样要求也变高并非带宽越大越好。在MATLAB中构造光源光谱我推荐直接在波数域均匀采样。因为后面做ifft的时候默认要求自变量是均匀间隔的。如果你在波长域均匀采样算出来的k是不均匀的直接fft会出问题。很多初学者在这个地方踩坑。直接在波数域定义的代码lambda_c 850e-9; % 中心波长 (m) delta_lambda 50e-9; % 光谱半高全宽 (m) N 4096; % 采样点数 % 波数域均匀采样 k_min 2*pi / (lambda_c 3*delta_lambda); k_max 2*pi / (lambda_c - 3*delta_lambda); k linspace(k_min, k_max, N); lambda 2*pi ./ k; dk k(2) - k(1); % 高斯型光源光谱 S_k exp( -((lambda - lambda_c) / (delta_lambda/2)).^2 );这里取3倍的半高全宽作为光谱范围已经覆盖了高斯光谱的绝大部分能量。S_k是一个长度为N的列向量取值范围0到1代表光源在每个波数分量上的强度。2.2 样品模型多层结构如何用代码表达样品模型就是定义一组深度-反射系数对。在频域OCT仿真中最常见的做法是把样品建模成离散的反射面集合。每个反射面用一个光程差d_n和一个反射系数r_n表示。需要注意的是这里的d_n是样品内部反射面与参考镜之间的光程差而不是样品的物理深度。因为光在组织里传播要考虑折射率光程 物理深度 × 折射率。而且OCT是反射式成像光要走一个来回所以实际光程差是2倍的单程光程。举个例子假设样品表面在零光程差位置和参考镜等光程内部深度z_n处有一个折射率n的界面那么d_n 2 × n × z_n如果n1.4z_n50μm那么d_n 2×1.4×50μm 140μm。在MATLAB中定义一个多层样品% 样品各反射面的光程差 (m) 和反射系数 d [30, 65, 120] * 1e-6; % 三个反射面 r [0.6, 0.3, 0.8]; % 对应的反射系数幅度这里r是反射系数的幅度取值0到1。实际生物组织中反射率通常很低比如0.01到0.1的量级但仿真中可以设置大一点方便观察结果。2.3 干涉信号计算与傅里叶重建有了光源、参考镜、样品之后就能算干涉光谱了。按照第一节的公式交叉干涉项是每个反射面的cos项加权求和% 参考镜反射系数 r_R 1; % 干涉光谱直流项 交叉干涉项 I_k S_k * r_R^2; % 直流项 for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2 * k * d(idx)); end算完之后I_k就是光谱仪接收到的干涉光谱。这是一个长度为N的实数数组。接下来就是关键的傅里叶重建% 傅里叶逆变换得到A-scan A_scan abs(ifft(I_k)); % 深度轴只取前N/2个点 z_axis (0:N/2-1) * (pi / (N * dk));这里有两个细节必须说明。第一为什么用ifft而不是fft从纯数学角度对实数信号做fft和ifft后取绝对值得到的结果几乎一样只差一个缩放因子。但从物理定义上I(k)是频域信号要变回空间域深度域对应的是傅里叶逆变换所以用ifft更严谨。这个细节在毕设答辩里如果被问到能答上来会加分。第二深度轴z_axis的计算。波数采样间隔是dk总采样点数是N那么最大可测深度由奈奎斯特采样定理决定z_max π / (2·dk)而z_axis的步长是π/(N·dk)对应FFT的最小频率分辨率。两者配合恰好覆盖从0到z_max的N/2个点。如果你改了N或者k的范围z_axis会自动跟着变不需要手动改。2.4 完整的A-scan仿真代码把上面几段拼起来就是一个完整的频域OCT单点仿真%% 频域OCT A-scan仿真 clear; clc; % ---- 光源参数 ---- lambda_c 850e-9; % 中心波长 (m) delta_lambda 50e-9; % 半高全宽 (m) N 4096; % 采样点数 % ---- 波数域采样 ---- k_min 2*pi / (lambda_c 3*delta_lambda); k_max 2*pi / (lambda_c - 3*delta_lambda); k linspace(k_min, k_max, N); lambda 2*pi ./ k; dk k(2) - k(1); % ---- 光源高斯光谱 ---- S_k exp( -((lambda - lambda_c)/(delta_lambda/2)).^2 ); % ---- 样品与参考镜 ---- r_R 1; d [30, 65, 120] * 1e-6; r [0.6, 0.3, 0.8]; % ---- 干涉光谱 ---- I_k S_k * r_R^2; for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2 * k * d(idx)); end % ---- 傅里叶逆变换 ---- A_scan abs(ifft(I_k)); % ---- 深度轴与显示 ---- z_axis (0:N/2-1) * (pi/(N*dk)); figure; plot(z_axis*1e6, A_scan(1:N/2), LineWidth, 1.2); xlabel(Depth (μm)); ylabel(Amplitude); title(FD-OCT A-scan); grid on;跑完这段代码你会看到三个峰值分别位于约30μm、65μm和120μm的位置峰的高度正比于反射系数r。这就是一个最基本的OCT深度剖面。3. 从一维A-scan到二维B-scan成像3.1 横向扫描的仿真思路单个A-scan只反映样品某一点沿深度方向的反射信息。要得到二维断层图像B-scan需要对样品进行横向扫描也就是让光束在样品表面移动在每个横向位置记录一条A-scan。横向扫描的结果拼在一起就是一张二维OCT图像。仿真里做横向扫描本质上是一个循环对每个横向位置x定义该位置对应的样品结构d和r计算A-scan存入B_scan矩阵。关键是让不同横向位置的样品结构有差异这样才能在图像里看到有意义的形态。比如模拟一个样品中间区域比两侧多了一层反射界面%% 二维B-scan仿真 W 256; % 横向扫描点数 B_scan zeros(N/2, W); % 预分配矩阵 for x 1:W % 根据横向位置定义样品结构 if x W/3 x 2*W/3 d [30, 65, 120]*1e-6; r [0.6, 0.3, 0.8]; else d [30, 65]*1e-6; r [0.6, 0.3]; end % 计算干涉光谱 I_k S_k * r_R^2; for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2*k*d(idx)); end % 取单边A-scan a abs(ifft(I_k)); B_scan(:, x) a(1:N/2); end横向分辨率主要由聚焦光斑大小决定仿真中可以通过直接调整不同横向位置的样品结构来模拟。如果横向扫描间距设得比横向分辨率小相邻位置的样品结构会互相覆盖这在仿真里可以做一个横向高斯平滑但初学者不用一上来就加这个复杂度。3.2 图像显示与后处理B-scan数据本身是强度值但OCT信号的动态范围很大从最强的直流项到最弱的样品信号可能差好几个数量级。直接显示线性强度弱信号会被完全淹没什么都看不见。正确的做法是转换成分贝刻度。以峰值强度作为0dB参考其他信号用相对分贝值表示% 分贝变换 B_scan_db 20 * log10(B_scan / max(B_scan(:)) eps); % 显示 figure; imagesc(1:W, z_axis(1:N/2)*1e6, B_scan_db, [-40, 0]); xlabel(Lateral position); ylabel(Depth (μm)); title(OCT B-scan (dB)); colormap(gray); axis xy; colorbar;这段代码里imagesc的第四个参数[-40,0]表示显示动态范围也就是只显示相对峰值衰减40dB以内的信号。低于-40dB的信号显示为黑色高于0dB的显示为白色。这个阈值可以根据实际效果调整。axis xy是必须加的因为imagesc默认的Y轴方向是从上到下增大加了axis xy之后深度0在图像顶部深度增大方向朝下这才符合OCT图像的常规显示习惯。3.3 后处理里最容易忽略的细节做B-scan显示的时候有几个细节如果不注意出来的图会很丑或者误导性强。第一直流项对应的零深度峰会非常亮它会在图像顶部形成一条亮线掩盖浅层信号。仿真里如果样品表面刚好在零光程差位置这个现象会很明显。解决办法是让样品表面稍微偏离零光程差位置比如把第一个反射面设在20μm而不是0μm。这样直流峰和样品信号在空间上分开了。第二FFT产生的镜像伪影。对实数干涉光谱做ifft之后信号关于零光程差位置对称所以你会在负深度方向看到一个镜像副本。显示的时候只取前N/2个点正深度部分就是为了避开镜像。但如果你把样品信号放在了靠近零光程差的位置镜像峰和实像峰可能发生重叠导致图像看起来多了一些结构。所以样品深度设定要避开零点附近。第三如果图像上的层结构看起来是弯曲的、倾斜的那不是仿真错了而是横向扫描过程中样品结构在变化。这个可以用来模拟非平坦表面的样品。4. 仿真结果的正确性验证与参数分析4.1 峰值位置验证第一步必须做的事写完仿真代码第一件事不是急着做二维成像而是验证A-scan的正确性。方法很简单设定已知的样品结构看重建的峰值位置和幅度是否与设定一致。比如你设置了d [30, 65, 120]μmr [0.6, 0.3, 0.8]那么A-scan里在30μm、65μm、120μm处应该有三个峰峰的比例接近0.6:0.3:0.8。需要说明的是由于ifft输出是复数取模而且光谱形状高斯包络影响着每个峰的幅度所以峰的绝对高度不严格等于r的值但相对比例基本保持。如果峰值位置偏了优先检查深度轴z_axis的计算是否正确如果峰值形状不对、出现展宽优先检查是不是直接在非均匀k域做了fft如果峰值出现折叠、跑到图像边缘说明最大成像深度不够需要减小dk也就是增加N或者扩大k范围。4.2 轴向分辨率与带宽的关系验证验证完位置和幅度下一步验证轴向分辨率。轴向分辨率的理论公式是δz ≈ 0.44·λ_c²/Δλ。这个公式可以用仿真来验证。具体操作保持中心波长不变分别用Δλ 25nm、50nm、100nm跑仿真在相同深度位置设置一个孤立反射面然后从A-scan里读取峰的半高全宽。几个带宽对应的理论分辨率半高全宽Δλ理论轴向分辨率半高全宽对应峰宽仿真量25nm约12.7μm本文还有配套的精品资源点击获取 频域OCT仿真MATLAB实现：从原理到代码全解析
简介本资源是一套面向电子信息、通信工程及生物医学工程等专业本科生的毕业设计实践材料聚焦光学相干断层扫描OCT成像原理的MATLAB仿真实现解决课程设计与毕设中缺乏可运行、可复现OCT系统建模案例的痛点。压缩包共135个文件含109个核心MATLAB脚本如光谱域OCT信号生成、干涉图重建、A/B-scan图像合成等关键模块、12张结果可视化PNG图、4个预置参数MAT数据、3份PDF格式的设计文档与研究报告以及配套说明文本与Markdown使用指南整体大小仅4.92MB轻量易部署。资源已获47人学习下载代码经实测稳定运行覆盖从光源建模、干涉信号仿真、频谱处理到图像重构的完整OCT链路附带清晰函数调用关系与参数注释特别适合初学者理解OCT物理机制也便于教师快速构建教学演示原型或学生开展二次开发拓展功能。 做OCT仿真毕设的人最开始几乎都会卡在同一个地方光学原理看了不少但打开MATLAB不知道从哪里写起。要么是搜了一堆讲OCT成像原理的PPT要么是找到了论文里的公式但公式和代码之间的路怎么走没人讲透。我当年做这个题目也一样光学教材翻了三天直觉上觉得该用干涉公式但一落到代码就蒙了。后来把整条链路拆成光源-干涉-重建三个模块思路才真正打开。这篇内容就把这条链路完整捋一遍频域OCT仿真在MATLAB里到底怎么建模、怎么实现、怎么验证。不管你是正在做相关毕业设计还是想用MATLAB实现光学仿真但找不到切入点都能按这套流程走通。我会把核心代码贴在正文里也会把那些跑出来结果不对时最容易被忽视的细节单独拎出来讲。至于文档结构和工程组织方式最后也提一下因为你拿到一个OCT仿真工程包的时候第一步就是看懂它的文件结构。1. 先搞清楚OCT仿真到底在“仿”什么1.1 物理过程拆解从迈克尔逊干涉说起OCT的中文名是光学相干层析成像本质上可以理解成光学版本的超声成像。超声成像靠的是声波在组织界面上的反射回波OCT靠的是光在样品内部不同折射率界面上的背向反射散射。光速太快没法像超声那样直接计时所以OCT用了一个非常巧妙的手段——低相干干涉。整个系统是一个迈克尔逊干涉仪宽谱光源发出的光经过分束器一路打到参考镜一路打到样品上。参考镜的反射光和样品内部的背向散射光在探测器位置相干叠加。只有当参考臂和样品臂的光程差小于光源的相干长度时两者才能发生明显的干涉。所以通过扫描参考镜的位置时域OCT或者分析干涉光谱的调制频率频域OCT就能还原出样品内部不同深度处的反射强度分布。在MATLAB仿真里我们不需要真的模拟光的传播过程也不需要做光线追迹只需要把这个干涉过程的数值模型写出来。核心就一个公式I(k) S(k) * |r_R Σ r_n e^(i·2k·d_n)|²展开之后会有三项参考镜的直流项、样品各层的自相干项、以及参考光与样品光之间的交叉干涉项。其中交叉干涉项才是真正包含深度信息的项它的形式是I_cross(k) 2·S(k)·r_R·Σ r_n·cos(2k·d_n)这里k是波数2π/λd_n是第n层与参考镜之间的光程差r_R是参考镜反射系数r_n是样品第n层的反射系数。这个cos项的振荡频率正比于深度d_n所以对I(k)做傅里叶变换就能把不同深度的反射信息在频率域也就是深度域分离出来。1.2 时域与频域为什么仿真首选频域方案OCT有两大类实现方式时域OCTTD-OCT和频域OCTFD-OCT。两者物理原理一样但信息获取方式完全不同。时域OCT需要机械移动参考镜在每个参考镜位置记录干涉信号的强度包络得到一个深度点。要在深度方向成像就得让参考镜快速往复扫描机械结构复杂成像速度慢。而且机械扫描的精度直接决定成像质量仿真的时候要模拟这个扫描过程代码里得嵌套大量循环效率低。频域OCT则不需要移动参考镜。它用光谱仪探测干涉光谱I(k)然后对I(k)做一次傅里叶逆变换就能得到整个深度方向的A-scan信号。也就是说一次测量就能获得所有深度信息速度快了几个数量级。灵敏度方面频域OCT也比时域OCT高很多这就是所谓的Fellgett优势。在MATLAB仿真里频域方案的优势更加明显没有循环扫描参考镜的过程只需要构造一个干涉光谱数组然后一次ifft完事。代码量少、容易调试、物理含义清晰。所以做毕设仿真除非题目明确要求实现时域OCT否则直接选频域方案。提示频域OCT在文献里有两种叫法光谱域OCTSD-OCT和扫频源OCTSS-OCT。前者用宽谱光源加光谱仪后者用扫频激光器加单点探测器。但在仿真层面两者数学模型几乎完全一样都是对I(k)做傅里叶变换。1.3 仿真模块划分光源、样品、干涉、重建在开始写代码之前先把整个仿真系统拆成五个模块。这一步非常重要因为毕设仿真最怕的就是把所有东西堆在一个脚本里后面调试起来根本无从下手。模块物理系统仿真建模内容MATLAB实现光源SLD超发光二极管 / 宽谱激光器高斯型光谱S(k)一个高斯函数分光与合束迈克尔逊干涉仪干涉光谱叠加向量运算样品被测组织 / 多层膜结构深度-反射系数分布数组d和r探测器光谱仪波数域采样k的等差数列重建信号处理傅里叶逆变换A-scan提取ifft这样拆完之后每个模块都是一个数组或者几个数组的运算代码逻辑一目了然。后面不管是换光源参数、改样品结构还是调采样点数都只需要改对应模块不用动其他代码。2. 频域OCT在MATLAB中的完整实现过程2.1 光源模型参数怎么定、代码怎么写光源模型是整个仿真中最简单的部分但参数的物理含义必须搞清楚。OCT中常用的光源是超辐射发光二极管SLD它的光谱近似为高斯形状。两个关键参数是中心波长λ_c和半高全宽Δλ。中心波长决定了成像的穿透深度和分辨率基准。生物组织OCT常用的波段是850nm和1300nm。850nm分辨率更好但组织散射强穿透浅1300nm穿透更深适合皮肤、血管等 imaging。毕设仿真里选850nm就好各方面参数都比较好算。半高全宽Δλ决定了轴向分辨率。轴向分辨率公式为δz (2·ln2 / π) · (λ_c² / Δλ) ≈ 0.44 · λ_c² / Δλ举个例子λ_c 850nmΔλ 50nm分辨率就是δz 0.44 × (850×10⁻⁹)² / (50×10⁻⁹) ≈ 6.36μm这个分辨率足以分辨细胞层面的结构。同样中心波长下带宽越宽轴向分辨率越高。但需要注意带宽变宽意味着光谱仪的采样要求也变高并非带宽越大越好。在MATLAB中构造光源光谱我推荐直接在波数域均匀采样。因为后面做ifft的时候默认要求自变量是均匀间隔的。如果你在波长域均匀采样算出来的k是不均匀的直接fft会出问题。很多初学者在这个地方踩坑。直接在波数域定义的代码lambda_c 850e-9; % 中心波长 (m) delta_lambda 50e-9; % 光谱半高全宽 (m) N 4096; % 采样点数 % 波数域均匀采样 k_min 2*pi / (lambda_c 3*delta_lambda); k_max 2*pi / (lambda_c - 3*delta_lambda); k linspace(k_min, k_max, N); lambda 2*pi ./ k; dk k(2) - k(1); % 高斯型光源光谱 S_k exp( -((lambda - lambda_c) / (delta_lambda/2)).^2 );这里取3倍的半高全宽作为光谱范围已经覆盖了高斯光谱的绝大部分能量。S_k是一个长度为N的列向量取值范围0到1代表光源在每个波数分量上的强度。2.2 样品模型多层结构如何用代码表达样品模型就是定义一组深度-反射系数对。在频域OCT仿真中最常见的做法是把样品建模成离散的反射面集合。每个反射面用一个光程差d_n和一个反射系数r_n表示。需要注意的是这里的d_n是样品内部反射面与参考镜之间的光程差而不是样品的物理深度。因为光在组织里传播要考虑折射率光程 物理深度 × 折射率。而且OCT是反射式成像光要走一个来回所以实际光程差是2倍的单程光程。举个例子假设样品表面在零光程差位置和参考镜等光程内部深度z_n处有一个折射率n的界面那么d_n 2 × n × z_n如果n1.4z_n50μm那么d_n 2×1.4×50μm 140μm。在MATLAB中定义一个多层样品% 样品各反射面的光程差 (m) 和反射系数 d [30, 65, 120] * 1e-6; % 三个反射面 r [0.6, 0.3, 0.8]; % 对应的反射系数幅度这里r是反射系数的幅度取值0到1。实际生物组织中反射率通常很低比如0.01到0.1的量级但仿真中可以设置大一点方便观察结果。2.3 干涉信号计算与傅里叶重建有了光源、参考镜、样品之后就能算干涉光谱了。按照第一节的公式交叉干涉项是每个反射面的cos项加权求和% 参考镜反射系数 r_R 1; % 干涉光谱直流项 交叉干涉项 I_k S_k * r_R^2; % 直流项 for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2 * k * d(idx)); end算完之后I_k就是光谱仪接收到的干涉光谱。这是一个长度为N的实数数组。接下来就是关键的傅里叶重建% 傅里叶逆变换得到A-scan A_scan abs(ifft(I_k)); % 深度轴只取前N/2个点 z_axis (0:N/2-1) * (pi / (N * dk));这里有两个细节必须说明。第一为什么用ifft而不是fft从纯数学角度对实数信号做fft和ifft后取绝对值得到的结果几乎一样只差一个缩放因子。但从物理定义上I(k)是频域信号要变回空间域深度域对应的是傅里叶逆变换所以用ifft更严谨。这个细节在毕设答辩里如果被问到能答上来会加分。第二深度轴z_axis的计算。波数采样间隔是dk总采样点数是N那么最大可测深度由奈奎斯特采样定理决定z_max π / (2·dk)而z_axis的步长是π/(N·dk)对应FFT的最小频率分辨率。两者配合恰好覆盖从0到z_max的N/2个点。如果你改了N或者k的范围z_axis会自动跟着变不需要手动改。2.4 完整的A-scan仿真代码把上面几段拼起来就是一个完整的频域OCT单点仿真%% 频域OCT A-scan仿真 clear; clc; % ---- 光源参数 ---- lambda_c 850e-9; % 中心波长 (m) delta_lambda 50e-9; % 半高全宽 (m) N 4096; % 采样点数 % ---- 波数域采样 ---- k_min 2*pi / (lambda_c 3*delta_lambda); k_max 2*pi / (lambda_c - 3*delta_lambda); k linspace(k_min, k_max, N); lambda 2*pi ./ k; dk k(2) - k(1); % ---- 光源高斯光谱 ---- S_k exp( -((lambda - lambda_c)/(delta_lambda/2)).^2 ); % ---- 样品与参考镜 ---- r_R 1; d [30, 65, 120] * 1e-6; r [0.6, 0.3, 0.8]; % ---- 干涉光谱 ---- I_k S_k * r_R^2; for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2 * k * d(idx)); end % ---- 傅里叶逆变换 ---- A_scan abs(ifft(I_k)); % ---- 深度轴与显示 ---- z_axis (0:N/2-1) * (pi/(N*dk)); figure; plot(z_axis*1e6, A_scan(1:N/2), LineWidth, 1.2); xlabel(Depth (μm)); ylabel(Amplitude); title(FD-OCT A-scan); grid on;跑完这段代码你会看到三个峰值分别位于约30μm、65μm和120μm的位置峰的高度正比于反射系数r。这就是一个最基本的OCT深度剖面。3. 从一维A-scan到二维B-scan成像3.1 横向扫描的仿真思路单个A-scan只反映样品某一点沿深度方向的反射信息。要得到二维断层图像B-scan需要对样品进行横向扫描也就是让光束在样品表面移动在每个横向位置记录一条A-scan。横向扫描的结果拼在一起就是一张二维OCT图像。仿真里做横向扫描本质上是一个循环对每个横向位置x定义该位置对应的样品结构d和r计算A-scan存入B_scan矩阵。关键是让不同横向位置的样品结构有差异这样才能在图像里看到有意义的形态。比如模拟一个样品中间区域比两侧多了一层反射界面%% 二维B-scan仿真 W 256; % 横向扫描点数 B_scan zeros(N/2, W); % 预分配矩阵 for x 1:W % 根据横向位置定义样品结构 if x W/3 x 2*W/3 d [30, 65, 120]*1e-6; r [0.6, 0.3, 0.8]; else d [30, 65]*1e-6; r [0.6, 0.3]; end % 计算干涉光谱 I_k S_k * r_R^2; for idx 1:length(d) I_k I_k 2 * S_k .* r_R .* r(idx) .* cos(2*k*d(idx)); end % 取单边A-scan a abs(ifft(I_k)); B_scan(:, x) a(1:N/2); end横向分辨率主要由聚焦光斑大小决定仿真中可以通过直接调整不同横向位置的样品结构来模拟。如果横向扫描间距设得比横向分辨率小相邻位置的样品结构会互相覆盖这在仿真里可以做一个横向高斯平滑但初学者不用一上来就加这个复杂度。3.2 图像显示与后处理B-scan数据本身是强度值但OCT信号的动态范围很大从最强的直流项到最弱的样品信号可能差好几个数量级。直接显示线性强度弱信号会被完全淹没什么都看不见。正确的做法是转换成分贝刻度。以峰值强度作为0dB参考其他信号用相对分贝值表示% 分贝变换 B_scan_db 20 * log10(B_scan / max(B_scan(:)) eps); % 显示 figure; imagesc(1:W, z_axis(1:N/2)*1e6, B_scan_db, [-40, 0]); xlabel(Lateral position); ylabel(Depth (μm)); title(OCT B-scan (dB)); colormap(gray); axis xy; colorbar;这段代码里imagesc的第四个参数[-40,0]表示显示动态范围也就是只显示相对峰值衰减40dB以内的信号。低于-40dB的信号显示为黑色高于0dB的显示为白色。这个阈值可以根据实际效果调整。axis xy是必须加的因为imagesc默认的Y轴方向是从上到下增大加了axis xy之后深度0在图像顶部深度增大方向朝下这才符合OCT图像的常规显示习惯。3.3 后处理里最容易忽略的细节做B-scan显示的时候有几个细节如果不注意出来的图会很丑或者误导性强。第一直流项对应的零深度峰会非常亮它会在图像顶部形成一条亮线掩盖浅层信号。仿真里如果样品表面刚好在零光程差位置这个现象会很明显。解决办法是让样品表面稍微偏离零光程差位置比如把第一个反射面设在20μm而不是0μm。这样直流峰和样品信号在空间上分开了。第二FFT产生的镜像伪影。对实数干涉光谱做ifft之后信号关于零光程差位置对称所以你会在负深度方向看到一个镜像副本。显示的时候只取前N/2个点正深度部分就是为了避开镜像。但如果你把样品信号放在了靠近零光程差的位置镜像峰和实像峰可能发生重叠导致图像看起来多了一些结构。所以样品深度设定要避开零点附近。第三如果图像上的层结构看起来是弯曲的、倾斜的那不是仿真错了而是横向扫描过程中样品结构在变化。这个可以用来模拟非平坦表面的样品。4. 仿真结果的正确性验证与参数分析4.1 峰值位置验证第一步必须做的事写完仿真代码第一件事不是急着做二维成像而是验证A-scan的正确性。方法很简单设定已知的样品结构看重建的峰值位置和幅度是否与设定一致。比如你设置了d [30, 65, 120]μmr [0.6, 0.3, 0.8]那么A-scan里在30μm、65μm、120μm处应该有三个峰峰的比例接近0.6:0.3:0.8。需要说明的是由于ifft输出是复数取模而且光谱形状高斯包络影响着每个峰的幅度所以峰的绝对高度不严格等于r的值但相对比例基本保持。如果峰值位置偏了优先检查深度轴z_axis的计算是否正确如果峰值形状不对、出现展宽优先检查是不是直接在非均匀k域做了fft如果峰值出现折叠、跑到图像边缘说明最大成像深度不够需要减小dk也就是增加N或者扩大k范围。4.2 轴向分辨率与带宽的关系验证验证完位置和幅度下一步验证轴向分辨率。轴向分辨率的理论公式是δz ≈ 0.44·λ_c²/Δλ。这个公式可以用仿真来验证。具体操作保持中心波长不变分别用Δλ 25nm、50nm、100nm跑仿真在相同深度位置设置一个孤立反射面然后从A-scan里读取峰的半高全宽。几个带宽对应的理论分辨率半高全宽Δλ理论轴向分辨率半高全宽对应峰宽仿真量25nm约12.7μm本文还有配套的精品资源点击获取