# Fourier-Domain Optical Coherence Tomography Signal Analysis and Numerical Modeling

Authors: J. Kalkman
Journal: International Journal of Optics
 (2017)

---

## Abstract

In this work the theory of the optical coherence tomography (OCT) signal after sampling, in dispersive media, with noise, and for a turbid medium is presented. The analytical theory is demonstrated with a one-dimensional numerical OCT model for (single) reflectors, dispersive media, and turbid media. For dispersive media the deterioration of the OCT axial resolution is quantified analytically and numerically. The OCT signal to noise ratio (SNR) is analyzed in the Fourier-domain and simulated with the numerical model. For an SNR based on the OCT intensity the conventional shot noise limited SNR is derived whereas for an SNR based on the OCT amplitude a 6.7 dB higher SNR is demonstrated. The OCT phase stability is derived in the Fourier-domain as 2SNR−1 and is validated using the numerical OCT model. The OCT single scattering model is simulated with the one-dimensional numerical model and applied to the simulation of an OCT image of a two-layered sample.


---

## Body

## 1. Introduction

Optical coherence tomography (OCT) is an optical imaging technique that is rapidly progressing into various application fields. Initially, OCT was invented for clinical diagnosis in the area of ophthalmology [1]. Currently it is used in medical application areas such as intravascular imaging and dermatology [2], as well as in various other nonmedical areas such as quality control [3], forensics [4], and biometry [5].The basic principle of OCT is to measure the time-of-flight of light echoes from tissue, which is done by creating an interference pattern between light propagating in the sample arm and light propagating in the reference arm of a Michelson interferometer. Initially, OCT depth scans were made by scanning the reference arm of the Michelson interferometer, so-called time-domain OCT. Later time-domain OCT was superseded by Fourier-domain OCT systems. Fourier-domain OCT systems are based on a measurement of the interference spectrum either in space on a spectrometer, this is called spectral-domain OCT (SD-OCT), or in time during the wavelength sweep of a rapidly tunable laser source, this is called swept-source OCT (SS-OCT). Subsequently, a depth scan is calculated by performing an inverse Fourier-transformation on the interference spectrum. Current state-of-the-art Fourier-domain OCT systems are capable of creating high quality images of in vivo tissue with micrometer resolution up to one to two millimeters deep. Individual depth scans can be acquired at multimegahertz rates.For both time-domain OCT and Fourier-domain OCT, in-depth knowledge of signal analysis and processing is of paramount importance for obtaining high quality images. There have been many reviews on OCT signal analysis and processing [6–9], but up till now no work has focused on the analytical theory in the z-k-domain and on showing the difference between intensity and amplitude based OCT signal analysis. Moreover, an easy-to-use numerical model of the OCT signal is lacking. Numerical simulations are an ideal tool to study the OCT signal as all Fourier-domain OCT processing steps are performed in the digital domain.In this study I give an overview of analytical theory of the OCT signal in Section2 that is followed by Section 3 in which the theoretical results are demonstrated with numerical simulation based on a discrete OCT signal model. In both chapters I discuss four topics. First, the OCT measurement process and sampling is described. Second, the effect of dispersion on the OCT signal is described. Third, the stochastic OCT signal and the determination of the signal to noise ratio and phase stability of the OCT signal is described. Fourth, the OCT signal in the single scattering approximation is described for a turbid medium.


## 2. OCT Theory

Figure1 shows a schematic of a Fourier-domain OCT system. The source has a spectral intensity distribution S(k) and emits electromagnetic waves into the interferometer. One-dimensional rectilinear propagation of plane scalar waves is assumed (i.e., polarization effects are neglected). The optical detection process averages over many optical cycles, hence the time dependence ejωt of the optical field is neglected. The incoming optical beam is split by an ideal beamsplitter with (intensity) reflection coefficient α and transmission coefficient 1-α. The reference arm field is reflected by the beamsplitter and travels a distance L in the reference arm before being reflected from the reference arm mirror. After propagation back over length L the field is transmitted through the beamsplitter and reaches the detector. The field transmitted by the beamsplitter travels in the sample arm to the zero-delay point (where both arms have equal length) and a further distance z where it is backscattered/reflected by the object. After traveling the same path back and being reflected by the beamsplitter, the field from the sample arm reaches the detector. The intensity as a function of wavenumber k is measured, as we assume here, with 100% efficiency on the detector. The detection of the intensity I(k) in space (SD-OCT) or in time (SS-OCT) is equivalent from a signal processing view.Figure 1
Schematic representation of the Fourier-domain OCT system. The length of the reference and sample arm is indicated withL and the unidirectional distance from the zero-delay point with z.


### 2.1. The Fourier-Domain OCT Signal

The source spectral intensity distributionS(k) is defined as a function of wavenumber k=2π/λ. The source spectral intensity distribution is normalized such that ∫S(k)dk=P0, where P0 is the source output power. The intensity launched into the interferometer is described by the product of a plane electromagnetic wave with its conjugate according to(1)Ik=SkejϕkSkejϕk∗,with ϕ(k) being a random phase and ∗ denoting complex conjugation. The phase ϕ(k) describes the initial randomly distributed phase of the radiators of the source. Assuming an ideal reference arm mirror with field reflectivity rr, the reference arm field Er(k) on the detector is given by(2)Erk=α1-αrrSkejϕkej2kL,where common paths traveled by both the sample and reference arm are neglected as their contribution is absent in the measured signal. The sample arm field on the detector is given by the summation of the field from all path lengths z~=2z of the sample (3)Esk=ejπα1-αSkejϕkej2kL·∫-∞∞az~2ejknzz~dz~,with a(z~/2)=a(z) being the complex valued depth dependent field reflection coefficient of the sample. The parameter z is defined as the unidirectional path length difference and is denoted as depth. The factor ejπ represents an additional phase factor for the sample arm traveling through the beamsplitter [10]. In this analysis it is initially assumed that n(z)=1 (i.e., the refractive index distribution of the sample is unity throughout the sample). The total intensity measured on the detector, (Er+Es)(Er+Es)∗, is composed of three contributions: the reference arm intensity, this is usually subtracted from the signal, the sample arm intensity, for turbid samples this is usually very small, and the interference intensity [11]. For the sake of clarity only the interference term is investigated, which is designated by I(k) as(4)Ik=ErkEsk∗+EskErk∗.Assuming that the sample is only located at z>0, then combining (4) with (2) and (3) results in(5)Ik=α1-αSkrre-jπ∫0∞a∗z~2e-jkz~dz~+ejπ∫0∞az~2ejkz~dz~.After some mathematical manipulation the intensity is found to be(6)Ik=α1-αSkrr∫-∞∞a~z~2e-jkz~dz~,with a~(z~)=e-jπa∗(z~/2)+ejπa(-z~/2) being a symmetric representation of the sample reflectivity. Equation (6) is recognized as the multiplication of the source spectrum S(k) with the Fourier transform of a~(z~/2) for the Fourier pair z~↔k. Hence, the complex valued depth dependent OCT signal i(z~) is obtained from an inverse Fourier transform of (6); that is,(7)iz~=Fk-1Ikz~=α1-αrrFk-1Skz~⊗a~z~2,where ⊗ indicates the convolution operation. Equation (7) states that the depth dependent OCT amplitude signal is given by the inverse Fourier transform of the source spectrum, the axial point spread function (PSF), convoluted with the sample reflectivity. After substitution of z~=2z the OCT signal is represented either by the amplitude |i(z)|=Re2⁡{i(z)}+Im2⁡{i(z)} or intensity i(z)2 of the complex valued signal.For a light source with arbitrary spectral shapeS(k), (7) can be used to calculate the expected OCT signal for an object with reflectivity a(z). The OCT axial point spread function (PSF) is given by the signal i(z) in response to a δ-function sample. For a Gaussian shaped spectrum S(k) with standard deviation σk and center wavenumber kc the OCT axial PSF is(8)iz=Fk-112πσk2exp⁡-k-kc22σk2z~z~=2z=e-jkc2z4πe-2z2σk2/2.From (8) it is determined that the full width at half maximum (FWHM) of the Gaussian point spread function in depth is FWHMz=σk-12ln⁡2. From FWHMk=2σk2ln⁡2 and the relation between the wavenumber k and wavelength λ the FWHMz is(9)FWHMz=2ln⁡2πλ2FWHMλ,which is generally known as the round trip coherence length [7]. Note that the expression in (9) not equal to the general definition of the coherence length [12]. In most papers the names axial resolution and coherence length are used interchangeably.


### 2.2. OCT Spectral Sampling and Detection

In Fourier-domain OCT the signal is either detected on spatially distributed pixels (spectral-domain OCT) or in time during a sweep of the laser wavelength (swept-source OCT). In both cases the optical signal is sampled and digitized. For both OCT modalities, spatial detection (SD-OCT) or temporal detection (SS-OCT) leads to an integration over wavenumber according to(10)Idetk=∫-∞∞Ik′Hk-k′dk′,with H(k) being the integration profile. Equation (10) is a convolution Idet(k)=I(k)⊗H(k).The detected continuous signalIdet(k) is subsequently sampled at regular intervals δk and digitized for storage in a computer. The sampled signal is represented as Isampled(k) and given by(11)Isampledk=∑m=-∞∞Ik⊗Hkδk-mδk,with m being an integer. The sampled OCT signal is given by inverse Fourier transform of (11). Sampling with spectral sampling rate δk leads to duplications of the z-domain signal shifted by multiples of δk [13]. Hence, the sampled OCT signal i(z) is(12)iz~=1δk∑m=-∞∞Fk-1Ikz~-m2πδk·Fk-1Hkz~-m2πδk.Sampling leads to a multitude of OCT depth scans i(z) that are separated by z=π/δk. To be able to separate these signals the maximum depth in the OCT signal i(z) should be halve this depth. Hence, the unidirectional maximum imaging depth is(13)zmax=π2δk.Using k=2π/λ and linearization, the well-known formula zmax=λ2/4δλ is found.Considering only the central part(m=0) of the inverse Fourier transform of (12), the OCT depth scan is(14)iz~=1δkFk-1Ikz~Fk-1Hkz~.The first inverse Fourier transform of the multiplication denotes the unsampled OCT signal conforming to (7). The second inverse Fourier transform describes the so-called roll-off of the OCT signal due to the detection process.In the case of spectral-domain OCT the signal is usually integrated over a square pixel of the camera with width and separationδk. The square pixel integration leads to a convolution in the k-domain with the filter(15)Hk=1δkrect⁡kδk,with “rect” being the rectangle function [14]. In case the pixel width is smaller than the pixel pitch δk, the function H(k) can be changed accordingly by substituting the spectral width of the pixel for δk in (15). In this case the roll-off of the OCT signal is reduced at the cost of less detected light. In addition, the spectrometer operates with a spectral resolution typically described by a Gaussian with standard deviation σr. The resolution integration leads to a convolution in the k-domain with the filter(16)Hk=12πσr2exp⁡-k22σr2.The filter operation in (15) and (16) in the k-domain is represented by a multiplication of the OCT signal in the z-domain by the functions(17)hz=12πsinc⁡δkzπ,(18)hz=12πexp⁡-2z2σr2,respectively. Note that the normalization with 2π-1 is due to the “physics” definition of the inverse Fourier transform that is used here.In SS-OCT the signal is integrated overk by a function determined by the detectors’ temporal response and the tuning speed. The spectral resolution is determined by the instantaneous line width of the tunable laser source [15].


#### 2.2.1. OCT Nonlinear Spectral Sampling

In OCT the interference signal is not measured in thek-domain directly, but is mapped to an intermediate coordinate, which is either space (SD-OCT) or time (SS-OCT). The interference signal is subsequently sampled at m discrete locations, with m an integer denoting the sample number ranging between 1 and M. Assuming that the m samples span a wavenumber range Mδk and that the wavenumber is linearly proportional with m the relation between wavenumber and sample index m can be written as(19)klinearm=k0+δkm.The intensities recorded Im at wavenumber klinearm are transformed to the z-domain using the discrete inverse Fourier transform. The spatial coordinate is determined from the linear k-domain and spans the range [-zmax,zmax] in M discrete steps.However, in general the relation betweenk and the intermediate coordinate is nonlinear. For simplicity it is assumed that the span over the detector remains fixed at Mδk. A simple nonlinear relation between the wavenumber k and the sample index m can be written as(20)knlm=k0+k1m+k2m2.The value of k2 is a parameter related to the amount of nonlinearity and can be freely chosen. The parameters k1 and k0 can be calculated1 such that the M samples span the bandwidth Mδk similar to the span for the spectrum sampled using (19). The spectrum I[m] sampled at linear k according to (19) is, most commonly, obtained from the measurements at the nonlinear k of (20) using interpolation.


### 2.3. OCT Signal in a Dispersive Medium

In the derivation of the OCT signal it was assumed that the refractive index of the sample is equal to unity. In general this is not the case and the refractive index is a function of both wavenumber and depth; that is,n=n(k,z). In case of a dispersive medium (6) changes to the form(21)Ik=α1-αrrSk∫-∞∞a~z~2e-jk∫0z~nk,z′dz′dz~,with n(k,z′) being the refractive index of the sample.Commonly the refractive index of a materialn(k) is expressed by the Sellmeier equation, which describes the refractive index with respect to the wavelength λ. To incorporate the refractive index in our OCT signal processing framework a polynomial expansion of n(k) around kc, the center wavenumber, is performed with(22)nk=∑p=0Pnpk-kcp.For mathematical clarity it is assumed that the refractive index is spatially invariant and completely described by the polynomial expansion of (22). Performing a Taylor expansion of the phase φ=kn(k)2z [16] (see Appendix A), the interference spectrum is(23)Ik=α1-αrrSk∫-∞∞e-jkcn0z~a~z~2∏p=0Pe-jnp+np+1z~k-kcp+1dz~.To gain insight into the effects of material dispersion on the OCT signal, the most simple sample is considered, a perfect reflector in a semi-infinite medium of refractive index n(k). The reflector is located at position d and described by a(z)=rsδ(z-d). From the reflectivity a~(z~/2) combined with (23) and after removal of common phase factors (these are absent after taking the amplitude of the inverse Fourier transform) the interference intensity is(24)Ik=α1-αrrrsSke-jkcn02d+ejkcn02d∏p=0Pe-jnp+np+12dk-kcp+1.Performing the inverse Fourier transform to z~ and using the shift in k property the OCT signal for z~>0 is(25)iz~=α1-αrrrsejkcz~Fk-1Skz~-n0+n12d⊗Fk-1∏p=1Pe-jnp+np+12dkp+1z~.The OCT signal peak is centered at position z=d(n0+n1); that is, the shift is proportional to the group index times the physical thickness d. Moreover, the width of the axial PSF increases due to a convolution with the inverse Fourier transform of a product of exponential functions. When the thickness is not zero (d≠0) and the material is dispersive (np≠0 for any p≥1) the convolution of the axial profile is not with a δ-function and consequently widens the OCT axial PSF. This effect is caused by the fact that the propagation speed of the wavelengths of the source are different in a dispersive material which results in their canceling effect for nonzero delays to be diminished.Assuming thatnp=0 for p≥2 (i.e., only linear dispersion (in k) is present), the last term of (25) is calculated analytically using the inverse Fourier transform of a complex valued Gaussian function(26)Fk-1e-jn12dk2z~z~=2z=12j2πn1dejz2/n12d.After performing the convolution and some algebra the width of the OCT axial PSF i(z) is(27)FWHMz=2ln⁡21σk2+16σk2n12d2.Comparing (27) with FWHMz=σk-12ln⁡2 derived previously, there is an additional term in the width due to the linear dispersion (n1) of the material, similar to the broadening factor described by others [17, 18]. Equation (27) demonstrates that even in a medium with only linear dispersion the OCT axial point spread function broadens.


### 2.4. OCT Signal to Noise Ratio

The OCT signals defined by the electric fields in (2) and (3) are deterministic signals. However, the total intensity emitted by the source is governed by random shot noise fluctuations. The presence of noise puts a limit on the OCT signal to noise ratio (SNR). The theoretical value for the OCT SNR is usually determined based on the following analysis.The OCT SNR is determined on a single depth scan basis with a mirror reflector in the sample arm. Using (6) for a single reflector a(z)=rsδ(z-z0), with rs being the mirror reflectivity, leads to a spectrum(28)Ik=α1-αSkrs2+rr2+2rsrrcos⁡k2z0.To theoretically determine the shot noise limited SNR an ideal rectangular spectrum S(k)=P0/Δkrect⁡k-kc/Δk is considered. Inserting this source spectrum into (28) and taking the inverse Fourier transform of both the signal and the noise, the OCT signal amplitude of the interference term is(29)iz=α1-αrsrrP02πsinc⁡Δkz±z0π.The factor (2π)-1 occurs for both the signal and the noise and can be disregarded. Hence, the peak OCT signal corresponds to a power P0α(1-α)rsrr for a total power received on the detector of P0α(1-α)(rs2+rr2).The fundamental noise limit in optical detection is determined by the shot noise of the photon arrival statistics. Considering only the effect of shot noise and assuming large number of photons, the shot noise statistics of the total number of photonsN is described by a Gaussian distributed random variable with mean N and standard deviation N. The peak signal in terms of detected number of photons is calculated as(30)Nsignal=P0τhνcα1-αrsrr,with νc=kcc/2π and τ being the detector integration time. In this definition the photon energy is calculated at the center frequency and not at every individual frequency. This is a good approximation since typically the optical bandwidth Δk is much smaller than the center wavenumber kc. The noise variance, which is equal to the number of photons, is estimated by summing the total optical power on the detector and converting it to the number of photons; that is,(31)Nnoise=σnoise2=P0τhνcα1-αrs2+rr2.Combining (30) and (31) and defining SNR as the square of the peak signal over the variance, the OCT shot noise limited SNR is(32)SNR=P0τhνcα1-αrs2rr2rs2+rr2.The numerator of (32) is proportional to the product of the number of photons from both the sample and the reference arm. In the denominator the variances of the various noise contributions (sample arm and reference arm) are added. In the limit of rr≫rs the SNR is equal to P0τ/hνcα(1-α)rs2; that is, it is only determined by the shot noise of the photon fluctuations in the sample arm [19].The OCT noise description presented here is based on the detected power, which is a valid description for the OCT intensity|i(z)|2. However, in many cases the OCT amplitude |i(z)| is used for analysis of the SNR and the phase stability analysis is performed on the complex OCT signal i(z). The amplitude based analysis is challenging to solve in the continuous signal description, therefore a discrete numerical OCT model is used to perform this analysis.


### 2.5. Single Scattering OCT Model

The single scattering OCT model is the most simple model of the OCT signal for a turbid medium sample [20, 21]. For a semi-infinite sample with the surface located at zs>0 the local unperturbed field at depth z, with z>zs, is determined by the attenuation of the optical field to position z. The transmission t(z) of the field to depth z, with z≥zs, is described by the modified Lambert-Beer law(33)tz=exp⁡-12μtz-zsuz-zs,with u(z) being the unit step function and the total attenuation coefficient μt=μs+μa equal to the sum of the scattering μs and absorption coefficient μa. In this description it is assumed that any scattering or absorption event removes the photon from the optical beam. The absorption and scattering coefficients are weakly dependent on k and this effect is ignored in the analysis. In the case of single scattering, a scattering event at depth z leads to some light being scattered back to the detector. The fraction of the total scattered intensity captured by the numerical aperture (NA) of the sample arm lens is described by the backscatter coefficient μb(34)μb=pNAμs=2πμs∫π-NAπpθsin⁡θdθ,where p(θ) is the (cylindrically symmetric) phase function. The field detected from this single scattering event is determined by the attenuation of the scattered field from depth z, through the sample, to the detector. Hence, the scattered field experiences a total transmission of t2(z)=exp⁡[-μt(z-zs)]. Consequently, the sample reflectivity for a semi-infinite turbid medium can be modeled as(35)az=μbexp⁡-μtz-zsuz-zs.The single scattering model can be extended by including confocal detection as a depth dependent light collection efficiency term [22].For typical OCT systems and samples it generally holds thatFWHMz≪μt-1. Consequently, the peak amplitude of the OCT signal is determined by setting the exponential term to unity. With this approximation, the convolution in (7) for an axial PSF well into the sample is an integration of the PSF multiplied by the backscatter coefficient. After performing this integration, the OCT signal i(z)2 for a semi-infinite turbid medium is(36)iz2=α21-α2rr2μbP028πσk2iz2=iz2=·exp⁡-2μtz-zsuz-zs.The height of the OCT signal is proportional to the square of the source power, proportional to the backscattering from the sample given by μb, and inversely proportional to the square of the coherence length.


## 2.1. The Fourier-Domain OCT Signal

The source spectral intensity distributionS(k) is defined as a function of wavenumber k=2π/λ. The source spectral intensity distribution is normalized such that ∫S(k)dk=P0, where P0 is the source output power. The intensity launched into the interferometer is described by the product of a plane electromagnetic wave with its conjugate according to(1)Ik=SkejϕkSkejϕk∗,with ϕ(k) being a random phase and ∗ denoting complex conjugation. The phase ϕ(k) describes the initial randomly distributed phase of the radiators of the source. Assuming an ideal reference arm mirror with field reflectivity rr, the reference arm field Er(k) on the detector is given by(2)Erk=α1-αrrSkejϕkej2kL,where common paths traveled by both the sample and reference arm are neglected as their contribution is absent in the measured signal. The sample arm field on the detector is given by the summation of the field from all path lengths z~=2z of the sample (3)Esk=ejπα1-αSkejϕkej2kL·∫-∞∞az~2ejknzz~dz~,with a(z~/2)=a(z) being the complex valued depth dependent field reflection coefficient of the sample. The parameter z is defined as the unidirectional path length difference and is denoted as depth. The factor ejπ represents an additional phase factor for the sample arm traveling through the beamsplitter [10]. In this analysis it is initially assumed that n(z)=1 (i.e., the refractive index distribution of the sample is unity throughout the sample). The total intensity measured on the detector, (Er+Es)(Er+Es)∗, is composed of three contributions: the reference arm intensity, this is usually subtracted from the signal, the sample arm intensity, for turbid samples this is usually very small, and the interference intensity [11]. For the sake of clarity only the interference term is investigated, which is designated by I(k) as(4)Ik=ErkEsk∗+EskErk∗.Assuming that the sample is only located at z>0, then combining (4) with (2) and (3) results in(5)Ik=α1-αSkrre-jπ∫0∞a∗z~2e-jkz~dz~+ejπ∫0∞az~2ejkz~dz~.After some mathematical manipulation the intensity is found to be(6)Ik=α1-αSkrr∫-∞∞a~z~2e-jkz~dz~,with a~(z~)=e-jπa∗(z~/2)+ejπa(-z~/2) being a symmetric representation of the sample reflectivity. Equation (6) is recognized as the multiplication of the source spectrum S(k) with the Fourier transform of a~(z~/2) for the Fourier pair z~↔k. Hence, the complex valued depth dependent OCT signal i(z~) is obtained from an inverse Fourier transform of (6); that is,(7)iz~=Fk-1Ikz~=α1-αrrFk-1Skz~⊗a~z~2,where ⊗ indicates the convolution operation. Equation (7) states that the depth dependent OCT amplitude signal is given by the inverse Fourier transform of the source spectrum, the axial point spread function (PSF), convoluted with the sample reflectivity. After substitution of z~=2z the OCT signal is represented either by the amplitude |i(z)|=Re2⁡{i(z)}+Im2⁡{i(z)} or intensity i(z)2 of the complex valued signal.For a light source with arbitrary spectral shapeS(k), (7) can be used to calculate the expected OCT signal for an object with reflectivity a(z). The OCT axial point spread function (PSF) is given by the signal i(z) in response to a δ-function sample. For a Gaussian shaped spectrum S(k) with standard deviation σk and center wavenumber kc the OCT axial PSF is(8)iz=Fk-112πσk2exp⁡-k-kc22σk2z~z~=2z=e-jkc2z4πe-2z2σk2/2.From (8) it is determined that the full width at half maximum (FWHM) of the Gaussian point spread function in depth is FWHMz=σk-12ln⁡2. From FWHMk=2σk2ln⁡2 and the relation between the wavenumber k and wavelength λ the FWHMz is(9)FWHMz=2ln⁡2πλ2FWHMλ,which is generally known as the round trip coherence length [7]. Note that the expression in (9) not equal to the general definition of the coherence length [12]. In most papers the names axial resolution and coherence length are used interchangeably.


## 2.2. OCT Spectral Sampling and Detection

In Fourier-domain OCT the signal is either detected on spatially distributed pixels (spectral-domain OCT) or in time during a sweep of the laser wavelength (swept-source OCT). In both cases the optical signal is sampled and digitized. For both OCT modalities, spatial detection (SD-OCT) or temporal detection (SS-OCT) leads to an integration over wavenumber according to(10)Idetk=∫-∞∞Ik′Hk-k′dk′,with H(k) being the integration profile. Equation (10) is a convolution Idet(k)=I(k)⊗H(k).The detected continuous signalIdet(k) is subsequently sampled at regular intervals δk and digitized for storage in a computer. The sampled signal is represented as Isampled(k) and given by(11)Isampledk=∑m=-∞∞Ik⊗Hkδk-mδk,with m being an integer. The sampled OCT signal is given by inverse Fourier transform of (11). Sampling with spectral sampling rate δk leads to duplications of the z-domain signal shifted by multiples of δk [13]. Hence, the sampled OCT signal i(z) is(12)iz~=1δk∑m=-∞∞Fk-1Ikz~-m2πδk·Fk-1Hkz~-m2πδk.Sampling leads to a multitude of OCT depth scans i(z) that are separated by z=π/δk. To be able to separate these signals the maximum depth in the OCT signal i(z) should be halve this depth. Hence, the unidirectional maximum imaging depth is(13)zmax=π2δk.Using k=2π/λ and linearization, the well-known formula zmax=λ2/4δλ is found.Considering only the central part(m=0) of the inverse Fourier transform of (12), the OCT depth scan is(14)iz~=1δkFk-1Ikz~Fk-1Hkz~.The first inverse Fourier transform of the multiplication denotes the unsampled OCT signal conforming to (7). The second inverse Fourier transform describes the so-called roll-off of the OCT signal due to the detection process.In the case of spectral-domain OCT the signal is usually integrated over a square pixel of the camera with width and separationδk. The square pixel integration leads to a convolution in the k-domain with the filter(15)Hk=1δkrect⁡kδk,with “rect” being the rectangle function [14]. In case the pixel width is smaller than the pixel pitch δk, the function H(k) can be changed accordingly by substituting the spectral width of the pixel for δk in (15). In this case the roll-off of the OCT signal is reduced at the cost of less detected light. In addition, the spectrometer operates with a spectral resolution typically described by a Gaussian with standard deviation σr. The resolution integration leads to a convolution in the k-domain with the filter(16)Hk=12πσr2exp⁡-k22σr2.The filter operation in (15) and (16) in the k-domain is represented by a multiplication of the OCT signal in the z-domain by the functions(17)hz=12πsinc⁡δkzπ,(18)hz=12πexp⁡-2z2σr2,respectively. Note that the normalization with 2π-1 is due to the “physics” definition of the inverse Fourier transform that is used here.In SS-OCT the signal is integrated overk by a function determined by the detectors’ temporal response and the tuning speed. The spectral resolution is determined by the instantaneous line width of the tunable laser source [15].


### 2.2.1. OCT Nonlinear Spectral Sampling

In OCT the interference signal is not measured in thek-domain directly, but is mapped to an intermediate coordinate, which is either space (SD-OCT) or time (SS-OCT). The interference signal is subsequently sampled at m discrete locations, with m an integer denoting the sample number ranging between 1 and M. Assuming that the m samples span a wavenumber range Mδk and that the wavenumber is linearly proportional with m the relation between wavenumber and sample index m can be written as(19)klinearm=k0+δkm.The intensities recorded Im at wavenumber klinearm are transformed to the z-domain using the discrete inverse Fourier transform. The spatial coordinate is determined from the linear k-domain and spans the range [-zmax,zmax] in M discrete steps.However, in general the relation betweenk and the intermediate coordinate is nonlinear. For simplicity it is assumed that the span over the detector remains fixed at Mδk. A simple nonlinear relation between the wavenumber k and the sample index m can be written as(20)knlm=k0+k1m+k2m2.The value of k2 is a parameter related to the amount of nonlinearity and can be freely chosen. The parameters k1 and k0 can be calculated1 such that the M samples span the bandwidth Mδk similar to the span for the spectrum sampled using (19). The spectrum I[m] sampled at linear k according to (19) is, most commonly, obtained from the measurements at the nonlinear k of (20) using interpolation.


## 2.2.1. OCT Nonlinear Spectral Sampling

In OCT the interference signal is not measured in thek-domain directly, but is mapped to an intermediate coordinate, which is either space (SD-OCT) or time (SS-OCT). The interference signal is subsequently sampled at m discrete locations, with m an integer denoting the sample number ranging between 1 and M. Assuming that the m samples span a wavenumber range Mδk and that the wavenumber is linearly proportional with m the relation between wavenumber and sample index m can be written as(19)klinearm=k0+δkm.The intensities recorded Im at wavenumber klinearm are transformed to the z-domain using the discrete inverse Fourier transform. The spatial coordinate is determined from the linear k-domain and spans the range [-zmax,zmax] in M discrete steps.However, in general the relation betweenk and the intermediate coordinate is nonlinear. For simplicity it is assumed that the span over the detector remains fixed at Mδk. A simple nonlinear relation between the wavenumber k and the sample index m can be written as(20)knlm=k0+k1m+k2m2.The value of k2 is a parameter related to the amount of nonlinearity and can be freely chosen. The parameters k1 and k0 can be calculated1 such that the M samples span the bandwidth Mδk similar to the span for the spectrum sampled using (19). The spectrum I[m] sampled at linear k according to (19) is, most commonly, obtained from the measurements at the nonlinear k of (20) using interpolation.


## 2.3. OCT Signal in a Dispersive Medium

In the derivation of the OCT signal it was assumed that the refractive index of the sample is equal to unity. In general this is not the case and the refractive index is a function of both wavenumber and depth; that is,n=n(k,z). In case of a dispersive medium (6) changes to the form(21)Ik=α1-αrrSk∫-∞∞a~z~2e-jk∫0z~nk,z′dz′dz~,with n(k,z′) being the refractive index of the sample.Commonly the refractive index of a materialn(k) is expressed by the Sellmeier equation, which describes the refractive index with respect to the wavelength λ. To incorporate the refractive index in our OCT signal processing framework a polynomial expansion of n(k) around kc, the center wavenumber, is performed with(22)nk=∑p=0Pnpk-kcp.For mathematical clarity it is assumed that the refractive index is spatially invariant and completely described by the polynomial expansion of (22). Performing a Taylor expansion of the phase φ=kn(k)2z [16] (see Appendix A), the interference spectrum is(23)Ik=α1-αrrSk∫-∞∞e-jkcn0z~a~z~2∏p=0Pe-jnp+np+1z~k-kcp+1dz~.To gain insight into the effects of material dispersion on the OCT signal, the most simple sample is considered, a perfect reflector in a semi-infinite medium of refractive index n(k). The reflector is located at position d and described by a(z)=rsδ(z-d). From the reflectivity a~(z~/2) combined with (23) and after removal of common phase factors (these are absent after taking the amplitude of the inverse Fourier transform) the interference intensity is(24)Ik=α1-αrrrsSke-jkcn02d+ejkcn02d∏p=0Pe-jnp+np+12dk-kcp+1.Performing the inverse Fourier transform to z~ and using the shift in k property the OCT signal for z~>0 is(25)iz~=α1-αrrrsejkcz~Fk-1Skz~-n0+n12d⊗Fk-1∏p=1Pe-jnp+np+12dkp+1z~.The OCT signal peak is centered at position z=d(n0+n1); that is, the shift is proportional to the group index times the physical thickness d. Moreover, the width of the axial PSF increases due to a convolution with the inverse Fourier transform of a product of exponential functions. When the thickness is not zero (d≠0) and the material is dispersive (np≠0 for any p≥1) the convolution of the axial profile is not with a δ-function and consequently widens the OCT axial PSF. This effect is caused by the fact that the propagation speed of the wavelengths of the source are different in a dispersive material which results in their canceling effect for nonzero delays to be diminished.Assuming thatnp=0 for p≥2 (i.e., only linear dispersion (in k) is present), the last term of (25) is calculated analytically using the inverse Fourier transform of a complex valued Gaussian function(26)Fk-1e-jn12dk2z~z~=2z=12j2πn1dejz2/n12d.After performing the convolution and some algebra the width of the OCT axial PSF i(z) is(27)FWHMz=2ln⁡21σk2+16σk2n12d2.Comparing (27) with FWHMz=σk-12ln⁡2 derived previously, there is an additional term in the width due to the linear dispersion (n1) of the material, similar to the broadening factor described by others [17, 18]. Equation (27) demonstrates that even in a medium with only linear dispersion the OCT axial point spread function broadens.


## 2.4. OCT Signal to Noise Ratio

The OCT signals defined by the electric fields in (2) and (3) are deterministic signals. However, the total intensity emitted by the source is governed by random shot noise fluctuations. The presence of noise puts a limit on the OCT signal to noise ratio (SNR). The theoretical value for the OCT SNR is usually determined based on the following analysis.The OCT SNR is determined on a single depth scan basis with a mirror reflector in the sample arm. Using (6) for a single reflector a(z)=rsδ(z-z0), with rs being the mirror reflectivity, leads to a spectrum(28)Ik=α1-αSkrs2+rr2+2rsrrcos⁡k2z0.To theoretically determine the shot noise limited SNR an ideal rectangular spectrum S(k)=P0/Δkrect⁡k-kc/Δk is considered. Inserting this source spectrum into (28) and taking the inverse Fourier transform of both the signal and the noise, the OCT signal amplitude of the interference term is(29)iz=α1-αrsrrP02πsinc⁡Δkz±z0π.The factor (2π)-1 occurs for both the signal and the noise and can be disregarded. Hence, the peak OCT signal corresponds to a power P0α(1-α)rsrr for a total power received on the detector of P0α(1-α)(rs2+rr2).The fundamental noise limit in optical detection is determined by the shot noise of the photon arrival statistics. Considering only the effect of shot noise and assuming large number of photons, the shot noise statistics of the total number of photonsN is described by a Gaussian distributed random variable with mean N and standard deviation N. The peak signal in terms of detected number of photons is calculated as(30)Nsignal=P0τhνcα1-αrsrr,with νc=kcc/2π and τ being the detector integration time. In this definition the photon energy is calculated at the center frequency and not at every individual frequency. This is a good approximation since typically the optical bandwidth Δk is much smaller than the center wavenumber kc. The noise variance, which is equal to the number of photons, is estimated by summing the total optical power on the detector and converting it to the number of photons; that is,(31)Nnoise=σnoise2=P0τhνcα1-αrs2+rr2.Combining (30) and (31) and defining SNR as the square of the peak signal over the variance, the OCT shot noise limited SNR is(32)SNR=P0τhνcα1-αrs2rr2rs2+rr2.The numerator of (32) is proportional to the product of the number of photons from both the sample and the reference arm. In the denominator the variances of the various noise contributions (sample arm and reference arm) are added. In the limit of rr≫rs the SNR is equal to P0τ/hνcα(1-α)rs2; that is, it is only determined by the shot noise of the photon fluctuations in the sample arm [19].The OCT noise description presented here is based on the detected power, which is a valid description for the OCT intensity|i(z)|2. However, in many cases the OCT amplitude |i(z)| is used for analysis of the SNR and the phase stability analysis is performed on the complex OCT signal i(z). The amplitude based analysis is challenging to solve in the continuous signal description, therefore a discrete numerical OCT model is used to perform this analysis.


## 2.5. Single Scattering OCT Model

The single scattering OCT model is the most simple model of the OCT signal for a turbid medium sample [20, 21]. For a semi-infinite sample with the surface located at zs>0 the local unperturbed field at depth z, with z>zs, is determined by the attenuation of the optical field to position z. The transmission t(z) of the field to depth z, with z≥zs, is described by the modified Lambert-Beer law(33)tz=exp⁡-12μtz-zsuz-zs,with u(z) being the unit step function and the total attenuation coefficient μt=μs+μa equal to the sum of the scattering μs and absorption coefficient μa. In this description it is assumed that any scattering or absorption event removes the photon from the optical beam. The absorption and scattering coefficients are weakly dependent on k and this effect is ignored in the analysis. In the case of single scattering, a scattering event at depth z leads to some light being scattered back to the detector. The fraction of the total scattered intensity captured by the numerical aperture (NA) of the sample arm lens is described by the backscatter coefficient μb(34)μb=pNAμs=2πμs∫π-NAπpθsin⁡θdθ,where p(θ) is the (cylindrically symmetric) phase function. The field detected from this single scattering event is determined by the attenuation of the scattered field from depth z, through the sample, to the detector. Hence, the scattered field experiences a total transmission of t2(z)=exp⁡[-μt(z-zs)]. Consequently, the sample reflectivity for a semi-infinite turbid medium can be modeled as(35)az=μbexp⁡-μtz-zsuz-zs.The single scattering model can be extended by including confocal detection as a depth dependent light collection efficiency term [22].For typical OCT systems and samples it generally holds thatFWHMz≪μt-1. Consequently, the peak amplitude of the OCT signal is determined by setting the exponential term to unity. With this approximation, the convolution in (7) for an axial PSF well into the sample is an integration of the PSF multiplied by the backscatter coefficient. After performing this integration, the OCT signal i(z)2 for a semi-infinite turbid medium is(36)iz2=α21-α2rr2μbP028πσk2iz2=iz2=·exp⁡-2μtz-zsuz-zs.The height of the OCT signal is proportional to the square of the source power, proportional to the backscattering from the sample given by μb, and inversely proportional to the square of the coherence length.


## 3. OCT Simulations

The OCT signal analysis in the continuousz-k-domain presented in Section 2 is implemented in discrete form using software written in MATLAB (MathWorks, R2016) (some examples of the simulations can be found at [23]). The simulations are performed using the parameters in Table 1 unless indicated otherwise. The input source spectrum is modeled as a Gaussian shaped discrete spectrum S[m] with M samples, total power P0, and normalization ∑m=1M‍S[m]=P0, with m an integer. The interferometer is represented by the intensity splitting ratio α and reference arm field reflectivity rr. The sample is represented by a mirror with field reflectivity rs. The center wavelength λc is converted to center wavenumber kc, the FWHM optical bandwidth on the detector Δλ is converted to optical bandwidth in wavenumber Δk. The spectra are sampled at M points between kc±1.5Δk. A quasi-continuous k-axis is calculated by upsampling the sampled k by a factor of 8. The sampled z-axis is calculated by distributing M points over the range [-zmax,zmax]. The interferometric signal is constructed according to (4) with the reference arm field given by (2) and the sample arm field given by (3).Table 1
Parameters used in the OCT simulations.
Parameter
Value
Description
λc
800 nm
Center wavelength
Δλ
50 nm
FWHM optical bandwidth
P0
1 mW
Input power
τ
1 ms
Integration time
α
0.5
Intensity splitting ratio
M
1024
Number of samples


### 3.1. OCT Spectral Sampling and Detection Modeling

A simulation of the effects of sampling and pixel integration is demonstrated in Figure2. The mirror reflector is placed at distances between z = 0.2 mm and z = 1.5 mm, with the maximum imaging depth equal to zmax = 1.2 mm. The quasi-continuous k-domain signal is convoluted using a discretized version of the filter in (15) and (18). Subsequently, the filtered k-domain signal is resampled at samples separated by δk=3Δk/M. The k-domain signals are subsequently transformed to the z-domain using the discrete inverse Fourier transform.Figure 2
The effect of sampling and detection on the OCT signal. (a) For pixel integration the fringe contrast in thek-domain is reduced relative to the input spectrum (blue dashed line) after convolution in the k-domain in the continuous (blue continuous line) and sampled signal (red circles). (b) Due to pixel integration the signal in the z-domain decreases in depth for the continuous (blue continuous line) and sampled (red circles) signal. The green curve indicates the OCT signal roll-off of (17). (c) Similar to (a) but for the effect of spectral integration with a resolution of σr=4δk. (d) Similar to (b) but for the effect of spectral integration.
(a)
(b)
(c)
(d)Figure2(a) shows the Gaussian shaped interference spectrum in (quasi-) continuous k and at the M sampled k points for a reflector at 0.85 mm distance and for the case of pixel integration and sampling. The fringe contrast in I(k) is reduced compared to the envelope of the spectrum of the ideal δ-function sampled signal. Figure 2(b) shows the OCT signal in the z-domain for multiple reflector positions. The OCT amplitude decreases for increasing optical path length due to the k-domain filtering of the signal. Also shown is the theoretical roll-off of the OCT signal with depth according to the sinc function of (17). The sampled and quasi-continuous signal follow the theoretical roll-off. However, when the mirror reflector is placed at a distance larger than zmax, calculated according to (13), aliasing takes place. This can be observed for the mirror at a distance of 1.5 mm in the quasi-continuous z-domain signal with the sampled OCT signal peak aliased to the distance of 0.9 mm.Simulations of the OCT signal for the effect of spectral resolution and sampling are shown in Figures2(c) and 2(d). In this case the k-domain signal I(k) is convoluted with a Gaussian filter according to (16) as shown in Figure 2(c), which leads to a Gaussian roll-off according to (18) as shown in Figure 2(d).


#### 3.1.1. OCT Spectral Resampling Modeling

In Figure3 the effect of nonlinear k sampling of the interference spectrum on the OCT signal is demonstrated. Figure 3(a) shows that for a spectrum that is sampled nonlinearly in k the OCT signal for a mirror, that is, the axial PSF, broadens in depth. This effect is quantified in Figure 3(b) showing the OCT axial FWHM in depth as compared to the bandwidth limited resolution that is obtained for a linear k sampled signal. Also shown are the OCT axial FWHM for resampled spectral data using different interpolation methods. Most of the broadening is corrected, however, at large depth the FWHM increases, especially for the more simple interpolation methods. This effect is attributed to small interpolation errors that occur at high frequency fringe modulations close to the maximum imaging depth.Figure 3
(a) OCT signal acquired with lineark (blue) and nonlinear k (red). (b) OCT axial PSF FWHM for different spectral interpolation algorithms. The lines are a guide to the eye.
(a)
(b)


### 3.2. OCT Dispersion Modeling

The effect of dispersion on the OCT signal is modeled by assuming that light in the sample arm propagates through a water layer with thicknessd where it is reflected by an ideal mirror. The refractive index of water is described by the Sellmeier equation(37)n2-1=∑p=1∞Apλ2λ2-λp2.The Sellmeier coefficients for water at a temperature of 20°C [24] are summarized in Table 2. The refractive index is fitted with a linear model n(k)=n0+n1(k-kc) resulting in n0 = 1.3309 and n1=9.6226·10-10 m rad-1. Figure 4(a) shows the refractive index of water and the linear fit.Table 2
Sellmeier parameters for water [24].
Parameter
Value
A1
5.684027565
⋅
1
0
-
1
A2
1.726177391
⋅
1
0
-
1
A3
2.086189578
⋅
1
0
-
2
λ12
5.101829712
⋅
1
0
-
3
μm2
λ22
1.821153936
⋅
1
0
-
2
μm2
λ32
2.620722293
⋅
1
0
-
2
μm2Figure 4
(a) Refractive index of water versus wavenumber (blue line) and linear fit (dashed red line). (b) OCT peak position versus thickness of water (red circles) and theory (blue line). Also indicated is the displacement in air (green dashed line). (c) FWHM peak width of the OCT signal versus water thickness (red circles) and theory (blue line). Also indicated is the dispersion-free FWHM peak width (green line).
(a)
(b)
(c)The effect of dispersion on the OCT signal is investigated using the linear model of the refractive index of water. After propagation of the light through the water and OCT signal construction, the peak of the OCT signal is fitted with a Gaussian function with center position and standard deviation as fit parameters. The OCT signal is upsampled by a factor of 8M in the k-domain signal to obtain a more accurate estimation of the center position and width of the OCT axial PSF. The fitted standard deviation is transformed to the FWHMz of the OCT axial PSF in air (n=1). Figure 4(b) shows the effect of the material refractive index on the center position of the OCT PSF and the comparison with the theory. The peak location is located at a depth z=(n0+n1)d; that is, the depth is equal to the group index multiplied with the physical water thickness, in accordance with (25). Figure 4(c) shows the effect of material dispersion on the FWHM of the OCT axial PSF. The OCT peak width broadens with increasing water thickness. The simulation is compared to the theoretical prediction of (27) and perfect agreement is observed.


### 3.3. OCT Signal to Noise Ratio Modeling

For the signal to noise analysis the OCT signal is modeled in thek-domain by a square spectrum with M channels that is normalized according to ∑m=1M‍Sm=P0; that is, the power in every channel is Sm=P0/M. Furthermore, without loss of generality, it is assumed that the sample is a mirror positioned such that the peak of z-domain signal is position q0; that is, the ratio M/q0 is an integer. The OCT interferometric signal is described in terms of detected number of photons at every channel for the interferometric part of the signal, which, according to (28), is(38)Nsignalm=P0τhνcMα1-αrsrr2cos⁡2πmq0M.Performing an inverse DFT of (38) results in a fully real signal with two peaks with amplitude(39)Nsignal=P0τhνcMα1-αrsrr;see Appendix B for a derivation. Next, the noise is considered in the absence of the interferometric signal. The detected power of the sample and reference arm in every channel, according to (28), is converted to mean number of noise photons(40)Nnoisem=P0τhνcMα1-αrs2+rr2.The shot noise in every detector channel is modeled as an independent white noise Gaussian distributed random variable N[m] with mean and variance Nnoise, with Nnoise the number of photons on the mth detector element described by (40). Multiple (5000) independent noise realizations of the k-domain OCT signal are generated. For every realization N[m], the field is calculated from the intensity and propagated through the interferometer similar to as described in Section 3.1. The simulation is performed for a fixed sample arm field reflectivity rs and varying reference arm field reflectivity rr. The OCT signal is obtained from the interference signal by a discrete inverse Fourier transform according to(41)iq=1M∑m=0M-1Nmej2πqm/M,with q being the integer valued depth variable. The OCT signal iq is a complex valued discrete random variable and the OCT signal is calculated by taking either the amplitude iq or the intensity iq2.ForN[m] a zero mean random variable, that is, with the mean number of background photons subtracted, it is shown in Appendix C that the real part, Re⁡{i[q]}, and the imaginary part, Im⁡{i[q]}, of the signal are uncorrelated Gaussian random variables. For a measurement in the absence of signal, the real and imaginary parts of iq have a variance equal to Nnoise/2M; see Appendix D.The construction of the OCT intensityi(z)2 from the real and imaginary part leads to a variable iq2 with a negative exponential distribution [25] with a standard deviation(42)σiq2=2σRe⁡iq2=NnoiseM.For an intensity based OCT signal iq2 the SNR is defined as the peak signal, that is, the square of (39), over the standard deviation of the noise intensity. Combining (38), (40), and (42), the SNR is equal to the usual expression of (32).The construction of the OCT amplitudei(z) from the real and imaginary part leads to a variable i[q] with a Rayleigh distribution [25] with a variance(43)σiq2=2-π2σRe⁡iq2=Nnoise2M2-π2.Consequently, for an OCT SNR defined from the amplitude in terms of the square of the peak amplitude iq2 over the noise variance of i[q], the fundamental OCT SNR limit is obtained by combining (38), (40), and (43) and results in(44)SNRA=P0τhνcα1-αrs2rr2rs2+rr211-π/4,where the subscript A indicates amplitude. Note that the small reduction of the difference between peak value and noise floor, due to the increase of the mean noise floor, is neglected. Hence, the SNRA is a factor ~4.7 higher than the usually stated SNR [19]. After subtraction of the noise signal with its mean and calculation of the amplitude, the noise is transformed to a Rayleigh distributed random variable which results in a reduction of the variance with a factor of 1-π/4. Based on the expression 10 10log SNRA, the SNR increases by approximately 6.7 dB compared to the SNR for the intensity based signal.Figure5 shows the performance of an OCT system analyzed for signal amplitude with only shot noise present. The sample arm reflectivity is fixed at rs=10-2 and the reference arm reflectivity is varied between rr=10-5 and 1. Figure 5(a) shows the square of the OCT peak amplitude and its comparison to the theoretical predication. Conforming to (39) the peak OCT signal increases linearly with the reference arm reflectivity rr2. Figure 5(b) shows the noise variance of the OCT amplitude for varying reference arm power for the simulated noise signal and for the analytical expression of (40). At small reference arm power the shot noise is entirely dominated by the contribution from the sample arm and hence is constant. At large reference arm powers the sample arm noise is negligible and the noise increases linearly with reference arm power rr2. In between these two regimes there is a transition region. Figure 5(c) shows the simulated SNR for the OCT amplitude signal and the theoretical prediction. The small discrepancy between the theory and the simulation is due to a small overestimation of the noise in the simulated data.Figure 5
(a) Simulation of the peak of the OCT signal (red circles) compared with the model (blue line). (b) Simulation of the noise of the OCT signal (red circles) compared with the model (blue line). (c) OCT signal to noise ratio simulation (red circles) and the model (blue line).
(a)
(b)
(c)The phase stability is determined from the OCT signal in the complexz-domain plane, as indicated in Figure 6. The OCT signal is represented with a vector with length equal to the peak amplitude of the OCT signal and angle 2πmq0/M. Because of the background subtraction the noise is located at the origin of the complex plane. The phase θ of the complex valued OCT signal i[q] is determined in the complex plane with(45)θ=arctan⁡Im⁡iqRe⁡iq.The real and imaginary part of the OCT signal with noise are given by Re⁡iq±σRe⁡iq and Im⁡iq±σIm⁡iq, with the noise variance given by Nnoise/2M. Generalized error analysis, as derived in Appendix E, is applied to calculate the phase uncertainty of (45). The variance of Var(θ) of the angle is(46)Var⁡θ=σRe⁡iq2Re2⁡iq+Im2⁡iq.Using the expression for the variance of the real part of the signal in (D.2) and the square of the signal in (39), the phase variance is derived as(47)Var⁡θ=rr2+rs22P0τ/hνcα1-αrr2rs2=12SNR.Hence, the phase variance of the OCT signal is equal to 2SNR−1 using the standard SNR definition of (32).Figure 6
Schematic of the construction of the OCT phase in thez-domain. The OCT signal is represented by the long arrow. A single noise realization is represented by the short arrow and is a vector from the Gaussian distributed noise around the origin. The OCT signal including the noise is constructed from the vectorial addition of the two contributions and leads to the variation Δθ of the angle θ.Figure7 shows the simulated phase variance as a function of SNR for the same data as in Figure 5. The simulation demonstrates the 2SNR−1 dependence of (47). The inset shows the Gaussian distribution of the OCT phase θ for rr=10-2.Figure 7
Simulation of the phase variance of the OCT signal (circles) as a function of SNR compared to the theoretical predication of (47) (solid line). The inset shows the histogram of the phase distribution of the OCT signal at rr=10-2.


### 3.4. Single Scattering OCT Signal

OCT deals mainly with imaging in turbid media and the process of light collection from a sample is a complicated 3D scattering problem. However, for sample volumes that are small and far away from the sample arm lens a 1D description can be made. In the simulations, the scattering volume is assumed to be a cube with sideslx=ly=l=25μm, and depth lz between 0.2 and 0.8 mm, see Figure 8. The three-dimensional scattering problem is converted to a one-dimensional model by considering the particles, located in the 3D volume l2lz, to be distributed randomly over a depth range lz at positions zi and described by δ-function reflectors. A numerical 1D OCT model is calculated for a sample consisting of stationary spherical particles with a radius of 1 μm, refractive index npart=1.4, and volume fraction of fv=0.01. The refractive index of the medium is nmed=1.33. Based on the volume fraction and particle radius the concentration of particles C is determined with the total number of particles in the simulation given by Npart=Cl2lz=1194. The scattering properties of the particles are calculated using Mie-scattering [26]. From the scattering efficiency Q the scattering coefficient μs=Qπr2C is calculated. The total power scattered Ps from the first particle is Ps=IinQπr2=PinQπr2/l2, with Iin and Pin the intensity and power incident on the first particle, respectively. Hence, the power transmission after a single scattering event is T=(P0-Ps)/P0. Consequently, the unscattered power remaining after propagation distance lz is the power after transmission through Npart particles; that is,(48)Plz=PinTNpart=Pin1-Qπr2l2Clzl2.Using the relation(49)limn→∞⁡1+xnn=ex, it can be approximated that for many particles P(lz)=Pinexp⁡(-μslz), in agreement with the single scattering model. Absorption can be treated in a similar fashion by multiplying the transmitted power with an additional factor exp⁡(-μalz).Figure 8
Schematic overview of the conversion of the three-dimensional distribution of particles (a) to a one-dimensional distribution in depth (b).
(a)
(b)The reflection of every particle is determined from the collected reflected intensityPr from every particle, which is Pr=PinpNAQπr2/l2, with Pin the power incident on the particle and pNA the fraction of the scattered intensity captured by the collection lens [27]. Hence, the intensity reflection coefficient per scattering event is(50)R=pNAQπr2l2. The field detected from a single particle in the sample arm is given by the field reaching the particle, multiplied by the particle’s reflection coefficient and the field transmitted through the sample back to the detector. Hence, the detected field from the pth particle is(51)Ei=PinTp-1RTp-1.In the model, every particle has a detected field given by the amplitude of (51) and a delay determined by its position zi. A simulation of the OCT intensity i(z)2 in the single scattering approximation is shown in Figure 9 where it is averaged over 50 independent realizations of the random particle positions. The OCT signal in Figure 9(a) is in agreement with the single scattering model as the OCT signal decays according to the single scattering model A02e-2μsz. The height of the numerical OCT signal is(52)A02=2P02M2α21-α22NpartδzlzpNAQπr2, with the factor M-2 from the discrete inverse Fourier transform. The factor 2 originates from the speckle statistics of the random phasor sum [25]. The number of particles (phasors) in a single depth bin is Npartδz/lz, with δz the width of a single depth bin. Figure 9(b) shows the OCT interference spectrum for the simulated sample and Figure 9(c) shows the OCT intensity distribution, evaluated at the point indicated in Figure 9(a), for 2000 realizations of the OCT signal. As expected the distribution of the OCT intensity is an exponential function with a contrast ratio of 1.0528, close to the theoretical limit of 1 for fully developed speckle.Figure 9
Simulation of the intensity OCT signal. (a) Simulated average OCT signal in depth. The dashed red line indicates the single scattering OCT model, and the arrow indicates the location where the distribution of the OCT signal is determined. (b) A spectrum for a single A-line simulation. (c) Distribution of OCT intensities.
(a)
(b)
(c)To demonstrate the use of the numerical OCT model, a two-layered piece of tissue is simulated; see Figure10. Layer 1 is 300 μm thick and consists of 500 nm radius particles at 4% volume fraction, with nmed=1.33 and npart=1.35. Layer 2 is 200 μm thick and consists of 1 μm particles at 2% volume fraction, with nmed=1.33 and npart=1.5.Figure 10
Simulation of the OCT signal intensity for a two-layer turbid sample. The intensity scale is logarithmic.


## 3.1. OCT Spectral Sampling and Detection Modeling

A simulation of the effects of sampling and pixel integration is demonstrated in Figure2. The mirror reflector is placed at distances between z = 0.2 mm and z = 1.5 mm, with the maximum imaging depth equal to zmax = 1.2 mm. The quasi-continuous k-domain signal is convoluted using a discretized version of the filter in (15) and (18). Subsequently, the filtered k-domain signal is resampled at samples separated by δk=3Δk/M. The k-domain signals are subsequently transformed to the z-domain using the discrete inverse Fourier transform.Figure 2
The effect of sampling and detection on the OCT signal. (a) For pixel integration the fringe contrast in thek-domain is reduced relative to the input spectrum (blue dashed line) after convolution in the k-domain in the continuous (blue continuous line) and sampled signal (red circles). (b) Due to pixel integration the signal in the z-domain decreases in depth for the continuous (blue continuous line) and sampled (red circles) signal. The green curve indicates the OCT signal roll-off of (17). (c) Similar to (a) but for the effect of spectral integration with a resolution of σr=4δk. (d) Similar to (b) but for the effect of spectral integration.
(a)
(b)
(c)
(d)Figure2(a) shows the Gaussian shaped interference spectrum in (quasi-) continuous k and at the M sampled k points for a reflector at 0.85 mm distance and for the case of pixel integration and sampling. The fringe contrast in I(k) is reduced compared to the envelope of the spectrum of the ideal δ-function sampled signal. Figure 2(b) shows the OCT signal in the z-domain for multiple reflector positions. The OCT amplitude decreases for increasing optical path length due to the k-domain filtering of the signal. Also shown is the theoretical roll-off of the OCT signal with depth according to the sinc function of (17). The sampled and quasi-continuous signal follow the theoretical roll-off. However, when the mirror reflector is placed at a distance larger than zmax, calculated according to (13), aliasing takes place. This can be observed for the mirror at a distance of 1.5 mm in the quasi-continuous z-domain signal with the sampled OCT signal peak aliased to the distance of 0.9 mm.Simulations of the OCT signal for the effect of spectral resolution and sampling are shown in Figures2(c) and 2(d). In this case the k-domain signal I(k) is convoluted with a Gaussian filter according to (16) as shown in Figure 2(c), which leads to a Gaussian roll-off according to (18) as shown in Figure 2(d).


### 3.1.1. OCT Spectral Resampling Modeling

In Figure3 the effect of nonlinear k sampling of the interference spectrum on the OCT signal is demonstrated. Figure 3(a) shows that for a spectrum that is sampled nonlinearly in k the OCT signal for a mirror, that is, the axial PSF, broadens in depth. This effect is quantified in Figure 3(b) showing the OCT axial FWHM in depth as compared to the bandwidth limited resolution that is obtained for a linear k sampled signal. Also shown are the OCT axial FWHM for resampled spectral data using different interpolation methods. Most of the broadening is corrected, however, at large depth the FWHM increases, especially for the more simple interpolation methods. This effect is attributed to small interpolation errors that occur at high frequency fringe modulations close to the maximum imaging depth.Figure 3
(a) OCT signal acquired with lineark (blue) and nonlinear k (red). (b) OCT axial PSF FWHM for different spectral interpolation algorithms. The lines are a guide to the eye.
(a)
(b)


## 3.1.1. OCT Spectral Resampling Modeling

In Figure3 the effect of nonlinear k sampling of the interference spectrum on the OCT signal is demonstrated. Figure 3(a) shows that for a spectrum that is sampled nonlinearly in k the OCT signal for a mirror, that is, the axial PSF, broadens in depth. This effect is quantified in Figure 3(b) showing the OCT axial FWHM in depth as compared to the bandwidth limited resolution that is obtained for a linear k sampled signal. Also shown are the OCT axial FWHM for resampled spectral data using different interpolation methods. Most of the broadening is corrected, however, at large depth the FWHM increases, especially for the more simple interpolation methods. This effect is attributed to small interpolation errors that occur at high frequency fringe modulations close to the maximum imaging depth.Figure 3
(a) OCT signal acquired with lineark (blue) and nonlinear k (red). (b) OCT axial PSF FWHM for different spectral interpolation algorithms. The lines are a guide to the eye.
(a)
(b)


## 3.2. OCT Dispersion Modeling

The effect of dispersion on the OCT signal is modeled by assuming that light in the sample arm propagates through a water layer with thicknessd where it is reflected by an ideal mirror. The refractive index of water is described by the Sellmeier equation(37)n2-1=∑p=1∞Apλ2λ2-λp2.The Sellmeier coefficients for water at a temperature of 20°C [24] are summarized in Table 2. The refractive index is fitted with a linear model n(k)=n0+n1(k-kc) resulting in n0 = 1.3309 and n1=9.6226·10-10 m rad-1. Figure 4(a) shows the refractive index of water and the linear fit.Table 2
Sellmeier parameters for water [24].
Parameter
Value
A1
5.684027565
⋅
1
0
-
1
A2
1.726177391
⋅
1
0
-
1
A3
2.086189578
⋅
1
0
-
2
λ12
5.101829712
⋅
1
0
-
3
μm2
λ22
1.821153936
⋅
1
0
-
2
μm2
λ32
2.620722293
⋅
1
0
-
2
μm2Figure 4
(a) Refractive index of water versus wavenumber (blue line) and linear fit (dashed red line). (b) OCT peak position versus thickness of water (red circles) and theory (blue line). Also indicated is the displacement in air (green dashed line). (c) FWHM peak width of the OCT signal versus water thickness (red circles) and theory (blue line). Also indicated is the dispersion-free FWHM peak width (green line).
(a)
(b)
(c)The effect of dispersion on the OCT signal is investigated using the linear model of the refractive index of water. After propagation of the light through the water and OCT signal construction, the peak of the OCT signal is fitted with a Gaussian function with center position and standard deviation as fit parameters. The OCT signal is upsampled by a factor of 8M in the k-domain signal to obtain a more accurate estimation of the center position and width of the OCT axial PSF. The fitted standard deviation is transformed to the FWHMz of the OCT axial PSF in air (n=1). Figure 4(b) shows the effect of the material refractive index on the center position of the OCT PSF and the comparison with the theory. The peak location is located at a depth z=(n0+n1)d; that is, the depth is equal to the group index multiplied with the physical water thickness, in accordance with (25). Figure 4(c) shows the effect of material dispersion on the FWHM of the OCT axial PSF. The OCT peak width broadens with increasing water thickness. The simulation is compared to the theoretical prediction of (27) and perfect agreement is observed.


## 3.3. OCT Signal to Noise Ratio Modeling

For the signal to noise analysis the OCT signal is modeled in thek-domain by a square spectrum with M channels that is normalized according to ∑m=1M‍Sm=P0; that is, the power in every channel is Sm=P0/M. Furthermore, without loss of generality, it is assumed that the sample is a mirror positioned such that the peak of z-domain signal is position q0; that is, the ratio M/q0 is an integer. The OCT interferometric signal is described in terms of detected number of photons at every channel for the interferometric part of the signal, which, according to (28), is(38)Nsignalm=P0τhνcMα1-αrsrr2cos⁡2πmq0M.Performing an inverse DFT of (38) results in a fully real signal with two peaks with amplitude(39)Nsignal=P0τhνcMα1-αrsrr;see Appendix B for a derivation. Next, the noise is considered in the absence of the interferometric signal. The detected power of the sample and reference arm in every channel, according to (28), is converted to mean number of noise photons(40)Nnoisem=P0τhνcMα1-αrs2+rr2.The shot noise in every detector channel is modeled as an independent white noise Gaussian distributed random variable N[m] with mean and variance Nnoise, with Nnoise the number of photons on the mth detector element described by (40). Multiple (5000) independent noise realizations of the k-domain OCT signal are generated. For every realization N[m], the field is calculated from the intensity and propagated through the interferometer similar to as described in Section 3.1. The simulation is performed for a fixed sample arm field reflectivity rs and varying reference arm field reflectivity rr. The OCT signal is obtained from the interference signal by a discrete inverse Fourier transform according to(41)iq=1M∑m=0M-1Nmej2πqm/M,with q being the integer valued depth variable. The OCT signal iq is a complex valued discrete random variable and the OCT signal is calculated by taking either the amplitude iq or the intensity iq2.ForN[m] a zero mean random variable, that is, with the mean number of background photons subtracted, it is shown in Appendix C that the real part, Re⁡{i[q]}, and the imaginary part, Im⁡{i[q]}, of the signal are uncorrelated Gaussian random variables. For a measurement in the absence of signal, the real and imaginary parts of iq have a variance equal to Nnoise/2M; see Appendix D.The construction of the OCT intensityi(z)2 from the real and imaginary part leads to a variable iq2 with a negative exponential distribution [25] with a standard deviation(42)σiq2=2σRe⁡iq2=NnoiseM.For an intensity based OCT signal iq2 the SNR is defined as the peak signal, that is, the square of (39), over the standard deviation of the noise intensity. Combining (38), (40), and (42), the SNR is equal to the usual expression of (32).The construction of the OCT amplitudei(z) from the real and imaginary part leads to a variable i[q] with a Rayleigh distribution [25] with a variance(43)σiq2=2-π2σRe⁡iq2=Nnoise2M2-π2.Consequently, for an OCT SNR defined from the amplitude in terms of the square of the peak amplitude iq2 over the noise variance of i[q], the fundamental OCT SNR limit is obtained by combining (38), (40), and (43) and results in(44)SNRA=P0τhνcα1-αrs2rr2rs2+rr211-π/4,where the subscript A indicates amplitude. Note that the small reduction of the difference between peak value and noise floor, due to the increase of the mean noise floor, is neglected. Hence, the SNRA is a factor ~4.7 higher than the usually stated SNR [19]. After subtraction of the noise signal with its mean and calculation of the amplitude, the noise is transformed to a Rayleigh distributed random variable which results in a reduction of the variance with a factor of 1-π/4. Based on the expression 10 10log SNRA, the SNR increases by approximately 6.7 dB compared to the SNR for the intensity based signal.Figure5 shows the performance of an OCT system analyzed for signal amplitude with only shot noise present. The sample arm reflectivity is fixed at rs=10-2 and the reference arm reflectivity is varied between rr=10-5 and 1. Figure 5(a) shows the square of the OCT peak amplitude and its comparison to the theoretical predication. Conforming to (39) the peak OCT signal increases linearly with the reference arm reflectivity rr2. Figure 5(b) shows the noise variance of the OCT amplitude for varying reference arm power for the simulated noise signal and for the analytical expression of (40). At small reference arm power the shot noise is entirely dominated by the contribution from the sample arm and hence is constant. At large reference arm powers the sample arm noise is negligible and the noise increases linearly with reference arm power rr2. In between these two regimes there is a transition region. Figure 5(c) shows the simulated SNR for the OCT amplitude signal and the theoretical prediction. The small discrepancy between the theory and the simulation is due to a small overestimation of the noise in the simulated data.Figure 5
(a) Simulation of the peak of the OCT signal (red circles) compared with the model (blue line). (b) Simulation of the noise of the OCT signal (red circles) compared with the model (blue line). (c) OCT signal to noise ratio simulation (red circles) and the model (blue line).
(a)
(b)
(c)The phase stability is determined from the OCT signal in the complexz-domain plane, as indicated in Figure 6. The OCT signal is represented with a vector with length equal to the peak amplitude of the OCT signal and angle 2πmq0/M. Because of the background subtraction the noise is located at the origin of the complex plane. The phase θ of the complex valued OCT signal i[q] is determined in the complex plane with(45)θ=arctan⁡Im⁡iqRe⁡iq.The real and imaginary part of the OCT signal with noise are given by Re⁡iq±σRe⁡iq and Im⁡iq±σIm⁡iq, with the noise variance given by Nnoise/2M. Generalized error analysis, as derived in Appendix E, is applied to calculate the phase uncertainty of (45). The variance of Var(θ) of the angle is(46)Var⁡θ=σRe⁡iq2Re2⁡iq+Im2⁡iq.Using the expression for the variance of the real part of the signal in (D.2) and the square of the signal in (39), the phase variance is derived as(47)Var⁡θ=rr2+rs22P0τ/hνcα1-αrr2rs2=12SNR.Hence, the phase variance of the OCT signal is equal to 2SNR−1 using the standard SNR definition of (32).Figure 6
Schematic of the construction of the OCT phase in thez-domain. The OCT signal is represented by the long arrow. A single noise realization is represented by the short arrow and is a vector from the Gaussian distributed noise around the origin. The OCT signal including the noise is constructed from the vectorial addition of the two contributions and leads to the variation Δθ of the angle θ.Figure7 shows the simulated phase variance as a function of SNR for the same data as in Figure 5. The simulation demonstrates the 2SNR−1 dependence of (47). The inset shows the Gaussian distribution of the OCT phase θ for rr=10-2.Figure 7
Simulation of the phase variance of the OCT signal (circles) as a function of SNR compared to the theoretical predication of (47) (solid line). The inset shows the histogram of the phase distribution of the OCT signal at rr=10-2.


## 3.4. Single Scattering OCT Signal

OCT deals mainly with imaging in turbid media and the process of light collection from a sample is a complicated 3D scattering problem. However, for sample volumes that are small and far away from the sample arm lens a 1D description can be made. In the simulations, the scattering volume is assumed to be a cube with sideslx=ly=l=25μm, and depth lz between 0.2 and 0.8 mm, see Figure 8. The three-dimensional scattering problem is converted to a one-dimensional model by considering the particles, located in the 3D volume l2lz, to be distributed randomly over a depth range lz at positions zi and described by δ-function reflectors. A numerical 1D OCT model is calculated for a sample consisting of stationary spherical particles with a radius of 1 μm, refractive index npart=1.4, and volume fraction of fv=0.01. The refractive index of the medium is nmed=1.33. Based on the volume fraction and particle radius the concentration of particles C is determined with the total number of particles in the simulation given by Npart=Cl2lz=1194. The scattering properties of the particles are calculated using Mie-scattering [26]. From the scattering efficiency Q the scattering coefficient μs=Qπr2C is calculated. The total power scattered Ps from the first particle is Ps=IinQπr2=PinQπr2/l2, with Iin and Pin the intensity and power incident on the first particle, respectively. Hence, the power transmission after a single scattering event is T=(P0-Ps)/P0. Consequently, the unscattered power remaining after propagation distance lz is the power after transmission through Npart particles; that is,(48)Plz=PinTNpart=Pin1-Qπr2l2Clzl2.Using the relation(49)limn→∞⁡1+xnn=ex, it can be approximated that for many particles P(lz)=Pinexp⁡(-μslz), in agreement with the single scattering model. Absorption can be treated in a similar fashion by multiplying the transmitted power with an additional factor exp⁡(-μalz).Figure 8
Schematic overview of the conversion of the three-dimensional distribution of particles (a) to a one-dimensional distribution in depth (b).
(a)
(b)The reflection of every particle is determined from the collected reflected intensityPr from every particle, which is Pr=PinpNAQπr2/l2, with Pin the power incident on the particle and pNA the fraction of the scattered intensity captured by the collection lens [27]. Hence, the intensity reflection coefficient per scattering event is(50)R=pNAQπr2l2. The field detected from a single particle in the sample arm is given by the field reaching the particle, multiplied by the particle’s reflection coefficient and the field transmitted through the sample back to the detector. Hence, the detected field from the pth particle is(51)Ei=PinTp-1RTp-1.In the model, every particle has a detected field given by the amplitude of (51) and a delay determined by its position zi. A simulation of the OCT intensity i(z)2 in the single scattering approximation is shown in Figure 9 where it is averaged over 50 independent realizations of the random particle positions. The OCT signal in Figure 9(a) is in agreement with the single scattering model as the OCT signal decays according to the single scattering model A02e-2μsz. The height of the numerical OCT signal is(52)A02=2P02M2α21-α22NpartδzlzpNAQπr2, with the factor M-2 from the discrete inverse Fourier transform. The factor 2 originates from the speckle statistics of the random phasor sum [25]. The number of particles (phasors) in a single depth bin is Npartδz/lz, with δz the width of a single depth bin. Figure 9(b) shows the OCT interference spectrum for the simulated sample and Figure 9(c) shows the OCT intensity distribution, evaluated at the point indicated in Figure 9(a), for 2000 realizations of the OCT signal. As expected the distribution of the OCT intensity is an exponential function with a contrast ratio of 1.0528, close to the theoretical limit of 1 for fully developed speckle.Figure 9
Simulation of the intensity OCT signal. (a) Simulated average OCT signal in depth. The dashed red line indicates the single scattering OCT model, and the arrow indicates the location where the distribution of the OCT signal is determined. (b) A spectrum for a single A-line simulation. (c) Distribution of OCT intensities.
(a)
(b)
(c)To demonstrate the use of the numerical OCT model, a two-layered piece of tissue is simulated; see Figure10. Layer 1 is 300 μm thick and consists of 500 nm radius particles at 4% volume fraction, with nmed=1.33 and npart=1.35. Layer 2 is 200 μm thick and consists of 1 μm particles at 2% volume fraction, with nmed=1.33 and npart=1.5.Figure 10
Simulation of the OCT signal intensity for a two-layer turbid sample. The intensity scale is logarithmic.


## 4. Discussion and Conclusion

In this article I presented a theoretical and numerical analysis of the most important signal processing steps in Fourier-domain OCT. This OCT analysis is based on a comparison of the signals in both thek- and z-domains.Linear spectral sampling and detection is theoretically described and numerically simulated. Good agreement is observed between the analytical model and the numerical simulations. The OCT signal roll-off described here has been demonstrated by numerous groups for SD-OCT (e.g., in [28]), or for SS-OCT (e.g., in [29]). In almost all cases the OCT interference spectrum is nonlinearly sampled. The resulting deterioration of the axial resolution can be removed using nonlinear fast Fourier transforms or, as is most common, linearization of the k-domain signal using numerical interpolation [30]. For nonlinear k-domain sampling, the Nyquist depth limit is dependent on the wavenumber as the spectral sampling rate varies over the spectrum. This partial aliasing effect [31] results in an OCT signal drop at large depths.In general, the effect of nonlinear sampling on thek-integration as presented for the case of pixel integration and spectral resolution is not addressed. In most OCT signal analyses the roll-off is described by a k-invariant convolution over the wavenumbers, which corresponds to a multiplication of the OCT signal in the z-domain. For small amounts of spectral nonlinearity this is, in general, seen as sufficient to characterize the OCT signal. The quantification of the full effect can be implemented with the presented numerical OCT model by applying a k-variant convolution.The effect of dispersion on the OCT axial PSF shows that even for a material with only linear dispersion ink, broadening of the OCT axial PSF takes place. Higher order dispersions are easily implemented in the numerical simulations by providing the full material dispersion as described by the Sellmeier equation.The simulations of the OCT SNR are in good correspondence with the analytical theory. The same behavior of OCT SNR versus reference arm reflectivity is demonstrated as, for example, measured by Grulkowski et al. [32] and Leitgeb et al. [19], although in experimental settings also other noise sources play a role. In contrast to the usual OCT SNR analysis, which is based on detected power, a k-domain representation of the signal to noise is presented for an (ideal) square source spectrum. Using the numerical model it is demonstrated how the shot noise of the light detected in the k-domain is transformed through the inverse Fourier-transformation to the intensity and amplitude of the complex z-domain OCT signal. In this analysis it is shown that for an SNR based on the OCT amplitude, the fundamental shot noise limit is a factor 1-π/4-1 higher than for an intensity based OCT signal analysis [19]. In case of an intensity based OCT signal, the experimentally determined SNR can be obtained close to the theoretical limit [33]. For a more realistic Gaussian shaped spectrum the OCT SNR is generally expected to be lower due to the less efficient distribution of optical noise over the detector elements.From thez-domain signal and noise description a rigorous derivation of the OCT phase stability is made. It is derived that the absolute OCT phase has a variance equal to 2SNR−1, where the SNR is defined based on the intensity. In this derivation no use was made of the approximation that the signal is much larger than the noise [34] or that the noise is orthogonal to the signal [35]. The derived result is similar to that obtained by Park et al. [34] and Vakoc et al. [36]. However, it differs from the result of Choma et al. [35], which seems to have an additional factor 1/2 in the signal component.The numerical model is developed and applied to simulate the OCT signal of a semi-infinite turbid medium. For a semi-infinite turbid medium the simulation matches the single scattering model. The origin of the exponential decay of the OCT signal is well reproduced by modeling the light in the sample after multiple transmission events. The OCT intensity has an exponential distribution [37], whereas the amplitude has a Rayleigh distribution, similar to what has been shown by [38]. The OCT signal intensity from the numerical model incorporates a factor 2 originating from the speckle distribution, which needs to be included in the analytical single scattering OCT model of (36).The one-dimensional OCT model accurately describes OCT measurements of low scattering media for the attenuation [20, 27] and the speckle statistics [39]. The model is simple to use and can easily be adapted for testing OCT attenuation quantification [40] or tissue segmentation algorithms [41]. Although it assumes light from the sample arm to be incident perpendicular to the sample, the effect of focusing and back-coupling efficiency [22] can be easily implemented by adding a depth dependent confocal back-coupling function. The numerical OCT model for a turbid medium does incorporate the effect of dependent scattering in its dependence on μs [42], however, multiple scattering effects are not incorporated. More elaborate analytical models [43] or Monte Carlo simulations [44] can be used to study the OCT signal in these cases [45]. The OCT model can be extended to incorporate time dependent scattering processes such as present in the case of Doppler OCT and speckle dynamics [46].


## 5. Conclusion

In conclusion, I presented an overview of analytical expressions for the Fourier-domain OCT signal after sampling, in dispersive media, with noise, and for a scattering medium. A numerical model is developed to simulate the OCT signal. Good agreement is observed between analytical and numerical results.



---

*Source: 9586067-2017-03-22.xml*