扫频源光学相干层析（SS-OCT）原理仿真（附matlab代码）（未完待更新）

2024-05-06 21:48
文章标签 代码 matlab 更新 仿真 原理 oct 光学 ss 相干 层析 完待 频源

本文主要是介绍扫频源光学相干层析（SS-OCT）原理仿真（附matlab代码）（未完待更新），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

   本文目的是用简单的matlab仿真语言将OCT涉及到的一些计算公式一步步还原。由简入难，从光的干涉，低相干，色散矫正，时域插值等部分进行讲解。

练习一、光的干涉。两束波长相同的单色光同向干涉。

光的波动公式可以表示为，U=A*cos(wt+p0); 其中p为初始相位。二者干涉，即为波动信号的直接叠加。

w = 5*pi;       % 角速度
t = -1:0.01:1;     % 采样时间
delta_p = 1/2 * pi; % 相位差
A1 = 2;         % 信号1的振幅
A2 = 1.5;       % 信号2的振幅
U1 = A1*cos(w*t);       % 信号1的波形
U2 = A2*cos(w*t+delta_p); % 信号2的波形
UI = U1 + U2;         %信号叠加后的波形figure(1);title('同一波长的干涉仿真');
hold on;grid on
plot(t,U1,'b');
plot(t,U2,'k');
plot(t,UI,'r');
xlabel('时间t');
ylabel('振幅A');
hold off

图1、仿真结果

可以看出，叠加后的信号，频率和源信号一致，但是幅值并不是源信号的直接相加，而是与相位差有关
修改 相位差delta_p的值，叠加后的振幅随之改变。改变规律，如下一部分练习。

练习二、同一波长在不同相位差下的干涉

p = linspace(0,5*pi,200);            %相位差
I0 = 1;                              %设两光路的强度均为1
I = I0 + I0 + 2*sqrt(I0*I0)*cos(p);  %叠加后的强度
figure;title('不用波长在同一光程差下的干涉强度');
plot(p,I);grid on
xlabel('相位差');
ylabel('干涉强度');


练习三、低相干的信号叠加。

低相干可以理解为，由于光源出射的光不是严格单色的，干涉条纹的涨落不是完全一致的。《光学原理》课本上解释，“多色的点光源出射准单色光，可以用分布在某一频率范围互不相干的许多单色成分的组合来表示，每个成分各产生一个上述的干涉图像，各点的总强度这些单色图样的强度之和.“

%%模拟低相干光源的输出
Np = 5000;
lamda0 = 1000;          %输出的中心波长
miu = 1.5;              %方差，也就是谱宽影响因子，
xx = linspace(800,1200,Np);
lamda = 1/(sqrt(2*pi)*miu) * exp(-(xx-lamda0).^2/miu^2);
%计算谱宽
valMax = max(lamda);
for idx = 2:Npif(lamda(idx-1)<=valMax/2 && lamda(idx)>valMax/2)nSt = idx;elseif(lamda(idx-1)>=valMax/2 && lamda(idx)<valMax/2)nEd = idx;end
end
delta_lamda = xx(nEd)-xx(nSt);
titllename = ['准单色光谱宽 ',num2str(delta_lamda),'nm'];
plot(xx,lamda),title(titllename)

Np = 5000;                         %采样点数
L = 1000e3;                        %光程差，nm
miu = 0.5;                         %扫频源中瞬时谱宽影响因子，也就是方差，在是%谱域系统中与光谱仪的分光能力有关
winG = gausswin(Np);               %高斯窗
lamda0 = linspace(1000,1100,Np);   %中心波长的分布
t = linspace(0,1,Np);              %一个扫频周期
lamda = zeros(Np,Np);%模拟光谱输出，强度（波长，时间）
for idx = 1:NplamdaT = 1/(sqrt(2*pi)*miu) * exp(-(lamda0-lamda0(idx)).^2/miu^2);lamda(:,idx) = winG'.*lamdaT;
endfigure
a = lamda(Np/2,:);
b = lamda(:,Np/2);
subplot(2,1,1),plot(lamda0,a);title('某一时刻输出的瞬时光谱分布')
subplot(2,1,2),plot(t,b);title('某一波长在整个扫频周期输出的强度分布')% 输出谱验证
I = zeros(Np,1);
for idx = 1:NpI0 = lamda(:,idx);I = I + I0; 
end
figure
plot(t,I);title('验证模拟的输出谱是否为高斯型')% 每一个时刻都有一系列波长在叠加干涉，也可以说一个波长在不同时刻产生叠加
I = zeros(Np,1);
for idx = 1:Npp = 2*pi*L./lamda0(idx);I0 = 2*lamda(:,idx)+2*lamda(:,idx)*cos(p);I = I + I0;
end
figure
plot(t,I);title('干涉信号波形')% 将每个瞬时输出，归为一个向量，并进行数字求和
I = zeros(Np,1);
for idx = 1:Npp = 2*pi*L./lamda0(idx);I0 = 2*lamda(:,idx)*cos(p);                % 省缺直流项I = I + I0;
end
figure
plot(t,I);title('平衡探测器下的干涉波形')
xlabel('时间周期');
ylabel('干涉强度');

可以修改"方差miu值"或者光程差L值，发现，（1）当miu或者L大于一定值的时候，干涉信号消失了。这就传说中光程差不能大于相干长度，否则不能探测到干涉信号。（2）L=0时fft后的镜面峰最大，随着L的增大，镜面峰的高度开始下降，下降到最大值的3dB位置即为相干长度，这时候就可以推算出瞬时线宽。

练习四：不用波长在同一光程差下的干涉强度 

光的干涉，基本都是同一光源的光分为两部分，然后分别传播一定距离后叠加，这时候两臂产生的距离差主要指空间距离，如果以空气为介质的话，那就是光程差了。如果光源输出的光不是宽带光，那在同一光程差下的干涉后强度分布可以表示如下：

S = 60e3;                            %光程差，60um
lamda = linspace(600,1000,2000);     %波长分布 100nm～200nm
p = 2*pi*S./lamda;                   %相位差
I0 = 1;                              %设两光路的强度均为1
I = I0 + I0 + 2*sqrt(I0*I0)*cos(p);  %叠加后的强度
figure;title('不用波长在同一光程差下的干涉强度');
plot(lamda,I);
xlabel('波长(nm)');
ylabel('干涉强度');

从上图可以看出，(1) 在相同光程差下，输出宽光谱的光时，干涉强度是呈现类正弦分布的。（2）分布规律不是标准的正弦波形，而是前快后慢的变化，该变化规律与波长的倒数有关。那是因为，叠加后的变化项是cos(p)，而p=2*pi*S./lamda.

练习五、SS-OCT的干涉信号仿真（波长随时间均匀分布），以样品中存在两个反射镜面为例。注：此时仍然认为单个波长的谱宽是无限窄的，即不同波长之间是不干扰的（实际上会有少量干扰，所以才叫低相干）

Np = 1376;
% 系统参数
lamda0 = 1060;              % 中心波长，nm
delata_lamda = 50;          % 半高宽(带宽的一半)，nm
lamda = linspace(lamda0-delata_lamda,lamda0+delata_lamda,Np); %波长分布
I = linspace(20,20,Np);     % 输出光强分布
winG = gausswin(Np);        % 高斯窗
I_Output = I.*winG';        % 理想条件下，输出光谱(强度)为高斯型
% plot(lamda,I_Output);title('理想的光谱光强输出')% 双层镜面干涉的仿真
L1 = 2000e3;                  % 设：参考臂和样品臂的距离差,镜面1
L2 = 5000e3;                  % 同上,镜面2
p1 = 2*pi*L1./lamda;          % 得：相位差与波长的关系,镜面1
p2 = 2*pi*L2./lamda;          % 同上,镜面2
I_ref = 1/2*I_Output;         % 参考臂（波长与强度关系）
I_samp1 = 1/4*I_Output;       % 样品臂,镜面1
I_samp2 = 1/5*I_Output;       % 样品臂,镜面2
I1 = I_ref+I_samp1+2*sqrt(I_ref.*I_samp1).*cos(p1);   %干涉强度
I2 = I_ref+I_samp2+2*sqrt(I_ref.*I_samp2).*cos(p2);   %干涉强度
I_Acq = I1 + I2;
plot(I_Acq),title('双层镜面反射的干涉信号');
wavDat = 2*sqrt(I_ref.*I_samp1).*cos(p1)+ 2*sqrt(I_ref.*I_samp2).*cos(p2); %平衡探测后的干涉强度
noise = rand(1,Np);           %加少量的噪声
wavDat = wavDat+noise;% 计算采样后的信号
winHn = hanning(Np);           % 汉宁窗
rawD = wavDat.*winHn'; % 加窗
ft = abs(fft(rawD));           % fft
Intensity = 20*log10(ft);      % 20log10的灰度映射
figure
plot(Intensity),title('双镜面反射的信号还原(波长随时间均匀分布)');
% TODO：
% (1)Np为采样点的个数，修改Np可以提高采样深度
% (2)从上述图像可以看出，随着深度的加深，信号展宽了，这是由于，采样是
%    根据波长随时间均匀分布采样的，实际要求应该是波数随时间均匀分布。

从上图可以看出,当输出波长是随时间均匀分布时，随着样品深度的加深，镜面峰宽度越来越大，也就是分辨率下降。

练习六、SS-OCT 干涉仿真二(波数随时间均匀分布) 
Np = 1376;
% 系统参数
lamdaSt = 1000;             % 起始波长，nm
lamdaEd = 1100;             % 结束波长，nm
k = linspace(2*pi/lamdaSt,2*pi/lamdaEd,Np); % 波数均匀分布
lamda = 2*pi./k;            % 由波数得波长分布
I = linspace(20,20,Np);     % 输出光强分布
winG = gausswin(Np);        % 高斯窗
I_Output = I.*winG';        % 理想条件下，输出光谱(强度)为高斯型
% plot(lamda,I_Output);title('理想的光谱光强输出')% 双层镜面干涉的仿真
L1 = 2000e3;                  % 设：参考臂和样品臂的距离差,镜面1
L2 = 5000e3;                  % 同上,镜面2
p1 = 2*pi*L1./lamda;          % 得：相位差与波长的关系,镜面1
p2 = 2*pi*L2./lamda;          % 同上,镜面2
I_ref = 1/2*I_Output;         % 参考臂（波长与强度关系）
I_samp1 = 1/4*I_Output;       % 样品臂,镜面1
I_samp2 = 1/5*I_Output;       % 样品臂,镜面2
I1 = I_ref+I_samp1+2*sqrt(I_ref.*I_samp1).*cos(p1);   %干涉强度
I2 = I_ref+I_samp2+2*sqrt(I_ref.*I_samp2).*cos(p2);   %干涉强度
I_Acq = I1 + I2;
plot(I_Acq),title('双层镜面反射的干涉信号');
wavDat = 2*sqrt(I_ref.*I_samp1).*cos(p1)+ 2*sqrt(I_ref.*I_samp2).*cos(p2); %平衡探测后的干涉强度
noise = rand(1,Np);           %加噪声
winHn = hanning(Np);          %加汉宁窗% 计算采样后的信号
rawD = (wavDat+noise).*winHn'; % 加窗
ft = abs(fft(rawD));           % fft
Intensity = 20*log10(ft);      % log10的灰度映射
figure
plot(Intensity),title('双镜面反射的信号还原(波数随时间均匀分布)');

当波数随时间均匀分布时，即波长的倒数分布是均匀的，镜面峰不会因为深度的加大导致展宽。根据练习四的结论，是因为cos(p)产生的干涉线是标准的正弦了，fft后频率比较单一。

练习七，扫频源OCT系统干涉信号的采集，从上述可以看出，按照波数均匀输出的信号不会出现镜面峰展宽，但是扫频源系统的波长一般都是按照时间均匀分布的。那么可以从采集入手，下面开始仿真采集卡不同采集模式下对信号的影响。

Nr = 10000;                       %探测器响应速度1GHz，扫频周期100kHz
FWHM = 50;                        %半高宽
lamda_center = 1050;              %中心波长
L = 1000e3;                       %光程差，nm
miu = 0.2;                        %扫频源中瞬时谱宽影响因子
lamda_St = lamda_center-FWHM;
lamda_Ed = lamda_center+FWHM;%% 构造探测器输出的干涉信号
winG = gausswin(Nr);              %高斯窗
t = linspace(0,1,Nr);             %一个扫频周期
lamda0 = linspace(lamda_St,lamda_Ed,Nr);  %中心波长的分布%模拟光谱输出，强度（波长，时间）
lamda = zeros(Nr,Nr);
for idx = 1:NrlamdaT = 1/(sqrt(2*pi)*miu) * exp(-(lamda0-lamda0(idx)).^2/miu^2);lamda(:,idx) = winG'.*lamdaT;
end% 每一个时刻都有一系列波长在叠加干涉，也可以说一个波长在不同时刻产生叠加
I_SS = zeros(Nr,1);
for idx = 1:Nrp = 2*pi*L./lamda0(idx);I0_SS = 2*lamda(:,idx)*cos(p);               % 省缺直流项I_SS = I_SS + I0_SS;
end
figure
plot(t,I_SS);title('平衡探测器输出的干涉信号波形')
xlabel('周期T')
ylabel('强度')

下面将用采集卡对探测器输出的信号进行采集。 

Np = 2048;
lamda_t = linspace(lamda_St,lamda_Ed,Np);       %中心波长的分布
k = linspace(2*pi/lamda_St,2*pi/lamda_Ed,Np);   % 波数均匀分布
lamda_k = 2*pi./k;
t = linspace(0,1,Np);figure,title('采样波长随时间的变化')
hold on
plot(t,lamda_t-lamda_t(1),'r')
plot(t,lamda_k-lamda_k(1),'b')
hold offLoc_t = (lamda_t-lamda_St)/(2*FWHM)*Nr;
Loc_k = (lamda_k-lamda_St)/(2*FWHM)*Nr;
Loc_t = round(Loc_t);
Loc_k = round(Loc_k);
Loc_t(1) = 1;
Loc_k(1) = 1;
Wave_t = I_SS(Loc_t);
Wave_k = I_SS(Loc_k);figure,title('采集卡采集的波形信号')
hold on
plot(t,Wave_t,'r');
plot(t,Wave_k,'k');
legend('时间均匀采集','波数均匀采集');
hold off


从采样波形可以看出，即使是同一个探测器输出的波形，采集卡用时间均匀采集和波数均匀采集的信号是略有不同的。下面对采集的信号进行处理对比。

noise = rand(1,Np);                     %加噪声
winHn = hanning(Np);                    %加汉宁窗
% 时间均匀
rawD_t = (Wave_t+noise').*winHn;       % 加窗
Amplft_t = abs(fft(rawD_t));           % fft
Intensity_t = 20*log10(Amplft_t);      % log10的灰度映射
% 波数均匀
rawD_k = (Wave_k+noise').*winHn;       % 加窗
Amplft_k = abs(fft(rawD_k));           % fft
Intensity_k = 20*log10(Amplft_k);      % log10的灰度映射figure,title('波长输出随时间的变化')
subplot(2,1,1),plot(Intensity_t+10,'r');title('时间均匀采集')
subplot(2,1,2),plot(Intensity_k-10,'k');title('波数均匀采集')

明显可得，波数均匀采集的信号，直接fft后镜面峰得到的分辨率要高。

练习八，SS-OCT 干涉仿真三(波数随时间均匀分布，但是带材料色散)

参考臂或者样品臂的光在各自的传输的过程中，除了光纤，还会遇到透镜，分束镜等材料，这些材料一般都是存在色散的，即折射率随波长改变。光程是空间距离与折射率的乘积，不同波长传输的空间距离肯定是相等的，但是折射率就不见得了。如果两臂的材料长度不相同，仿真信号的变化如下：

%   Sellmeier公式：n^2-1 = B1*lamda^2/(lamda^2-C1) 
%                        + B2*lamda^2/(lamda^2-C2)
%                        + B3*lamda^2/(lamda^2-C3);
Np = 1376;
% 系统参数
lamdaSt = 1000;             % 起始波长，nm
lamdaEd = 1100;             % 结束波长，nm
k = linspace(2*pi/lamdaSt,2*pi/lamdaEd,Np); % 波数均匀分布
lamda = 2*pi./k;            % 由波数得波长分布
I = linspace(20,20,Np);     % 输出光强分布
winG = gausswin(Np);        % 高斯窗
I_Output = I.*winG';        % 理想条件下，输出光谱(强度)为高斯型% 计算色散分布，以材料 N-SF6HT为例     
lamda_disp = lamda/1000;    %换算为单位 um
B1=1.77931763;B2=0.338149866;B3=2.08734474;%色散常数，查表(以um为单位）
C1=0.0133714182;C2=0.0617533621;C3=174.01759;
N2_SF6HT = B1*(lamda_disp.^2)./(lamda_disp.^2-C1)...+ B2*(lamda_disp.^2)./(lamda_disp.^2-C2)...+ B3*(lamda_disp.^2)./(lamda_disp.^2-C3)...+ 1;
N_SF6HT = sqrt(N2_SF6HT);
L = 40e6;                       %材料长度 10mm
p_disp = 2*pi*L*N_SF6HT./lamda; %附加相位
p_disp = p_disp - p_disp(Np/2); %多的配长，调试时会被抵消掉
plot(lamda,p_disp);% 双层镜面干涉的仿真
L1 = 1000e3;                  % 设：参考臂和样品臂的距离差,镜面1
L2 = 2000e3;                  % 同上,镜面2
p1 = 2*pi*L1./lamda+p_disp;   % 得：相位差与波长的关系,镜面1(加入材料色散)
p2 = 2*pi*L2./lamda+p_disp;   % 同上,镜面2
I_ref = 1/2*I_Output;         % 参考臂（波长与强度关系）
I_samp1 = 1/4*I_Output;       % 样品臂,镜面1
I_samp2 = 1/5*I_Output;       % 样品臂,镜面2
wavDat = 2*sqrt(I_ref.*I_samp1).*cos(p1)+ 2*sqrt(I_ref.*I_samp2).*cos(p2); %平衡探测后的干涉强度
noise = rand(1,Np);           %加噪声
winHn = hanning(Np);          %加汉宁窗% 计算采样后的信号
rawD = (wavDat+noise).*winHn'; % 加窗
ft = abs(fft(rawD));           % fft
Intensity = 20*log10(ft);      % log10的灰度映射
plot(Intensity),title('双镜面反射的信号还原(带材料色散)');

可以看出，即便是波数随时间均匀分布，但是由于材料色散的存在，镜面峰也出现展宽（该展宽与深度无关）。由于不是低相干仿真，所以展宽变化不是很明显。

这篇关于扫频源光学相干层析（SS-OCT）原理仿真（附matlab代码）（未完待更新）的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/965473。

相关文章

C++使用栈实现括号匹配的代码详解

C++使用栈实现括号匹配的代码详解

《C++使用栈实现括号匹配的代码详解》在编程中,括号匹配是一个常见问题,尤其是在处理数学表达式、编译器解析等任务时,栈是一种非常适合处理此类问题的数据结构,能够精确地管理括号的匹配问题,本文将通过C+... 目录引言问题描述代码讲解代码解析栈的状态表示测试总结引言在编程中，括号匹配是一个常见问题，尤其是在
阅读更多...
Java调用DeepSeek API的最佳实践及详细代码示例

Java调用DeepSeek API的最佳实践及详细代码示例

《Java调用DeepSeekAPI的最佳实践及详细代码示例》：本文主要介绍如何使用Java调用DeepSeekAPI,包括获取API密钥、添加HTTP客户端依赖、创建HTTP请求、处理响应、... 目录1. 获取API密钥2. 添加HTTP客户端依赖3. 创建HTTP请求4. 处理响应5. 错误处理6.
阅读更多...
使用 sql-research-assistant进行 SQL 数据库研究的实战指南(代码实现演示)

使用 sql-research-assistant进行 SQL 数据库研究的实战指南(代码实现演示)

《使用sql-research-assistant进行SQL数据库研究的实战指南(代码实现演示)》本文介绍了sql-research-assistant工具,该工具基于LangChain框架,集... 目录技术背景介绍核心原理解析代码实现演示安装和配置项目集成LangSmith 配置（可选）启动服务应用场景
阅读更多...
Python中顺序结构和循环结构示例代码

Python中顺序结构和循环结构示例代码

《Python中顺序结构和循环结构示例代码》：本文主要介绍Python中的条件语句和循环语句,条件语句用于根据条件执行不同的代码块,循环语句用于重复执行一段代码,文章还详细说明了range函数的使... 目录一、条件语句（1）条件语句的定义（2）条件语句的语法（a）单分支 if（b）双分支 if-else（
阅读更多...
MySQL数据库函数之JSON_EXTRACT示例代码

MySQL数据库函数之JSON_EXTRACT示例代码

《MySQL数据库函数之JSON_EXTRACT示例代码》：本文主要介绍MySQL数据库函数之JSON_EXTRACT的相关资料,JSON_EXTRACT()函数用于从JSON文档中提取值,支持对... 目录前言基本语法路径表达式示例示例 1: 提取简单值示例 2: 提取嵌套值示例 3: 提取数组中的值注意
阅读更多...
CSS3中使用flex和grid实现等高元素布局的示例代码

CSS3中使用flex和grid实现等高元素布局的示例代码

《CSS3中使用flex和grid实现等高元素布局的示例代码》：本文主要介绍了使用CSS3中的Flexbox和Grid布局实现等高元素布局的方法，通过简单的两列实现、每行放置3列以及全部代码的展示，展示了这两种布局方式的实现细节和效果，详细内容请阅读本文，希望能对你有所帮助... 过往的实现方法是使用浮动加
阅读更多...
JAVA调用Deepseek的api完成基本对话简单代码示例

JAVA调用Deepseek的api完成基本对话简单代码示例

《JAVA调用Deepseek的api完成基本对话简单代码示例》：本文主要介绍JAVA调用Deepseek的api完成基本对话的相关资料,文中详细讲解了如何获取DeepSeekAPI密钥、添加H... 获取API密钥首先，从DeepSeek平台获取API密钥，用于身份验证。添加HTTP客户端依赖使用Jav
阅读更多...
Java实现状态模式的示例代码

Java实现状态模式的示例代码

《Java实现状态模式的示例代码》状态模式是一种行为型设计模式,允许对象根据其内部状态改变行为,本文主要介绍了Java实现状态模式的示例代码,文中通过示例代码介绍的非常详细,需要的朋友们下面随着小编来... 目录一、简介1、定义2、状态模式的结构二、Java实现案例1、电灯开关状态案例2、番茄工作法状态案例
阅读更多...
nginx-rtmp-module模块实现视频点播的示例代码

nginx-rtmp-module模块实现视频点播的示例代码

《nginx-rtmp-module模块实现视频点播的示例代码》本文主要介绍了nginx-rtmp-module模块实现视频点播,文中通过示例代码介绍的非常详细,对大家的学习或者工作具有一定的参考学习... 目录预置条件Nginx点播基本配置点播远程文件指定多个播放位置参考预置条件配置点播服务器 192.
阅读更多...
MySQL中的MVCC底层原理解读

MySQL中的MVCC底层原理解读

《MySQL中的MVCC底层原理解读》本文详细介绍了MySQL中的多版本并发控制（MVCC）机制,包括版本链、ReadView以及在不同事务隔离级别下MVCC的工作原理,通过一个具体的示例演示了在可重... 目录简介ReadView版本链演示过程总结简介MVCC（Multi-Version Concurr
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2