光学相干层析成像OCT的介绍与模拟（附Python和Matlab代码）

wolf

写在前面的废话

本人正在学习OCT成像系统的有关知识，会陆陆续续发表一些总结在这个专栏中，希望借助这个机会对自己学到的东西做一个总结。这篇文章主要参考了汪立宏教授的神作《Biomedical Optics: principles and imaging》，Fujimoto的《Optical Coherence Tomography ：Technology and Applications》以及几篇综述文章。Matlab代码为汪教授书中源码，我自己按照思路用Python写了一遍。另外，我会在引用块说一些题外话，觉得罗嗦可以忽略不看哈^_^。欢迎大家一起讨论，如果有不对的地方，希望在评论区告知我，十分感谢！
<hr/>背景介绍

OCT是一种具有非侵入性、高分辨率特性的三维成像系统，其在医学领域特别是眼科发挥着极其重要的作用。第一台时域OCT（Time-domain OCT，TD-OCT）由华人David Huang以及他的导师James G. Fujimoto于1991年发明。1997年S. R. Chinn和其导师Jame G. Fujimoto发明了第一台扫频OCT(Optical coherence tomography using a frequency-tunable optical source)。2002年Maciej Wojtkowski和其导师Adolf F. Fercher发明了第一台傅立叶域OCT(Fourier-domain OCT，FD-OCT)。之后OCT的发展主要是围绕傅立叶域OCT展开。
瞻仰一下OCT领域的大佬们



2017 Russ Prize,颁发给在OCT领域做出突出贡献的科学家们

OCT主要分为时域OCT(TD-OCT)和傅立叶域OCT(FD-OCT)，傅立叶域OCT又分为光谱域OCT（Spectral-domain OCT, SD-OCT）和扫频OCT（Swept Source OCT, SS-OCT）。



OCT分类

相比其他的光学成像系统，OCT具有自己独特的优势。从下图可见，OCT弥补了从共聚焦显微镜到超声成像之间的一个空缺，其分辨率在1-15μm，穿透深度为1-2mm。



OCT与其他成像仪器的对比

<hr/>迈克尔逊干涉仪

OCT是基于迈克尔逊干涉仪发明的，其光路与迈克尔逊干涉仪基本一致，区别是将迈克尔逊干涉仪的M2平面镜换成了生物组织等其他被观察对象。所以我们先来介绍一下Michelson interferometer。

一说到这就感慨万分。伟大发明离我们很近啊有木有？大二做基础物理实验时就做过迈克尔逊干涉仪，我们只需要把一块平面镜换成生物组织，然后移动另一面平面镜，在观察点记录光强变化，就能知道生物组织不同深度的信息，以此成像！第一代OCT就这么发明了哈^_^虽然有点想当然，但想法真的很简单啊！

迈克尔逊干涉仪

由图可知，单色光从光源射出后，被分束器分成了两个光束，一个射向平面镜然后反射，在OCT中我们把它叫做Reference arm；另一束穿过分束器射向样本(在迈克尔逊干涉仪中为平面镜)然后反射,叫做Sample arm。两束光经反射后发生干涉，然后射向探测器。

值得注意的是，射向Reference平面镜的光束三次穿过了分束器，而射向Sample的光只穿过一次分束器。

所以在真正实验中，会在Sample的光路上加一个光路补偿器。OCT系统中也有类似的补偿，例如在Sample arm中有透镜，那么在Reference arm中就会加光学器件进行补偿。

经过反射回到分束器处的两束光的电场可以描述为下面的表达式：
\begin{equation}\begin{array}{l} E_{R}=E_{R 0} \exp \left(i\left(2 k_{R} l_{R}-\omega t\right)\right) \\ E_{S}=E_{S 0} \exp \left(i\left(2 k_{S} l_{S}-\omega t\right)\right) \end{array}\\\end{equation}
E为振幅，k为传播系数，l为臂长（从分束器分光点开始），数字2表示一来一回这样一个传播过程。叠加后的电场如下：
I(t)=\underline{\underline{E_{R 0}^{2}+E_{S 0}^{2}}}+\underline{\underline{2 E_{R 0} E_{S 0} \cos \left(2 k_{S} l_{S}-2 k_{R} l_{R}\right)}}\\ DC\ term \hspace{ 2.8cm }AC\ term

探测器探测到的是一个反应周期（皮秒或者纳秒）内的平均值: \begin{equation}I(t)=\left\langle\left|E_{R}+E_{S}\right|^{2}\right\rangle\\ \end{equation}

它由直流电(Direct Current, DC)和交流电(Alternating Current, AC)两部分组成，若令：
k_{R}=k_{S}=k=\frac{2 \pi n }{\lambda_{0}}\\ \Delta \phi=2 k\left(l_{S}-l_{R}\right)=2 \pi \frac{2 n \Delta l}{\lambda_{0}}\\ I(t)=I_{const}+2E_{R0}E_{S0}cos\left(\Delta\phi\right)  
\Delta \phi 随着 \begin{equation}\Delta l\end{equation} 变化便产生了圆形干涉条纹。


当移动Reference平面镜时， \begin{equation}\Delta l\end{equation} 也会发生改变，发生吞吐现象。这就是迈克尔逊干涉仪的原理。
<hr/>相干长度和相干时间

OCT是根据低相干光成像，所以有必要知道光源对相干性的影响。

例如双缝干涉，如果光源不是单色光，而是有一定频谱宽度，这时候不同波长的光产生的干涉条纹叠加在一起，随着条纹级数 K 的增大，干涉条纹越来越不清晰，只有当K\le\frac{\lambda_{0}^{2}}{\Delta \lambda}时我们才可以清晰的分辨中心波长产生的干涉条纹。

相干时间\tau_{c}定义为光的电场 E(t) 的自相干函数 G(\tau) 的半峰全宽FWHM。只有两束相干光时间差小于相干时间\tau_{c}时才能产生可以观测的相干信号。而相干长度也类似，相干长度
G(\tau)=\int_{-\infty}^{+\infty} E(t) E(t+\tau) d t\\
l_{c}=c \tau_{c} ， c 为光速。只有两束相干光之间的光程差小于相干长度时，才可以产生可观测的相干信号。根据 Wiener–Khintchine 定理，光源的频谱函数 S(w) 和光的电场的自相干函数G(\tau)满足下面关系：（其中 E(w) 为E(t)的傅立叶变换）
S(\omega)=|E(\omega)|^{2}=\int_{-\infty}^{+\infty} G(\tau) \exp (i \omega \tau) d \tau\\
如果 S(w) 为高斯型，那么我们可以求出相干长度 l_c :
l_{c}=\frac{4 \ln 2}{\pi} \frac{\lambda_{0}^{2}}{\Delta \lambda}\\



高斯光

其中 \lambda_0 为中心波长， \Delta\lambda 为FWHM。所以，中心波长越大，FWHM越小，相干长度越大。所以在OCT中，只有Reference arm反射的光和Sample arm反射的光之间的光程差小于相干长度l_c时，我们才可以观察到相干信号。那么我们是希望相干长度 l_c 越大越好呢，还是越小越好呢？答案是尽量小一些，原因是我们只希望获取观察区域的信息，如果相干长度太长，受其他地方的干扰，我们获取的相干信号就不能够精准的对应观察区域了。这里的分析其实也暗示着，OCT的纵向分辨率会与相干长度有关。
<hr/>TD-OCT的介绍与模拟

OCT光路可以基于光直接在空间中传播，也可以使用光纤，现在一般实验室常用光纤。但使用光纤时需要使用单模光纤，因为多模光纤中光的色散使不同频率的光射出光纤时间不一样，最后会导致OCT系统纵向分辨率降低。TD-OCT中，光源使用的是超发光二极管（Superluminescent diode, SLD），通常实验室里使用近红外光，原因是其波长长，不容易被散射，反射效果更好。



TD-OCT

Reference arm处还是使用平面镜进行反光，但注意这里加了一个振荡器，目的我们后面再谈；Sample arm处则放一块生物组织样本。TD-OCT的探测器为光电二极管，将光信号转换为电信号，然后进行信号处理。
TD-OCT中，光源有一定的频率宽带，每一种频率的光都会在分束成Sample arm和Reference arm之后，反射回来叠加产生一个相干信号：
I(\omega)=\left|E_{R}(\omega)+E_{S}(\omega)\right|^{2}=\underline{\underline{\left|E_{R}(\omega)\right|^{2}+\left|E_{S}(\omega)\right|^{2}}}+\underline{\underline{2 \operatorname{Re}\left\{E_{R}(\omega) E_{S}^{*}(\omega)\right\}}}\\  \hspace{4.5cm}DC \ term \hspace{2cm}AC\ term\\
有用的信号为交流电AC部分。然后不同频率的光经过这一过程产生的相干信号进行叠加：
I_{\mathrm{AC}}=2 \operatorname{Re}\left\{\int_{-\infty}^{\infty} E_{R}(\omega) E_{S}^{*}(\omega) d \omega\right\}\\
Sample arm处相当于有无数个反射面，每个反射面都会造成光的延迟，与Reference arm产生光程差。由上面提到的相干长度部分的知识可以知道，如果保持Reference arm和Sample arm静止不动，那么只有光程差小于相干长度 l_c 的层反射后产生的叠加干涉信号很明显，而远处的信号在叠加后几乎为0。所以理论上，我们在探测器处可以接收到一个如下图一样的信号：


中间部分信号强度最大，原因是无论入射光频率是多少（不考虑折射率），其光程差为0。但实际上我们的探测器由于灵敏度不够，接收到的信号是一个较长时间内的平均信号，也就是光电二极管只感受到了一个恒定的光强，上图这些波动或者这个信号的包络是观察不到的。那如何探测信号的包络呢？还记得上面说的Reference arm上的振荡器吗？它的作用就是通过周期性振动来产生一个Doppler shift,由于Reference arm移动也可以改变光程差，周期性移动就会带来探测器处光强的变化，经过信号处理后将信号复原，我们就得到了Sample上一个点附近信号的包络。然后再移动Reference arm平面镜的位置，进行震荡，得到第二个点的信号包络，以此类推，最后可以得到一个二维或者三维的图。
如果我们光源使用\lambda_0=830nm,\ \Delta\lambda=60nm 的高斯光，经过公式推导(由知友提问，遂将详细过程写出来，这是我自己的推导过程，和书上有所不一样，书上用了群速相速等知识）得到交流AC部分的公式：

I_{\mathrm{AC}}=2 \operatorname{Re}\left\{\int_{-\infty}^{\infty} E_{R}(\omega) E_{S}^{*}(\omega) d \omega\right\}\\ \because \ S(\omega) \propto E_{R 0}(\omega) E_{S 0}^{*}(\omega)\\ I_{\mathrm{AC}} \propto \operatorname{Re}\left\{\int_{-\infty}^{\infty} S(\omega) \exp (-i \Delta \phi(\omega)) d \omega\right\}\\ I_{\mathrm{AC}} \propto \operatorname{Re}\left\{\int_{-\infty}^{\infty} S(\omega) \exp (-i\frac{2w\Delta l}{c}) d \omega\right\}\\   I_{\mathrm{AC}} \propto \operatorname{Re}\left\{\frac{1}{\sqrt{2 \pi} \sigma_{\omega}} \int_{-\infty}^{\infty} \exp \left(-\frac{\left(\omega-\omega_{0}\right)^{2}}{2 \sigma_{\omega}^{2}}\right)\dot \ \exp (-i\frac{2w\Delta l}{c}) d \omega\right\}\\   I_{\mathrm{AC}} \propto \operatorname{Re}\left\{\frac{1}{\sqrt{2 \pi} \sigma_{\omega}} \int_{-\infty}^{\infty} \exp \left(-\left(\frac{\left(\omega-\omega_{0}\right)^{2}}{2 \sigma_{\omega}^{2}} +i\frac{2\Delta l(w-w_0)}{c}-i\frac{2w_0\Delta l}{c}\right)\right) d \omega\right\}\\     
套用如下计算公式：
\int_{-\infty}^{\infty} \exp \left(-\left(a x^{2}+b x+c\right)\right) d x=\sqrt{\frac{\pi}{a}} \exp \left(\frac{b^{2}-4 a c}{4 a}\right)\\
我们得到：
I_{\mathrm{AC}} \propto exp\left(\frac{2\ \sigma_w^2 \ \Delta l^2 }{c^2}\right)cos\frac{2w_0\Delta l}{c}\\
高斯光源标准差 \sigma_w 以及相干长度存 l_c 在如下关系：
\sigma_w = \frac{\pi c\Delta \lambda}{\lambda_0^2 \sqrt{2 \ln 2}}\\ l_{c}=\frac{4 \ln 2}{\pi} \frac{\lambda_{0}^{2}}{\Delta \lambda}\\ \therefore \ \frac{\sigma_w}{c}=\frac{2\sqrt{2ln2}}{l_c}\\ 另外: k=\frac{w}{c}
代入可得：
\\I_{\mathrm{AC}}=\exp \left(-16 \ln 2\left(\frac{\Delta l}{l_{\mathrm{c}}}\right)^{2}\right) \cos \left(2 k_{0} \Delta l\right)\\
最后模拟获取包络的过程如下：


另外，还可以推导出OCT系统中的分辨率公式。
纵向分辨率：\Delta z=\frac{l_{c}}{2}=\frac{1}{2} t_{c} c=\frac{2 \ln 2}{\pi } \frac{\lambda_{0}^{2}}{\Delta \lambda}
横向分辨率： \Delta x=\frac{2 }{\pi} \frac{\lambda_{0}}{\mathrm{NA}}
所以，要想得到较好的纵向分辨率，可以使用频带较宽的光源，以及尽量小的中心波长。提高横向分辨率的方法也是使用较小的中心波长以及较大数值孔径的物镜。但是波长也不是越小越好，相反，实验室常常使用中心波长较长的近红外光，原因是波长较短时发生散射比较严重，并且无论波长多少，光吸收也是一个需要考虑的因素，最后综合这些因素后，选择近红外光为中心波长，以达到最大的成像深度。
<hr/>FD-OCT的介绍与模拟

TD-OCT最大的缺点就是Referene arm需要机械移动，会造成成像速度慢。而FD-OCT比较好的解决了这一个问题，FD-OCT不需要通过移动Reference arm来获得深度信息，而是能够一次扫描直接获取深度信息，所以成像速度和灵敏度都大大提高。



FD-OCT

对于任一频率 w ，Sample arm和Reference arm的反射信号如下：
\begin{aligned} &E_{R}(\omega)=E_{0}(\omega) r_{R} \exp \left(i\left(2 k_{R}(\omega) l_{R}-\omega t\right)\right)\\ &E_{S}(\omega)=E_{0}(\omega) \int_{-\infty}^{+\infty} r_{S}^{\prime}\left(l_{S}\right) \exp \left(i\left(2 k_{S}(\omega) l_{S}-\omega t\right)\right) d l_{S} \end{aligned}\\
r 为反射率，对于Sample arm,其相当于有无数个反射面，每个反射面拥有不同的反射率，然后进行累加。最后来自两臂的光进行叠加：
令：\frac{k_{R}}{n_{R}}=\frac{k_{S}}{n_{S}}=k=\frac{\omega}{c}\\
\begin{aligned} I(k)=\left|E_{R}(k c)+E_{S}(k c)\right|^{2}=& S(k) r_{R}^{2}\hspace{1cm}    (1)\\ &+2 S(k) r_{R} \int_{-\infty}^{+\infty} r_{S}^{\prime}\left(l_{S}\right) \cos \left(2 k\left(n_{S} l_{S}-l_{R}\right)\right) d l_{S}\ (2)\\ &+S(k)\left|\int_{-\infty}^{+\infty} r_{S}^{\prime}\left(l_{S}\right) \exp \left(i 2 k\left(n_{S} l_{S}\right)\right) d l_{S}\right|^{2}\ (3) \end{aligned}\\
（1）DC term 反映的是Reference平面镜的反射强度，它可以通过使反射率 r_{S}^{\prime}=0 来测量，具体的方法为用一块黑布将Sample arm的光路挡住，得到背景信息。
（2）Cross-correlation term 反映的是互相关信号，也是有效部分。每一层的反射率  r_{S}^{\prime}\left(l_{S}\right)  都被编码为有不同频率 2\left(n_{S} l_{S}-l_{R}\right) 的余弦函数的振幅，而且，频率随着深度 l_S 增加而变高，所以如果我们使用一个傅里叶变换，就可以得到不同深度的反射率，以此成像。

之所以这里频率为 2\left(n_{S} l_{S}-l_{R}\right) ，是因为函数自变量为 k ,而不是指光的频率。
同时，这里频率也会随着 l_R 增大而减小，所以在做实验时，有时也会调整一下Reference arm的平面镜的位置。

（3）Auto-correlation term 反映的是Sample中不同层反射的信号叠加，通常这个信号频率会很小（光程差小），而且强度弱（组织的反射率比平面镜小很多，平方之后更小）。要获取这一部分的信息也很简单，我们将Reference arm的光路用黑布挡住，使 r_R=0 。
我们最关心的还是（2）部分产生的信号，所以可以使用上述黑布挡光的方法得到并去除（1）（3）两部分的信息，也可以移动样品使相位产生 \pi 的变化得到另一个干涉信号，这时这个干涉信号的（2）部分与原来的正好符号相反，两个信号相减就去掉了（1）（3）项。
所以最后我们可以将干涉信号看做是一个以 k 为自变量的，频率为2\left(n_{S} l_{S}-l_{R}\right)的无穷级数的和。最后我们将干涉信号通过一个光谱仪进行色散（相当于一个傅立叶变换），将不同频率的光分开，用一个CCD接受信号，每个频率的光的振幅对应着Sample不同深度的反射率，这样就可以成像了。

Reference arm的镜面位置必须在Sample外面，也就是使得 2\left(n_{S} l_{S}-l_{R}\right) 恒为正数，不然如果镜面位置在Sample里面，那么会造成混乱，因为正频率和负频率经过色散后都是表现为正频率，但是他们各自振幅代表着不同深度的反射率。

光谱仪

在模拟中，光源依旧使用高斯光。设平面镜的镜面位置为原点，Sample全部在正轴方向，同时，我们只选取2个反射面作为示例，每个反射面有不同的反射系数。模拟了干涉信号，它的傅里叶变换，傅里叶变换的高斯滤波，以及去除掉（1）（3）两项且经过高斯滤波之后的信号。


对第二个图的解释如下：



Sample arm具有两个反射面

换在一个实际的例子上（眼睛）：


以上FD-OCT的光源使用的是宽带光源，所以又叫Spectra-domain OCT;另一种FD-OCT为Swept source OCT，其使用高速的扫频光源，每次只射出一种频率的光，如下：


它的算法和上面很相似，只不过探测器不是CCD，换成了光电二极管，也没有了光谱仪。 I(k) 由每次直接射的不同频率k的光的光强组成。然后将这一序列进行傅里叶变换即可（相当于SD-OCT中光谱仪的作用）。后面的处理就一样了。
SD-OCT的分辨率取决于中心波长以及FWHM；SS-OCT的纵向分辨率取决于光谱包络的FWHM和中心波长
<hr/>代码汇总

TD-OCT(Python)
&#34;&#34;&#34;
Created on Tue Mar 17 15:39:03 2020

@author: Fuwei
&#34;&#34;&#34;
import numpy as np
import matplotlib.pyplot as plt

pi = np.pi
ln2 = np.log(2)
lambda_0 = 830e-9 #center wavelength
d_lambda = 60e-9 #bandwidth (delta lambda)
c = 3e8 #speed of light
lc = 4*ln2/pi*lambda_0**2/d_lambda #coherence length
Number_of_periods = 0.5*lc/(lambda_0/2) #of periods in FWHM

plt.figure(1)
plt.suptitle(&#34;TD-OCT Simulation&#34;,fontsize=25,fontweight=800)
plt.subplot(4,1,1) #interferogram
N = 2**12 #number of sampling points
dl = lc*np.linspace(-2, 2,N) #array for Delta_1
k_0 = 2*pi/lambda_0 #propagation constant
I_ac = np.exp(-16*ln2*((dl/lc)**2))*np.cos(2*k_0*dl)
plt.plot(dl/lc, I_ac, &#39;k&#39;)
plt.title(&#34;(a) Interferogram&#34;)
plt.xlabel(r&#34;$\mathregular{\Delta{l}/l_c}$&#34;)
plt.ylabel(&#34;Signal&#34;)
plt.axis([-0.6, 0.6, -1, 1.2])

plt.subplot(4,1,2) #rectified interferogram
I_rec = np.abs(I_ac) #absolute value of I_ac

plt.plot(dl/lc, I_rec, &#39;k&#39;)
plt.title(&#34;(b) Rectified interferogram&#34;)
plt.xlabel(r&#34;$\mathregular{\Delta{l}/l_c}$&#34;)
plt.ylabel(&#34;Signal&#34;)
plt.axis([-0.6, 0.6, -1, 1.2])

plt.subplot(4,1,3) #spectrum of the rectified interferogram

Frec_1 = np.fft.fft(I_rec)/np.sqrt(N)# order of frequencies: 0,1...(N/2-1),-N/2,-(N/2-1)...-1
Frec_2 =np.fft.fftshift(Frec_1)#shifted order of frequencies: -N/2,-(N/2-1)...-1, 0,1...(N/2-1)
d_freq = 1/(4*lc)  #4=2-(-2) frequency bin size = 1/sampling range
freq=d_freq*np.arange(-N/2,N/2) #frequency array
plt.plot(freq*lambda_0, np.abs(Frec_2), &#39;k&#39;)
plt.title(&#34;(c) Spectrum of the rectified interferogram&#34;)
plt.xlabel(r&#34;Frequency $\mathregular{(1/\lambda_0)}$&#34;)
plt.ylabel(&#34;Amplitude&#34;)
plt.axis([-10,10,0,5.5])


plt.subplot(4,1,4) #envelope
freq_cut = 1/lambda_0/2 #cut-off frequency for filtering
i_cut = int(round(freq_cut/d_freq)) #convert freq_cut to an array index
F_filt = Frec_1.copy() #initialize array
F_filt[i_cut:N-i_cut+1]=0 #filter
I_filt = np.abs(np.fft.ifft(F_filt))*np.sqrt(N) #amplitude of inverse FFT
I_ac_en =  np.exp(-16*ln2*((dl/lc)**2))
plt.plot(dl/lc, I_filt/I_filt.max(), &#39;k&#39;,label=&#34;Demodulted&#34;)
plt.plot(dl[np.arange(0,N,int(N/32))]/lc, I_ac_en[np.arange(0,N,int(N/32))], &#39;ro&#39;,label=&#34;original&#34;)
plt.title(&#34;( d ) Envelopes&#34;)
plt.xlabel(r&#34;$\mathregular{\Delta{l}/l_c}$&#34;)
plt.ylabel(&#34;Signals&#34;)
plt.axis([-0.6,0.6,-1,1.2])
plt.legend()

plt.subplots_adjust(hspace=0.6)

plt.show()
# plt.tight_layout(h_pad=2)TD-OCT(Matlab)
% Use SI units throughout
lambda0 = 830E-9; % center wavelength
dlambda = 60E-9; % bandwidth (delta lambda)
c = 3E8; % speed of light
lc=4*log(2)/pi*lambda0^2/dlambda % coherence length
Number_of_periods = 0.5*lc/(lambda0/2) % # of periods in FWHM
figure(1);
N = 2^12; % number of sampling points
dl = lc*linspace(-2,2, N); % array for Delta_l
k0 = 2*pi/lambda0; % propagation constant
subplot(4, 1 , 1 ) % interferogram
Iac = exp(-16*log(2)*(dl/lc).^2) .* cos(2*k0 * dl);
plot(dl/lc, Iac, &#39;k&#39;)
title(&#39;(a) Interferogram&#39;)
xlabel(&#39;\Deltal/l_c&#39;)
ylabel(&#39;Signal&#39;)
axis([-0.6, 0.6, -1, 1])
subplot(4, 1 , 2 ) % rectified interferogram
Irec = abs(Iac);
plot(dl/lc, Irec, &#39;k&#39;)
title(&#39;(b) Rectified interferogram&#39;)
xlabel(&#39;\Deltal/l_c&#39;)
ylabel(&#39;Signal&#39;)
axis([-0.6, 0.6, -1, 1])
subplot(4, 1 , 3 ) % spectrum of the rectified interferogram
Frec1 = fft(Irec)/sqrt(N);
% order of frequencies: 0,1...(N/2-1),-N/2,-(N/2-1)...-1
Frec2 = fftshift(Frec1);
% shifted order of frequencies: -N/2,-(N/2-1)...-1, 0,1...(N/2-1)
dfreq = 1/(4*lc); % freq bin size = 1/sampling range
freq = dfreq*(-N/2:N/2-1); % frequency array
plot(freq*lambda0, abs(Frec2), &#39;k&#39;)
title(&#39;(c) Spectrum of the rectified interferogram&#39;)
xlabel(&#39;Frequency (1/\lambda_0)&#39;)
ylabel(&#39;Amplitude&#39;)
axis([-10, 10, 0, 5])
subplot(4, 1 , 4 ) % envelope
freq_cut = 1/lambda0/2; % cut-off frequency for filtering
i_cut = round(freq_cut/dfreq); % convert freq_cut to an array index
Ffilt = Frec1; % initialize array
Ffilt(i_cut:N-i_cut+1) = 0; % filter
Ifilt = abs(ifft(Ffilt))*sqrt(N); % amplitude of inverse FFT
plot(dl/lc, Ifilt/max(Ifilt), &#39;k&#39;)
Iac_en = exp(-16*log(2)*(dl/lc).^2); % envelope
hold on;
plot(dl(1:N/32:N)/lc, Iac_en(1:N/32:N), &#39;ko&#39;)
hold off;
title ( &#39; ( d ) Envelopes&#39;)
xlabel ( &#39; \Deltal/l_c &#39; )
ylabel(&#39;Signals&#39;)
axis ([ - 0.6 , 0.6, - 1 , 1])
legend(&#39;Demodulated&#39;,&#39;Original&#39;)FD-OCT(Python)
&#34;&#34;&#34;
Created on Wed Mar 18 20:25:19 2020

@author: Fuwei
&#34;&#34;&#34;
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.ticker import (AutoMinorLocator, MultipleLocator)
pi = np.pi
ln2 = np.log(2)
lambda_0 = 830e-9 # center wavelength of source
d_lambda = 20e-9 # FWHM wavelength bandwidth of source
ns=1.38 # refractive index of sample
l_s1 = 100e-6 # location of backscatterer 1
l_s2 = 150e-6 # location of backscatterer 2
r_s1 = 0.5 # reflectivity of backscatterer 1
r_s2 = 0.25 # reflectivity of backscatterer 2
k_0=2*pi/lambda_0 # center propagation constant
delta_k=2*pi*d_lambda/lambda_0**2 # FWHM bandwidth of k
sigma_k = delta_k/np.sqrt(2*ln2) # standard deviation of k
N=2**10 # number of sampling points
nsigma = 5 # number of standard deviations to plot on each side of kO

plt.figure(1)
plt.suptitle(&#34;FD-OCT Simulation&#34;,fontsize=25,fontweight=800)
plt.subplot(4,1,1) # Generate the interferogram
k = k_0 + sigma_k*np.linspace(-nsigma,nsigma, N) # array for k
S_k = np.exp(-(1/2)*((k-k_0)**2)/(sigma_k**2)) # Gaussian source PSD
E_s1 = r_s1*np.exp(2j*k*ns*l_s1) # sample electric field from scatter 1
E_s2 = r_s2*np.exp(2j*k*ns*l_s2) # sample electric field from scatter 2
I_k1 = S_k * np.abs(1 + E_s1 + E_s2)**2 # interferogram (r_R = 1)
plt.plot(k/k_0,I_k1/I_k1.max(), &#39;k&#39;)
plt.title(&#39;(a)Interferogram&#39;)
plt.xlabel(r&#39;Propagation constant $\mathregular{k/k_0}$&#39;)
plt.ylabel(&#39;Normalized intensity&#39;)
plt.axis([0.9,1.1,0,1.1])

plt.subplot(4,1,2) # Inverse Fourier transform (IFT) of the interferogram
spec1 =np.abs(np.fft.fftshift(np.fft.ifft(I_k1)))/np.sqrt(N)
d_ls_prime = 1/(2*nsigma*sigma_k/(2*pi)) # bin = 1/sampling range
ls_prime = d_ls_prime*np.arange(-N/2,N/2)# frequency array
plt.plot(ls_prime/(2*ns),spec1/spec1.max(), &#39;k&#39;) # scale the frequency
plt.title(&#39;(b)IFT of the interferogram&#39;)
plt.xlabel(r&#34;Depth $\mathregular{I_s}$ (m)&#34;)
plt.ylabel(&#39;Relative reflectivity&#39;)
plt.gca().xaxis.get_major_formatter().set_powerlimits((0,1))
plt.axis([-2*l_s2,2*l_s2,0,1.1])

plt.subplot(4,1,3); # IFT of the deconvolved interferogram
spec1_norm =np.abs(np.fft.fftshift(np.fft.ifft(I_k1/S_k)) )/np.sqrt(N)
d_ls_prime = 1/(2*nsigma*sigma_k/(2*pi)) # bin size = 1/sampling range
ls_prime = d_ls_prime*np.arange(-N/2,N/2) # frequency array
plt.plot(ls_prime/(2*ns),spec1_norm/spec1_norm.max(), &#39;k&#39;)
plt.title(&#39;IFT of the deconvolved interferogram&#39;)
plt.xlabel(r&#39;Depth $\mathregular{I_s}$ (m)&#39;)
plt.ylabel(&#39;Relative reflectivity&#39;)
plt.axis([-2*l_s2,2*l_s2,0,1.1])
plt.gca().xaxis.get_major_formatter().set_powerlimits((0,1))

plt.subplot(4,1,4) # IFT of the deconvolved differential interferogram
I_k2 = S_k * np.abs(-1 + E_s1 + E_s2)**2 # interferogram
delta_I_k = I_k1 - I_k2
spec2=np.abs(np.fft.fftshift(np.fft.ifft(delta_I_k/S_k)))/np.sqrt(N)
plt.plot(ls_prime/(2*ns),spec2/spec2.max(), &#39;k&#39;)
plt.title(&#39;IFT of the deconvolved differential interferogram&#39;)
plt.xlabel(&#34;Depth $\mathregular{I_s}$ (m)&#34;)
plt.ylabel(&#39;Relative reflectivity&#39;)
plt.axis([-2*l_s2,2*l_s2,0,1.1])
plt.gca().xaxis.get_major_formatter().set_powerlimits((0,1))
plt.subplots_adjust(hspace=0.6)
plt.show()FD-OCT(Matlab)
% Use SI units throughout
lambda0 = 830E-9; % center wavelength of source
dlambda = 20E-9; % FWHM wavelength bandwidth of source
ns=1.38; % refractive index of sample
ls1 = 100E-6; % location of backscatterer 1
ls2 = 150E-6; % location of backscatterer 2
rs1 = 0.5; % reflectivity of backscatterer 1
rs2 = 0.25; % reflectivity of backscatterer 2
k0=2*pi/lambda0; % center propagation constant
delta_k=2*pi*dlambda/lambda0^2; % FWHM bandwidth of k
sigma_k = delta_k/sqrt(2*log(2)); % standard deviation of k
N=2^10; % number of sampling points
nsigma = 5; % number of standard deviations to plot on each side of kO
subplot(4,1,1); % Generate the interferogram
k = k0 + sigma_k*linspace(-nsigma,nsigma, N); % array for k
S_k = exp(-(1/2)*(k-k0).^2/sigma_k^2); % Gaussian source PSD
E_s1 = rs1*exp(i*2*k*ns*ls1); % sample electric field from scatter 1
E_s2 = rs2*exp(i*2*k*ns*ls2); % sample electric field from scatter 2
I_k1 = S_k .* abs(1 + E_s1 + E_s2).^2; % interferogram (r_R = 1)
plot(k/k0,I_k1/max(I_k1), &#39;k&#39;);
title(&#39;Interferogram&#39;);
xlabel(&#39;Propagation constant k/k_0&#39;);
ylabel(&#39;Normalized intensity&#39;);
axis([0.9 1.1 0 1]);
subplot(4,1,2); % Inverse Fourier transform (IFT) of the interferogram
spec1 =abs(fftshift(ifft(I_k1)))/sqrt(N);
dls_prime = 1/(2*nsigma*sigma_k/(2*pi)); % bin = 1/sampling range
ls_prime = dls_prime*(-N/2:N/2-1); % frequency array
plot(ls_prime/(2*ns),spec1/max(spec1), &#39;k&#39;); % scale the frequency
title(&#39;IFT of the interferogram&#39;);
xlabel(&#39;Depth I_s (m)&#39;);
ylabel(&#39;Relative reflectivity&#39;);
axis([-2*ls2 2*ls2 0 1]);
subplot(4,1,3); % IFT of the deconvolved interferogram
spec1_norm =abs(fftshift(ifft(I_k1./S_k)) )/sqrt(N);
dls_prime = 1/(2*nsigma*sigma_k/(2*pi)); % bin size = 1/sampling range
ls_prime = dls_prime*(-N/2:N/2-1); % frequency array
plot(ls_prime/(2*ns),spec1_norm/max(spec1_norm), &#39;k&#39;);
title(&#39;IFT of the deconvolved interferogram&#39;);
xlabel(&#39;Depth I_s (m)&#39;);
ylabel(&#39;Relative reflectivity&#39;);
axis([-2*ls2 2*ls2 0 1]);
subplot(4,1,4); % IFT of the deconvolved differential interferogram
I_k2 = S_k .* abs(-1 + E_s1 + E_s2).^2; % interferogram
delta_I_k = I_k1 - I_k2;
spec2=abs(fftshift(ifft(delta_I_k./S_k)))/sqrt(N);
plot(ls_prime/(2*ns),spec2/max(spec2), &#39;k&#39;);
title(&#39;IFT of the deconvolved differential interferogram&#39;);
xlabel(&#34;Depth I_s ( m )&#34;);
ylabel(&#39;Relative reflectivity&#39;) ;
axis([-2*ls2 2*ls2 0 1]);<hr/>推荐书籍


原文地址：https://zhuanlan.zhihu.com/p/114015396