扫频源光学相干层析(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:Np
    if(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:Np
    lamdaT = 1/(sqrt(2*pi)*miu) * exp(-(lamda0-lamda0(idx)).^2/miu^2);
    lamda(:,idx) = winG'.*lamdaT;
end

figure
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:Np
    I0 = lamda(:,idx);
    I = I + I0; 
end
figure
plot(t,I);title('验证模拟的输出谱是否为高斯型')

% 每一个时刻都有一系列波长在叠加干涉,也可以说一个波长在不同时刻产生叠加
I = zeros(Np,1);
for idx = 1:Np
    p = 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:Np
    p = 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:Nr
    lamdaT = 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:Nr
    p = 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 off

Loc_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('双镜面反射的信号还原(带材料色散)');

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

雷达系统设计matlab仿真代码具体如下: ```matlab % 创建雷达信号 tx_signal = randn(1, 1000); % 随机生成1000个点的发送信号 % 设置雷达参数 fc = 10e9; % 雷达中心频率为10 GHz fs = 100e6; % 采样频率为100 MHz T_sweep = 100e-6; % 雷达扫频时间为100 us bw = 1e6; % 扫频带宽为1 MHz nm = 2; % 设置噪声系数 % 生成雷达接收信号 rx_signal = radar_receiver(tx_signal, fc, fs, T_sweep, bw, nm); % 处理雷达接收信号 processed_signal = radar_processing(rx_signal); % 雷达接收信号处理函数 function processed_signal = radar_processing(rx_signal) % 雷达信号处理算法 % 对接收信号进行去噪、调频去斜、目标检测等处理 % 这里只做简单的处理示例,去除直流分量和方波滤波 % 去除直流分量 rx_signal = rx_signal - mean(rx_signal); % 方波滤波 b = [1 -1]; a = 1; filtered_signal = filter(b, a, rx_signal); processed_signal = filtered_signal; end % 雷达接收信号生成函数 function rx_signal = radar_receiver(tx_signal, fc, fs, T_sweep, bw, nm) % 雷达接收信号生成算法 % 对发射信号进行正弦频率调制、加性高斯白噪声等处理 % 这里只做简单的处理示例,加入高斯白噪声 % 正弦频率调制 t = (0:length(tx_signal)-1) / fs; rx_signal = tx_signal .* exp(1j * 2 * pi * fc * t); % 加入高斯白噪声 rx_signal = rx_signal + nm * randn(size(rx_signal)); end ``` 以上是一个简单的雷达系统设计matlab仿真代码,包括雷达接收信号生成函数和雷达接收信号处理函数。在生成接收信号时,对发射信号进行了正弦频率调制,并加入了高斯白噪声。在处理接收信号时,去除了直流分量并进行了方波滤波。具体的算法实现可以根据具体需求进行修改和扩展。
评论 12
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值