扫频源光学相干层析（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。

相关文章

深入探索协同过滤：从原理到推荐模块案例

深入探索协同过滤：从原理到推荐模块案例

文章目录 前言一、协同过滤1. 基于用户的协同过滤（UserCF）2. 基于物品的协同过滤（ItemCF）3. 相似度计算方法 二、相似度计算方法1. 欧氏距离2. 皮尔逊相关系数3. 杰卡德相似系数4. 余弦相似度 三、推荐模块案例1.基于文章的协同过滤推荐功能2.基于用户的协同过滤推荐功能 前言     在信息过载的时代，推荐系统成为连接用户与内容的桥梁。本文聚焦于
阅读更多...
poj3468(线段树成段更新模板题)

poj3468(线段树成段更新模板题)

题意：包括两个操作：1、将[a.b]上的数字加上v；2、查询区间[a,b]上的和 下面的介绍是下解题思路： 首先介绍  lazy－tag思想：用一个变量记录每一个线段树节点的变化值，当这部分线段的一致性被破坏我们就将这个变化值传递给子区间，大大增加了线段树的效率。 比如现在需要对[a,b]区间值进行加c操作，那么就从根节点[1,n]开始调用update函数进行操作，如果刚好执行到一个子节点，
阅读更多...
hdu1394(线段树点更新的应用)

hdu1394(线段树点更新的应用)

题意：求一个序列经过一定的操作得到的序列的最小逆序数 这题会用到逆序数的一个性质，在0到n-1这些数字组成的乱序排列，将第一个数字A移到最后一位，得到的逆序数为res-a+(n-a-1) 知道上面的知识点后，可以用暴力来解 代码如下： #include<iostream>#include<algorithm>#include<cstring>#include<stack>#in
阅读更多...
hdu1689(线段树成段更新)

hdu1689(线段树成段更新)

两种操作：1、set区间[a,b]上数字为v；2、查询[ 1 , n ]上的sum 代码如下： #include<iostream>#include<algorithm>#include<cstring>#include<stack>#include<queue>#include<set>#include<map>#include<stdio.h>#include<stdl
阅读更多...
hdu4407(容斥原理)

hdu4407(容斥原理)

题意：给一串数字1,2，......n，两个操作:1、修改第k个数字，2、查询区间[l,r]中与n互质的数之和。 解题思路：咱一看，像线段树，但是如果用线段树做，那么每个区间一定要记录所有的素因子，这样会超内存。然后我就做不来了。后来看了题解，原来是用容斥原理来做的。还记得这道题目吗？求区间[1,r]中与p互质的数的个数，如果不会的话就先去做那题吧。现在这题是求区间[l,r]中与n互质的数的和
阅读更多...
活用c4d官方开发文档查询代码

活用c4d官方开发文档查询代码

当你问AI助手比如豆包，如何用python禁止掉xpresso标签时候，它会提示到 这时候要用到两个东西。https://developers.maxon.net/论坛搜索和开发文档 比如这里我就在官方找到正确的id描述 然后我就把参数标签换过来
阅读更多...
poj 1258 Agri-Net（最小生成树模板代码）

poj 1258 Agri-Net（最小生成树模板代码）

感觉用这题来当模板更适合。 题意就是给你邻接矩阵求最小生成树啦。~ prim代码：效率很高。172k...0ms。 #include<stdio.h>#include<algorithm>using namespace std;const int MaxN = 101;const int INF = 0x3f3f3f3f;int g[MaxN][MaxN];int n
阅读更多...
hdu 1754 I Hate It（线段树，单点更新，区间最值）

hdu 1754 I Hate It（线段树，单点更新，区间最值）

题意是求一个线段中的最大数。 线段树的模板题，试用了一下交大的模板。效率有点略低。 代码： #include <stdio.h>#include <string.h>#define TREE_SIZE (1 << (20))//const int TREE_SIZE = 200000 + 10;int max(int a, int b){return a > b ? a :
阅读更多...
AI行业应用(不定期更新)

AI行业应用(不定期更新)

ChatPDF 可以让你上传一个 PDF 文件，然后针对这个 PDF 进行小结和提问。你可以把各种各样你要研究的分析报告交给它，快速获取到想要知道的信息。https://www.chatpdf.com/
阅读更多...
GIS图形库更新2024.8.4-9.9

GIS图形库更新2024.8.4-9.9

更多精彩内容请访问 dt.sim3d.cn ，关注公众号【sky的数孪技术】，技术交流、源码下载请添加微信：digital_twin123 Cesium 本期发布了1.121 版本。重大新闻，Cesium被Bentley收购。 ✨ 功能和改进 默认启用 MSAA，采样 4 次。若要关闭 MSAA，则可以设置scene.msaaSamples = 1。但是通过比较，发现并没有多大改善。
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2