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.

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 [69], 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.

2. 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. 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 of the optical field is neglected. The incoming optical beam is split by an ideal beamsplitter with (intensity) reflection coefficient and transmission coefficient . 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.

2.1. The Fourier-Domain OCT Signal

The source spectral intensity distribution is defined as a function of wavenumber . The source spectral intensity distribution is normalized such that , where 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 towith 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 bywhere 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 of the sample with 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 (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 asAssuming that the sample is only located at , then combining (4) with (2) and (3) results inAfter some mathematical manipulation the intensity is found to bewith being a symmetric representation of the sample reflectivity. Equation (6) is recognized as the multiplication of the source spectrum with the Fourier transform of for the Fourier pair . Hence, the complex valued depth dependent OCT signal 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 the OCT signal is represented either by the amplitude or intensity 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 isFrom (8) it is determined that the full width at half maximum (FWHM) of the Gaussian point spread function in depth is . From and the relation between the wavenumber and wavelength the iswhich 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 towith being the integration profile. Equation (10) is a convolution .

The detected continuous signal is subsequently sampled at regular intervals and digitized for storage in a computer. The sampled signal is represented as and given bywith 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 isSampling 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 isUsing and linearization, the well-known formula is found.

Considering only the central part of the inverse Fourier transform of (12), the OCT depth scan isThe 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 . The square pixel integration leads to a convolution in the -domain with the filterwith “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 filterThe filter operation in (15) and (16) in the -domain is represented by a multiplication of the OCT signal in the -domain by the functionsrespectively. Note that the normalization with 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].

2.2.1. 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 asThe intensities recorded at wavenumber are transformed to the -domain using the discrete inverse Fourier transform. The spatial coordinate is determined from the linear -domain and spans the range 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 asThe value of is a parameter related to the amount of nonlinearity and can be freely chosen. The parameters and can be calculated1 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.

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, . In case of a dispersive medium (6) changes to the formwith 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 withFor 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 [16] (see Appendix A), the interference spectrum isTo 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 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 isPerforming the inverse Fourier transform to and using the shift in property the OCT signal for isThe OCT signal peak is centered at position ; 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 () and the material is dispersive ( for any ) 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 for (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 functionAfter performing the convolution and some algebra the width of the OCT axial PSF isComparing (27) with derived previously, there is an additional term in the width due to the linear dispersion () 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 , with being the mirror reflectivity, leads to a spectrumTo theoretically determine the shot noise limited SNR an ideal rectangular spectrum 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 isThe factor occurs for both the signal and the noise and can be disregarded. Hence, the peak OCT signal corresponds to a power for a total power received on the detector of .

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 aswith 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 isThe 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 ; 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 . 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.

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 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 lawwith 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 . Consequently, the sample reflectivity for a semi-infinite turbid medium can be modeled asThe 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 . 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 for a semi-infinite turbid medium isThe 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.

3. 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 , and normalization , 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 . 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 . The interferometric signal is constructed according to (4) with the reference arm field given by (2) and the sample arm field given by (3).

3.1. 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 = 1.2 mm. The quasi-continuous -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 . 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 quasi-continuous signal follow the theoretical roll-off. However, when the mirror reflector is placed at a distance larger than , 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).

3.1.1. 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.

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 thickness where it is reflected by an ideal mirror. The refractive index of water is described by the Sellmeier equationThe 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 resulting in = 1.3309 and m rad. 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 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 of the OCT axial PSF in air (). 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 ; 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 the -domain by a square spectrum with channels that is normalized according to ; that is, the power in every channel is . Furthermore, without loss of generality, it is assumed that the sample is a mirror positioned such that the peak of -domain signal is position ; that is, the ratio 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), isPerforming an inverse DFT of (38) results in a fully real signal with two peaks with amplitudesee 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 photonsThe shot noise in every detector channel is modeled as an independent white noise Gaussian distributed random variable with mean and variance , with the number of photons on the 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 towith being the integer valued depth variable. The OCT signal is a complex valued discrete random variable and the OCT signal is calculated by taking either the amplitude or the intensity .

For 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, , and the imaginary part, , of the signal are uncorrelated Gaussian random variables. For a measurement in the absence of signal, the real and imaginary parts of have a variance equal to ; see Appendix D.

The construction of the OCT intensity from the real and imaginary part leads to a variable with a negative exponential distribution [25] with a standard deviationFor an intensity based OCT signal 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 with a Rayleigh distribution [25] with a varianceConsequently, for an OCT SNR defined from the amplitude in terms of the square of the peak amplitude over the noise variance of , the fundamental OCT SNR limit is obtained by combining (38), (40), and (43) and results inwhere 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 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 . Based on the expression 10 10log , 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 and the reference arm reflectivity is varied between 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 . 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 . 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 . Because of the background subtraction the noise is located at the origin of the complex plane. The phase of the complex valued OCT signal is determined in the complex plane withThe real and imaginary part of the OCT signal with noise are given by and , with the noise variance given by . Generalized error analysis, as derived in Appendix E, is applied to calculate the phase uncertainty of (45). The variance of Var() of the angle isUsing 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 asHence, 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 .

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 sides , 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 , 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 , and volume fraction of . The refractive index of the medium is . 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 . The scattering properties of the particles are calculated using Mie-scattering [26]. From the scattering efficiency the scattering coefficient is calculated. The total power scattered from the first particle is , with and the intensity and power incident on the first particle, respectively. Hence, the power transmission after a single scattering event is . Consequently, the unscattered power remaining after propagation distance is the power after transmission through particles; that is,Using the relation it can be approximated that for many particles , in agreement with the single scattering model. Absorption can be treated in a similar fashion by multiplying the transmitted power with an additional factor .

The reflection of every particle is determined from the collected reflected intensity from every particle, which is , with the power incident on the particle and 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 particle isIn 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 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 . The height of the numerical OCT signal is with the factor 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 , with 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.

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 and . Layer 2 is 200 m thick and consists of 1 m particles at 2% volume fraction, with and .

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 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 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 -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 OCT amplitude, the fundamental shot noise limit is a factor 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 back-coupling 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].

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.

Appendix

A. Taylor Expansion of the OCT Signal Phase

Inserting the polynomial expansion of (22) up to quadratic order in the variable ,is obtained. The first term of the Taylor expansion isThe second term of the Taylor expansion isThe third term isHence, the total phase is as used in (26).

B. Peak Value of the OCT Signal

A cosine sampled at points and period points is given by . The inverse DFT of this cosine isThe peak value for positive occurs when the denominator of the second term is zero, that is, when . Using L’Hôpital’s rule it is found that Hence, the signal described by (38) and the peak amplitude of are given by . With a similar approach it can be demonstrated that .

C. Correlation of Real and Imaginary Part of the OCT Signal i

From the definition of the discrete Fourier transform (41), the correlation between and is calculated as follows:which, for yields zero. Hence and are uncorrelated.

D. Variance of Real and Imaginary Part of OCT Signal i

For a zero mean Gaussian random variable with variance , the variance of (and equivalently ) is calculated asUsing the identities and the variance of isFollowing a similar derivation it can be demonstrated that the imaginary part has an identical variance.

E. Variance of the OCT Phase

The phase of the OCT signal in the complex plane iswith the real and imaginary parts of the signal with noise. Defining the argument of the arctan as the variation of the angle is Consider the variable to be a random variable with real and imaginary parts having mean and variance . Then, using standard error propagation [47], isThe variance of the phase is derived asWhen the real and imaginary variances are equal, iswhich is equal to (46).

Competing Interests

The author declares that they have no competing interests.

Acknowledgments

The author thanks M. van Roosmalen, A. K. Trull, J. F. de Boer, and L. J. van Vliet for useful discussions. This research is supported by the Dutch Technology Foundation STW, which is part of the Netherlands Organisation for Scientific Research (NWO) and which is partly funded by the Ministry of Economic Affairs.

Endnotes

  1. It can be shown that, for varying amount of nonlinearity, denoted by , the parameter isand that is