Fourier-Domain Optical Coherence Tomography Signal Analysis and Numerical Modeling

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.


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-offlight 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.Fourierdomain 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 Fouriertransformation on the interference spectrum.Current stateof-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][7][8][9], but up till now no work has focused on the analytical theory in the --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 Section 2 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.

OCT Theory
Figure 1 shows a schematic of a Fourier-domain OCT system.The source has a spectral intensity distribution () and emits electromagnetic waves into the interferometer.Onedimensional 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 e  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  in the reference arm before being reflected from the reference arm mirror.After propagation back over length  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  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  is measured, as we assume here, with 100% efficiency on the detector.The detection of the intensity () in space (SD-OCT) or in time (SS-OCT) is equivalent from a signal processing view.

The Fourier-Domain OCT Signal.
The source spectral intensity distribution () is defined as a function of wavenumber  = 2/.The source spectral intensity distribution is normalized such that ∫ ()d =  0 , where  0 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 with () being a random phase and * denoting complex conjugation.The phase () describes the initial randomly distributed phase of the radiators of the source.Assuming an ideal reference arm mirror with field reflectivity   , the reference arm field   () on the detector is given by 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 = 2 of the sample with (z/2) = () being the complex valued depth dependent field reflection coefficient of the sample.The parameter  is defined as the unidirectional path length difference and is denoted as depth.The factor   represents an additional phase factor for the sample arm traveling through the beamsplitter [10].In this analysis it is initially assumed that () = 1 (i.e., the refractive index distribution of the sample is unity throughout the sample).The total intensity measured on the detector, (  +   )(  +   ) * , 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 () as Assuming that the sample is only located at  > 0, then combining (4) with ( 2) and (3) results in After some mathematical manipulation the intensity is found to be with ã(z) =  −  * (z/2) +   (−z/2) being a symmetric representation of the sample reflectivity.Equation ( 6) is recognized as the multiplication of the source spectrum () with the Fourier transform of ã(z/2) for the Fourier pair z ↔ .Hence, the complex valued depth dependent OCT signal (z) is obtained from an inverse Fourier transform of (6); that is, 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 = 2 the OCT signal is represented either by the amplitude |()| = √Re 2 {()} + Im 2 {()} or intensity |()| 2 of the complex valued signal.For a light source with arbitrary spectral shape (), (7) can be used to calculate the expected OCT signal for an object with reflectivity ().The OCT axial point spread function (PSF) is given by the signal () in response to a -function sample.For a Gaussian shaped spectrum () with standard deviation   and center wavenumber   the OCT axial PSF is .
From ( 8) it is determined that the full width at half maximum (FWHM) of the Gaussian point spread function in depth is From FWHM  = 2  √ 2 ln 2 and the relation between the wavenumber  and wavelength  the FWHM  is 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.

OCT Spectral Sampling and Detection.
In Fourierdomain 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 with () being the integration profile.Equation ( 10) is a convolution  det () = () ⊗ ().
The detected continuous signal  det () is subsequently sampled at regular intervals  and digitized for storage in a computer.The sampled signal is represented as  sampled () and given by with  being an integer.The sampled OCT signal is given by inverse Fourier transform of (11).Sampling with spectral sampling rate  leads to duplications of the -domain signal shifted by multiples of  [13].Hence, the sampled OCT signal () is Sampling leads to a multitude of OCT depth scans () that are separated by  = /.To be able to separate these signals the maximum depth in the OCT signal () should be halve this depth.Hence, the unidirectional maximum imaging depth is Using  = 2/ and linearization, the well-known formula  max =  2 /4 is found.
Considering only the central part ( = 0) of the inverse Fourier transform of ( 12), the OCT depth scan is 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 rolloff 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 .The square pixel integration leads to a convolution in the -domain with the filter with "rect" being the rectangle function [14].In case the pixel width is smaller than the pixel pitch , the function () can be changed accordingly by substituting the spectral width of the pixel for  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   .The resolution integration leads to a convolution in the domain with the filter The filter operation in (15) and ( 16) in the -domain is represented by a multiplication of the OCT signal in the domain by the functions 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 over  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].

OCT Nonlinear Spectral Sampling.
In OCT the interference signal is not measured in the -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  discrete locations, with  an integer denoting the sample number ranging between 1 and .Assuming that the  samples span a wavenumber range  and that the wavenumber is linearly proportional with  the relation between wavenumber and sample index  can be written as The intensities recorded [] at wavenumber  linear [] are transformed to the -domain using the discrete inverse Fourier transform.The spatial coordinate is determined from the linear -domain and spans the range [− max ,  max ] in  discrete steps.However, in general the relation between  and the intermediate coordinate is nonlinear.For simplicity it is assumed that the span over the detector remains fixed at .A simple nonlinear relation between the wavenumber  and the sample index  can be written as The value of  2 is a parameter related to the amount of nonlinearity and can be freely chosen.The parameters  1 and  0 can be calculated 1 such that the  samples span the bandwidth  similar to the span for the spectrum sampled using (19).The spectrum [] sampled at linear  according to ( 19) is, most commonly, obtained from the measurements at the nonlinear  of (20) using interpolation.

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,  = (, ).In case of a dispersive medium (6) changes to the form with (,   ) being the refractive index of the sample.
Commonly the refractive index of a material () 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 () around   , the center wavenumber, is performed with 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  = ()2 [16] (see Appendix A), the interference spectrum is 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 ().
The reflector is located at position  and described by () =   ( − ).From the reflectivity ã(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 Performing the inverse Fourier transform to z and using the shift in  property the OCT signal for z > 0 is The OCT signal peak is centered at position  = ( 0 +  1 ); that is, the shift is proportional to the group index times the physical thickness .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 ( ̸ = 0) and the material is dispersive (  ̸ = 0 for any  ≥ 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 that   = 0 for  ≥ 2 (i.e., only linear dispersion (in ) is present), the last term of ( 25) is calculated analytically using the inverse Fourier transform of a complex valued Gaussian function After performing the convolution and some algebra the width of the OCT axial PSF |()| is Comparing (27) with FWHM  =  −1  √ 2 ln 2 derived previously, there is an additional term in the width due to the linear dispersion ( 1 ) 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.

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 () =   ( −  0 ), with   being the mirror reflectivity, leads to a spectrum To theoretically determine the shot noise limited SNR an ideal rectangular spectrum () = ( 0 /Δ)rect(( −   )/Δ) 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 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  0 (1 − )    for a total power received on the detector of  0 (1 − )( 2  +  2  ).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 photons  is described by a Gaussian distributed random variable with mean  and standard deviation √ .The peak signal in terms of detected number of photons is calculated as with ]  =   /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 Δ is much smaller than the center wavenumber   .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, Combining ( 30) and ( 31) and defining SNR as the square of the peak signal over the variance, the OCT shot noise limited SNR is 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   ≫   the SNR is equal to ( 0 /ℎ]  )(1 − ) 2  ; 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 |()| 2 .However, in many cases the OCT amplitude |()| is used for analysis of the SNR and the phase stability analysis is performed on the complex OCT signal ().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.

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   > 0 the local unperturbed field at depth , with  >   , is determined by the attenuation of the optical field to position .The transmission () of the field to depth , with  ≥   , is described by the modified Lambert-Beer law with u() being the unit step function and the total attenuation coefficient   =   +   equal to the sum of the scattering   and absorption coefficient   .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  and this effect is ignored in the analysis.In the case of single scattering, a scattering event at depth  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 where () 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 , through the sample, to the detector.Hence, the scattered field experiences a total transmission of  2 () = exp[−  ( −   )].Consequently, the sample reflectivity for a semi-infinite turbid medium can be modeled as 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 that FWHM  ≪  −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 |()| 2 for a semi-infinite turbid medium is The height of the OCT signal is proportional to the square of the source power, proportional to the backscattering from the sample given by   , and inversely proportional to the square of the coherence length.

OCT Simulations
The OCT signal analysis in the continuous --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 [] with  samples, total power  0 , and normalization ∑  =1 [] =  0 , with  an integer.The interferometer is represented by the intensity splitting ratio  and reference arm field reflectivity   .The sample is represented by a mirror with field reflectivity   .The center wavelength   is converted to center wavenumber   , the FWHM optical bandwidth on the detector Δ is converted to optical bandwidth in wavenumber Δ.The spectra are sampled at  points between   ± 1.5Δ.A quasi-continuous -axis is calculated by upsampling the sampled  by a factor of 8.The sampled -axis is calculated by distributing  points over the range [− max ,  max ].The interferometric signal is constructed according to (4) with the reference arm field given by ( 2) and the sample arm field given by (3).

OCT Spectral Sampling and Detection Modeling.
A simulation of the effects of sampling and pixel integration is demonstrated in Figure 2. The mirror reflector is placed at distances between  = 0.2 mm and  = 1.5 mm, with the maximum imaging depth equal to  max = 1.2 mm.The quasicontinuous -domain signal is convoluted using a discretized version of the filter in (15) and (18).Subsequently, the filtered -domain signal is resampled at samples separated by  = 3Δ/.The -domain signals are subsequently transformed to the -domain using the discrete inverse Fourier transform.
Figure 2(a) shows the Gaussian shaped interference spectrum in (quasi-) continuous  and at the  sampled  points for a reflector at 0.85 mm distance and for the case of pixel integration and sampling.The fringe contrast in () is reduced compared to the envelope of the spectrum of the ideal -function sampled signal.Figure 2(b) shows the OCT signal in the -domain for multiple reflector positions.The OCT amplitude decreases for increasing optical path length due to the -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 quasicontinuous signal follow the theoretical roll-off.However, when the mirror reflector is placed at a distance larger than  max , 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 -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 Figures 2(c) and 2(d).In this case the -domain signal () 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).

OCT Spectral Resampling Modeling.
In Figure 3 the effect of nonlinear  sampling of the interference spectrum on the OCT signal is demonstrated.Figure 3(a) shows that for a spectrum that is sampled nonlinearly in  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  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.

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 thickness  where it is reflected by an ideal mirror.The refractive index of water is described by the Sellmeier equation  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 () =  0 +  1 ( −   ) resulting in  0 = 1.3309 and  1 = 9.6226 ⋅ 10 −10 m rad −1 .Figure 4(a) shows the refractive index of water and the linear fit.
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 International Journal of Optics  Table 2: Sellmeier parameters for water [24].
deviation as fit parameters.The OCT signal is upsampled by a factor of 8 in the -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 FWHM  of the OCT axial PSF in air ( = 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  = ( 0 +  1 ); 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.

OCT Signal to Noise Ratio Modeling.
For the signal to noise analysis the OCT signal is modeled in the -domain by a square spectrum with  channels that is normalized according to ∑  =1 [] =  0 ; that is, the power in every channel is [] =  0 /.Furthermore, without loss of generality, it is assumed that the sample is a mirror positioned such that the peak of -domain signal is position  0 ; that is, the ratio / 0 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 Performing an inverse DFT of ( 38) results in a fully real signal with two peaks with amplitude 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 The shot noise in every detector channel is modeled as an independent white noise Gaussian distributed random variable N[] with mean and variance  noise , with  noise the number of photons on the th detector element described by (40).Multiple (5000) independent noise realizations of the -domain OCT signal are generated.For every realization [], 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   and varying reference arm field reflectivity   .The OCT signal is obtained from the interference signal by a discrete inverse Fourier transform according to with  being the integer valued depth variable.The OCT signal i[] is a complex valued discrete random variable and the OCT signal is calculated by taking either the amplitude |i For N[] 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[]}, and the imaginary part, Im{i[]}, of the signal are uncorrelated Gaussian random variables.For a measurement in the absence of signal, the real and imaginary parts of i[] have a variance equal to  noise /2; see Appendix D.
The construction of the OCT intensity |()| 2 from the real and imaginary part leads to a variable |i[]| 2 with a negative exponential distribution [25] with a standard deviation For an intensity based OCT signal |i[]| 2 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 amplitude |()| from the real and imaginary part leads to a variable |i[]| with a Rayleigh distribution [25] with a variance Consequently, for an OCT SNR defined from the amplitude in terms of the square of the peak amplitude |[]| 2 over the noise variance of |i[]|, the fundamental OCT SNR limit is obtained by combining (38), (40), and (43) and results in where the subscript  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 SNR  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 10 log SNR  , the SNR increases by approximately 6.7 dB compared to the SNR for the intensity based signal.
Figure 5 shows the performance of an OCT system analyzed for signal amplitude with only shot noise present.The sample arm reflectivity is fixed at   = 10 −2 and the reference arm reflectivity is varied between   = 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  2  .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  2  .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.
The phase stability is determined from the OCT signal in the complex -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 0 /.Because of the background subtraction the noise is located at the origin of the complex plane.The phase  of the complex International Journal of Optics valued OCT signal i[] is determined in the complex plane with The real and imaginary part of the OCT signal with noise are given by Re{[]} ±  Re{i[]} and Im{[]} ±  Im{i[]} , with the noise variance given by  noise /2.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 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 Hence, the phase variance of the OCT signal is equal to 2SNR −1 using the standard SNR definition of (32). Figure 7 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   = 10 −2 .

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 sides ℓ  = ℓ  = ℓ = 25 m, and depth ℓ  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 ℓ 2 ℓ  , to be distributed randomly over a depth range ℓ  at positions   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  part = 1.4,and volume fraction of  V = 0.01.The refractive index of the medium is  med = 1.33.Based on the volume fraction and particle radius the concentration of particles  is determined with the total number of particles in the simulation given by  part = ℓ 2 ℓ  = 1194.The scattering properties of the particles are calculated using Mie-scattering [26].From the scattering efficiency  the scattering coefficient   =  2  is calculated.The total power scattered   from the first particle is with  in and  in the intensity and power incident on the first particle, respectively.Hence, the power transmission after a single scattering event is  = ( 0 −   )/ 0 .Consequently, the unscattered power remaining after propagation distance ℓ  is the power after transmission through  part particles; that is, Using the relation lim it can be approximated that for many particles (ℓ  ) =  in exp(−  ℓ  ), 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(−  ℓ  ).The reflection of every particle is determined from the collected reflected intensity   from every particle, which is   =  in  NA  2 /ℓ 2 , with  in the power incident on the particle and  NA the fraction of the scattered intensity captured by the collection lens [27].Hence, the intensity reflection coefficient per scattering event is 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 th particle is In the model, every particle has a detected field given by the amplitude of (51) and a delay determined by its position   .A simulation of the OCT intensity |()| 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  2 0  −2   .The height of the numerical OCT signal is with the factor  −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  part /ℓ  , with  the width of a single depth bin.To demonstrate the use of the numerical OCT model, a two-layered piece of tissue is simulated; see Figure 10.Layer 1 is 300 m thick and consists of 500 nm radius particles at 4% volume fraction, with  med = 1.33 and  part = 1.35.Layer 2 is 200 m thick and consists of 1 m particles at 2% volume fraction, with  med = 1.33 and  part = 1.5.

Discussion and Conclusion
In this article I presented a theoretical and numerical analysis of the most important signal processing steps in Fourierdomain OCT.This OCT analysis is based on a comparison of the signals in both the and -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 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 -domain signal using numerical interpolation [30].For nonlinear -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 the 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 -invariant convolution over the wavenumbers, which corresponds to a multiplication of the OCT signal in the -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 -variant convolution.The effect of dispersion on the OCT axial PSF shows that even for a material with only linear dispersion in , 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 -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 -domain is transformed through the inverse Fourier-transformation to the intensity and amplitude of the complex -domain OCT signal.In this analysis it is shown that for an SNR based on the International Journal of Optics 13 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 the -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 backcoupling function.The numerical OCT model for a turbid medium does incorporate the effect of dependent scattering in its dependence on   [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].

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.

2 InternationalFigure 1 :
Figure 1: Schematic representation of the Fourier-domain OCT system.The length of the reference and sample arm is indicated with  and the unidirectional distance from the zero-delay point with .

Figure 2 :
Figure 2: The effect of sampling and detection on the OCT signal.(a) For pixel integration the fringe contrast in the -domain is reduced relative to the input spectrum (blue dashed line) after convolution in the -domain in the continuous (blue continuous line) and sampled signal (red circles).(b) Due to pixel integration the signal in the -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   = 4.(d) Similar to (b) but for the effect of spectral integration.

Figure 3 :
Figure 3: (a) OCT signal acquired with linear  (blue) and nonlinear  (red).(b) OCT axial PSF FWHM for different spectral interpolation algorithms.The lines are a guide to the eye.

Figure 4 :
Figure4: (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).

Figure 5 :
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).

Figure 6 :
Figure6: Schematic of the construction of the OCT phase in the 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 .

Figure 7 :
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   = 10 −2 .

Figure 8 :
Figure 8: Schematic overview of the conversion of the three-dimensional distribution of particles (a) to a one-dimensional distribution in depth (b).

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 Figure9(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 :Figure 10 :
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.

Table 1 :
Parameters used in the OCT simulations.