1. Introduction

Dynamic light scattering optical coherence tomography (DLS-OCT) relies on the measurement of fluctuations of scattered light and coherence gating to obtain simultaneous depth-resolved information about diffusive and translational motion of particles. Initially, DLS-OCT was used for measuring diffusion coefficients [1] and flow speeds of particle suspensions [2–5]. Later, several improvements have been suggested for increasing the DLS-OCT flow velocity dynamic range [6,7].

DLS-OCT has the advantage over phase-resolved Doppler OCT that a flow can be measured for zero Doppler angle. It can also be used for particle sizing where the particle size is determined from the estimated diffusion coefficient using the Stokes-Einstein relation. The combined flow and diffusion estimation is particularly interesting for in-line particle sizing during process control [8]. Sensitivity and precision of phase-resolved Doppler OCT has been widely studied and reported in literature [9–12]. However, there is very little information available about the precision and bias of DLS-OCT diffusion and flow measurements. These are crucial for reliable measurements, especially in medical and pharmaceutical applications. Effects of noise and bias in OCT on the measured autocorrelation function have been reported [13], but their influence on the underlying parameters remains unclear.

In this work we perform simulations, measurements, and a theoretical analysis to quantify the precision and bias of diffusion and flow measurements using DLS-OCT. In our analysis, we consider both the first and second-order normalized autocovariance functions, $g_1(z,\tau )$ and $g_2(z,\tau )$, diffusive particle motion, and, for both dilute and non-dilute suspensions, translational particle motion. We derive analytical expressions for the highest attainable precision using the Cramer-Rao bound and compare them with the obtained results. We also assess the bias for different averaging techniques.

2. Theory

The geometry for OCT diffusion and flow measurements is described in [6,7]. The propagation of the optical beam is described by a Gaussian beam along the $z$-direction. We assume that the scattering process is stationary.

2.1 Correlation functions for particle diffusion

For a non-flowing particle suspension, the normalized depth-dependent autocovariance of the OCT complex-valued signal in a backscattering geometry, including the effect of SNR, is given by [3,4,13–15]

The autocorrelation function of the mean-subtracted OCT signal intensity, also known as the Pearson correlation or autocovariance, decays twice as fast than $g_1(z,\tau )$ [1,15] and can be expressed with the normalized second-order autocovariance using the Siegert relation [16,17]

2.2 Correlation functions for particle flow

For diffusing and flowing particle suspensions the first-order normalized autocovariance function magnitude is

For flow measurements we focus on the second-order normalized autocovariance function, $g_2(z, \tau )$, that does not depend on the phase, is easier to implement, and can also be implemented in phase-unstable OCT systems. Here, we differentiate between very dilute and non-dilute sample regimes. For the flowing non-dilute particle suspensions, where the number of particles in the scattering volume, $N_s$, is much greater than one, $g_2(z, \tau )$ is again obtained using the Siegert relation [3,4,6,14,15]

In dilute suspensions, the expected number of particles in the scattering volume can be significantly lower than 1. Therefore, for stationary particles certain depth voxels would give zero signal. However, due to the translational particle motion, the scattering signal is obtained from every voxel during the acquisition time. When the number of particles in the scattering volume is very low, the scattering process becomes non-Gaussian, and the Siegert relationship does not apply anymore [18,21,22]. In the dilute case with $N\lesssim 1$, the second-order normalized autocovariance is [7,13,19]

Here, $A_3(z)$ is given by

2.3 Precision of diffusion and flow estimation

Precision is defined as the spread of the measurement around the mean and can be quantified by the square root of the measurement variance. For DLS-OCT the precision is governed by random errors, i.e., unpredictable changes in the measured intensity. Since the particle diffusion coefficient and the flow speed are determined by fitting the models from Sec. 2.1 and Sec. 2.2 to the measured $g_1(z,\tau )$ and $g_2(z,\tau )$, the precision of the fit parameters is determined by the random errors of the normalized autocovariance function at different time delays. The maximum obtainable precision in the fit parameters is determined by the Cramer-Rao lower bound (CRLB). It is computed by inverting the Fisher information matrix [23]. For calculating the Cramer-Rao lower bound the probability distribution functions of $g_1(z,\tau )$ and $g_2(z,\tau )$ must be known. The correlation coefficients from a long time series data are approximately normally distributed [24–26]. Typical diffusion or flow measurements contain at least thousands of temporal sampling points, which is more than enough to assume that $g_1(z,\tau )$ and $g_2(z,\tau )$ are normally distributed at every time delay. In general, errors (residuals) at different time delays can be correlated. Calculating the CRLB becomes challenging because we don’t know the error correlations and their impact on the model bias. To address this, we assume uncorrelated residuals. In this case, for $M$ observations $g(z,\tau _1),g(z,\tau _2),\ldots g(z,\tau _{M})$ with $N$ model (fit) parameters $p_1,p_2,\ldots p_N$, the Fisher matrix $F$ is an $N\times N$ symmetric matrix given by

The covariance matrix of the model parameters is obtained by inverting the Fisher information matrix. Consequently, the variances of our model parameters, $\sigma ^2_{p_j}$, are represented by the diagonal entries of the covariance matrix. This relationship is described by the equation

2.3.1 Precision of estimating the diffusion coefficient

For a non-flowing particle suspension the diffusion coefficient is determined by fitting Eq. (1) to the real part of $g_1(z,\tau _m)$, with $A_1(z)$ and $D$ being the fit (model) parameters. The lowest attainable variance CRLB in the fitted diffusion coefficient, assuming uncorrelated noise in the Fisher information matrix, is given by

A similar analysis can be performed with $g_2(z,\tau _m)$. In this case we use Eq. (2) for fitting with model parameters $D$ and $A_2(z)$. The Cramer-Rao lower bound on the variance of the diffusion coefficient, when using the second-order normalized autocovariance function, is

Here the higher order (3^rd and 4^th) intensity autocorrelation functions were calculated using theory derived by Lemieux et al. [30]. This analysis is limited to the Gaussian scattering process. Equation (14) incorporates the Siegert relation and relies on the fact that the diffusive $g_1(z,\tau _m)$ is real-valued.

2.3.2 Precision of velocity estimation

The velocity of a non-dilute flowing particle suspension can be obtained by fitting Eq. (4) at every depth $z$ to the measured $g_2(z,\tau _m)$. In this case $v(z)$ and $A_2(z)$ are the model parameters, while $D$, $w_0$, $w_z$, and $\theta$ are assumed to be exactly known a-priori (calibrated). The minimum achievable variance of the flow speed $v$ is given by

For a dilute suspension the velocity can be determined by fitting Eq. (5) to the measured $g_2(z,\tau _m)$. In this case $v(z)$ and $A_3(z)$ are the model parameters, while $D$, $w_0$, $w(z)$, $w_z$, $N_s(z)$ and $\theta$ are assumed to be known in advance. The minimum achievable variance of the flow speed is given by

For comparison, the minimum achievable variance of the flow speed when using $\big |g_1(z,\tau _m)\big |$ from Eq. (3) is independent of the particle concentration and given by

2.4 Bias of diffusion and flow estimation

Bias in DLS-OCT is determined by systematic errors in the correlation function and is a measure of the accuracy of the method. The bias is not random but leads to a consistent over or under estimation of the fit parameter. The systematic errors arise when the model parameters in DLS-OCT are determined by fitting the models from Sec. 2 to the measured $g_1(z,\tau _m)$ and $g_2(z,\tau _m)$. These error sources are:

3. OCT signal simulations for flow and diffusion

In simulations we generate a complex OCT signal scattered from an ensemble of particles with random (diffusion) and directional (flow) motion. The signal intensity and the normalized autocovariance functions are subsequently calculated. This is repeated $N_b$ times which allows us to compute the autocorrelation variances at every time delay.

3.1 Simulation of diffusion

Noiseless time-dependent complex field scattered from $N_p$ diffusing particles is simulated using

3.2 Simulation of flow

The scattered field from suspensions flowing perpendicular to the beam optical axis is simulated by replacing $E_0(t)$ in Eq. (26) with

As we solely model the transverse flow, leading to identical decorrelation regardless of its direction within the transverse plane, we consider only one lateral dimension, denoted as $x$. This allows us to reduce the computational complexity and the number of required particles to be modeled. The particles are simulated with random, uniformly distributed $x$-positions between $-L(z)$ and $+L(z)$ given by

Simulations using Eq. (27) are valid for arbitrary particle concentrations. For non-dilute particle suspensions with $N_s\gg 1$, the normalized flow autocovariance functions from Sec. 2.2 depend only on the beam waist in focus and are independent of $w(z)$ and $R(z)$. In this case Eq. (27) simplifies into

4. Materials and methods

4.1 OCT system

The experiments were performed using a Thorlabs GANYMEDE II HR series spectral domain OCT system, which has been described in detail in our previous works [6,7]. The acquisition rate was 5.5 kHz for diffusion and dilute flow measurements (low-speed), and 36 kHz for non-dilute flow measurements (high-speed). The acquired signal spectrum was measured with a spectrometer with 2048 pixels. The maximum imaging depth in air is 1.87 µm. After acquisition, the measured spectrum was first resampled to a linearly-sampled wavenumber domain and then apodized using a Gaussian filter. After the apodization, the measured coherence function waist in sample was $w_z=2.11$ µm. We have neglected the effect of a gradient of the axial velocity on the autocovariance function for two reasons [20]. First, the Doppler angle in this work is essentially zero $(\theta \approx 0)$. Second, our optical resolution is high both in axial and transverse directions compared to the flow channel dimensions. Hence, the flow velocity within PSF can be assumed to be constant.

The OCT system is operated with a scan lens (LSM04-BB, Thorlabs) in a confocal setup with a manufacturer provided focal spot size of $w_0=6$ µm in air which was validated by axial confocal response measurements [6,7]. The system has $\text {NA}=0.05$. Depending on the angle of incidence, refractive index contrast and Gaussian beam parameters, $w_0$ and $w(z)$ vary somewhat because of the passage of the beam through the various interfaces [35]. Therefore, for flow measurements $w(z)$ and $w_0$ were calibrated using the procedures described in [6,7]. Since for the given OCT setup the coherence length is small and the NA is very low, it can be assumed that the scattering angle is $180^{\circ }$ and the scattering wavenumber $q$ in the correlation analysis is constant at $q=2nk_0$.

4.2 Averaging strategies

We acquire $N_b$ separate OCT interference time traces, $S(k_0,t_m)$, each having a length of $M$, to compute $N_b$ normalized autocovariance functions. Measurements conducted under identical conditions can be averaged. The processing and averaging steps of DLS-OCT are illustrated in Fig. 1. Non-linear curve fitting is utilized, truncating the autocorrelation functions when they reach zero and are fully decorrelated. From $N_b$ measurements, three distinct approaches are used to derive the average model parameters:

For a single time trace measurement the precision is given by the standard deviation $(\sigma _D$ and $\sigma _v)$, but when averaging multiple measurements the precision is given by the standard error of the mean $(\sigma _{\text {average}\:D}$ and $\sigma _{\text {average}\:v})$. Averaging $N_b$ uncorrelated measurements lowers the variance by a factor of $N_b$. For positively correlated measurements, the decrease is slower.

 

Fig. 1. Data processing steps for obtaining $D$ and $v$ from OCT data with different averaging techniques.

Download Full Size | PDF

 

Fig. 2. Schematic overview of the calculation of standard and mixed autocorrelation functions.

Download Full Size | PDF

When averaging the fitted model parameters, $\sigma ^2_D$ and $\sigma ^2_v$ are reduced upon averaging, whereas when averaging the autocorrelation functions, $\sigma ^2_{g_1}$ and $\sigma ^2_{g_2}$ are scaled at all time delays and the noise is reduced. As Eq. (11,13,15,17) show, multiplying $\sigma ^2_{g_1}$ or $\sigma ^2_{g_2}$ by a constant $N_b^{-1}$ at every $\tau _m$ results in the identical scaling of $\sigma ^2_D$ and $\sigma ^2_v$. This holds true for both uncorrelated and correlated random variables. As the averaging sequence does not influence the noise correlation, for unbiased or even slightly biased estimators identical overall precision is expected for the first and second approaches. In the third approach, if the autocorrelation noise at different time delays is correlated, we anticipate that the mixing process will eliminate the correlation in the noise.

Bias cannot be improved by averaging. However, one of the sources of the bias is the estimation bias which, according to Sec. 2.4, exhibits some random nature. Hence, mixing or averaging the correlation functions can reduce this variability. Therefore, we study the bias for all three averaging schemes. The bias is generally very small compared to the measurement uncertainty, and due to a lack of a-priori information it can only be quantified using simulations.

4.3 Diffusion measurements

For diffusion measurements the acquisition rate was 5.5 kHz which is the highest sensitivity setting of our OCT system. All measurements are performed using monodisperse 50 nm radius silica particles with a volume fraction $f_v\approx 1{\% }$, provided by CWK (Chemiewerk Bad Köstritz GmbH, Germany). In total 1100 depth-resolved measurements were performed with a time series length of 4096 points covering $T=0.74$ s. From this data 1100 correlation functions $g_1(z,\tau _m)$ and $g_2(z,\tau _m)$ were calculated at every depth. For diffusion analysis only the real part of $g_1(z,\tau _m)$ was used. Since the OCT sensitivity decreases in depth [36], a usable depth-range of 0.92 mm was chosen where the fitted diffusion coefficients were constant and where the single scattering regime holds. The reference diffusion coefficient was calculated by fitting Eq. (1) to the mean $\Re \big (g_1(z,\tau _m)\big )$ and averaging the resulting diffusion coefficients over the chosen depth range, resulting in $D_0=4.02\pm 0.02$ µm$^2$/s which corresponded to a mean particle radius of 53 nm. The SNR at every depth was calculated based on Eq. (1) using the fitted $A_1(z)$. Simulations were performed using the same parameters as input. In total 10000 simulations were performed for each SNR. Simulated variances were used for calculating the Cramer-Rao lower bounds using Eq. (11) and Eq. (13). Precision as a function of SNR for both diffusion measurements and simulations was determined by computing the standard deviation and the standard error of the mean in the fitted diffusion coefficient both when using $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$. The bias for different averaging schemes was assessed using simulations.

4.4 Flow measurements

The flow was generated using a syringe pump with a variable discharge rate (Fusion 100, Chemyx, Inc.). The flow geometry was aligned so that $\theta \approx 0$. This minimized the velocity gradient effect on the correlation [20] and simplified the beam waist calibration. Other than that, the effect of a nonzero Doppler angle on DLS-OCT is minimal because of similar optical resolutions in axial and transverse directions.

4.4.1 Non-dilute flow measurements

For non-dilute flow measurements we used 20 mL syringe (Terumo Europe NV) and OCT acquisition rate of 36 kHz. Diluted Intralipid solutions were used as a sample with a particle volume fraction $f_v \approx 0.6{\% }$. The flow passes through a quartz rectangular flow cell with internal dimensions of 0.2 mm depth and 10 mm width (type 45-F, Starna Scientific). M-scan 1D measurements were performed at the center width of the flow cell as a function of depth with discharge rates of 1.5 and 2 mL/min. The beam focus was moved away from the center depth of the sample to increase the number of particles in the scattering volume and reduce the effect of number fluctuations [18,19]. The particle diffusion coefficient and the beam waist in focus were calibrated using a $g_2(z,\tau _m)$ based on lateral scanning over the static sample [6] and were found to be $D=1.40\pm 0.09$ µm$^2$/s and $w_0=8.08\pm 0.31$ µm, respectively. The flow speed was determined by fitting Eq. (3,4) to the measured or simulated first and second-order autocovariance functions. The depth-dependent SNR was calculated from the fitted autocovariance magnitude. In total 1000 measurements and 10000 simulations were performed with a time series length of 4096. The average depth-dependent flow velocities and SNR values from measurements were used as the input for simulations. Simulated variances were used for calculating the Cramer-Rao lower bounds both for $ |g_1(z,\tau _m) |$ and $g_2(z,\tau _m)$. The bias for different averaging schemes was also determined using simulations.

4.4.2 Dilute flow measurements

Dilute flow measurements were performed with the second-order normalized autocovariance function using the number fluctuations analysis. The OCT acquisition rate was 5.5 kHz and a 5 mL syringe (BD Plastipak) was used. Experiments were performed using a monodisperse suspension of polystyrene particles with volume fraction $f_v\approx 0.006{\% }$ and expected particle radius of $230-250$ nm, provided by InProcess-LSP. The flow passes through a quartz rectangular flow cell with internal dimensions of 0.5 mm depth and 10 mm width (type 45-F, Starna Scientific). M-scan 1D measurements were performed at the center width of the flow cell as a function of depth at discharge rates of 0.15 and 0.225 mL/min. In contrast to our previous work, where only the number fluctuation term was used [7], here we use the full model from Eq. (5).

The particle diffusion coefficient was determined (calibrated) in the stationary suspension using $\Re \big (g_1(z,\tau _m)\big )$ and Eq. (1). The obtained diffusion coefficient of $D_0=1.00\pm 0.03$ µm$^2$/s corresponded to an average particle radius of 220 nm. The beam shape calibration procedure was identical to that in [7], but the analysis was slightly different because of using the full autocorrelation model. The SNR at every depth was calculated by fitting Eq. (3) to the measured $\big |g_1(z,\tau _m)\big |$. The kurtosis was determined using Eq. (8). Then $A_3(z)$ was found based on Eq. (7) using the known SNR and kurtosis. Finally, Eq. (6) was fitted to the measured $g_2(z,\tau _m)$ with $N_s(z)$ and $w(z)$ being the fit parameters and using $A_3(z)$, SNR, and the measured $\big |g_1(z,\tau _m)\big |$. The obtained beam shape $w(z)$ was fitted using $w(z)=w_0\sqrt {1+(z-z_0)^2/z^2_R}$, with $w_0$, $z_0$, and $z_R$ being the free parameters equaling $6.45\pm 0.02$ µm, $187\pm 3$ µm and $293\pm 3$ µm, respectively. The particle volume fraction $f_v$ was computed using $w(z)$, $D_0$, $w_z$ and $N_s(z)$ according to [7]. It is given in Fig. 3(a,b) along with the obtained $w(z)$.

 

Fig. 3. Obtained (a) beam shape $w(z)$ and (b) particle volume fraction $f_v$ as a function of depth.

Download Full Size | PDF

The obtained $w_z(z)$ and $N_s(z)$ were used in subsequent number fluctuation flow measurements. In total, 750 measurements were performed each with a time series length of 8192. The flow speed was determined by fitting Eq. (5) to the measured second-order autocovariance function with $A_3(z)$ and $v(z)$ being the free parameters. Simulations for the dilute regime were not performed due to the adverse effects of boundary conditions and long computational times. Therefore, the measured variances were used for calculating the Cramer-Rao bounds. Even though number fluctuations are only present in $g_2(z,\tau _m)$, for comparison the CRLB for $ |g_1(z,\tau _m) |$ was also determined.

5. Results

5.1 DLS-OCT diffusion measurement

Depth dependent diffusion measurements were performed on a static particle suspension. Examples of measured $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ for $\text {SNR}=100$ are shown in Fig. 4(a). Observe here that the autocorrelation functions do not oscillate around the mean, but, instead, consistently remain above or below it for some time. From $N_b$ measured normalized autocovariance curves (signal ACF), the normalized autocovariances of errors (error ACF) in the measured and simulated $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ were calculated using

 

Fig. 4. Example (a) standard and (b) mixed autocorrelation functions (signal ACF) for $\text {SNR}=100$. (c,d) Measured (points) and simulated (lines) error ACF, $\rho _e(z,\tau _e)$.

Download Full Size | PDF

Figure 4(c,d) illustrates that in standard autocorrelation functions, errors at different time delays exhibit strong correlation. This means that the magnitude and direction of errors in $\Re \big (g_1(z,\tau _m)\big )$ or $g_2(z,\tau _m)$ at small $\tau _m$ significantly affect the errors at larger $\tau _m$, which can be observed from the fact that any single autocorrelation function is consistently above or below the mean autocorrelation function. As described in Sec. 4.2, we implement a novel way of reducing these correlations or even completely eliminating them by mixing different autocorrelation functions at every $\tau _m$. The number of independent autocorrelation functions used for mixing is denoted by $N_\text {mix}$. To fully get rid of the correlations $N_\text {mix}$ must be greater or equal to the number of sampling points it takes for $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ to go to zero. For $g_2(z,\tau _m)$ it takes fewer realizations to eliminate correlations because it decays twice as fast.

Figure 4(b-d) show that for a single fully mixed $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ with sufficiently large $N_\text {mix}$ the errors at different $\tau _m$ are completely uncorrelated. As a result, the measured mixed autocorrelation data oscillate approximately around the mean as a function of $\tau _m$. With standard (unmixed) measurements, as shown in Fig. 4(a), error ACF is nonzero and every measured $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ deviates from the mean autocorrelation function. Figure 4(c,d) show that once the autocorrelation functions are sufficiently mixed (e.g. $N_\text {mix}=32$), error ACF becomes a delta function. At intermediate mixing ratios $(N_\text {mix}=4)$ error ACF becomes periodic as the correlations reappear at later time points.

5.1.1 Precision in diffusion estimation

Next, we look at the diffusion estimation precision from a single autocorrelation measurement. To estimate the precision, the autocorrelation variance at every time delay needs to be known. Figure 5(a,b) show measured, simulated, and analytical variances from Eq. (12), (14) for the diffusive $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ at different SNR values. All three match relatively well with each other. The variances are higher for larger SNR values, but, in this case the correlation coefficients are also larger. Even though the variances of $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ are comparable at larger time delays, $\sigma ^2_{g_2}(z,\tau _m)$ is significantly larger at small $\tau _m$. The different shape of $\sigma ^2_{g_2}(z,\tau _m)$ and its sharp increase at small $\tau _m$ is due to the mean subtraction from the intensity time series.

 

Fig. 5. Measured, simulated, and analytical variances of diffusion only (a) $g_1(z,\tau _m)$ and (b) $g_2(z,\tau _m)$.

Download Full Size | PDF

Figure 6(a,b) show measured and simulated standard deviations in the fitted diffusion coefficient from a single correlation function, as well as the Cramer-Rao lower bound calculated with Eq. (11) and (13, where the $z$-dependence is converted to SNR. The results are based on 10000 simulations and 1100 measurements. Here $D_0$ was calculated from the measured $\Re \big (g_1(z,\tau _m)\big )$ as described in Sec. 4.3 and subsequently used as input for simulations. The obtained $\sigma _D$ strongly depends on whether the errors in our normalized autocovariance functions are correlated. The obtained $\sigma _D$ from the standard measurements is several factors larger than the Cramer-Rao lower bound. After mixing the resulting $\sigma _D$ matches very well with the Cramer-Rao bound and we achieve the most precise determination of the diffusion coefficient possible. The measured and simulated $\sigma _D$ overlap each other except for $g_2(z,\tau _m)$ at very low SNR values. Here the simulations underestimates the error. The lowest spread in the fitted $D$, given by $\sigma _D$, is obtained when mixing the normalized autocovariance functions for removing the error correlations. However, this requires many measurements, as shown in Fig. 6(a,b), and is therefore less relevant for practical applications with few measurements. Typically, the number of DLS-OCT measurements is on the order of 5-10.

 

Fig. 6. Measured (points) and simulated (lines) standard deviation for estimating the diffusion coefficient from a single autocorrelation function and the corresponding Cramer-Rao lower bound for (a) $g_1(z,\tau _m)$ and (b) $g_2(z,\tau _m)$.

Download Full Size | PDF

Figure 7(a,b) show the standard deviation (spread) in the diffusion coefficient obtained from fitting a single correlation function as a function of the number of measurements $N_b$. For conventional DLS-OCT measurements without mixing, $\sigma _D$ is constant as expected for independent measurements. When employing autocorrelation mixing we see that $\sigma _D$ decreases with increasing $N_b$ until reaching the CRLB and after which it remains constant with $N_b$. For $g_2(z,\tau _m)$ it takes fewer measurements for $\sigma _D$ to reach the CRLB.

 

Fig. 7. (a,b) Standard deviation and (c,d) standard error of the fitted diffusion coefficient as a function of the number of averaged correlation functions for $\text {SNR}=100$. Points correspond to DLS-OCT measurements and lines denote simulations.

Download Full Size | PDF

In DLS-OCT measurements are averaged to improve the precision. In this case the measurement precision is given by the standard error of the mean, $\sigma _{\text {average}\:D}$. Figure 7(c) shows simulated and measured standard error of the mean diffusion coefficient, $\sigma _{\text {average}\:D}$, as a function of $N_b$ for different averaging techniques. The standard error of the mean is identical for all three averaging methods, even though $\sigma _D$ itself is different depending on error correlations. The decrease is proportional to ${N_b}^{-1/2}$. Figure 7(d) shows how $\sigma _{\text {average}\:D}$ is related to $\sigma _{D}$ for the three different averaging techniques. For the standard analysis we see that $\sigma _{\text {average}\:D}=\frac {\sigma _{D}}{\sqrt {N_b}}$ for both averaging methods. However, when mixing the normalized autocovariance functions, we notice that $\sigma _{\text {average}\:D}>\frac {\sigma _{D}}{\sqrt {N_b}}$ as the fitted diffusion coefficients now become statistically interdependent. Here $\sigma _{\text {average}\:D}$ initially decreases very slowly with $N_b$ and then starts to reduce faster above a certain $N_b$. This coincides with the $N_b$ at which the CRLB is reached in Fig. 7(a,b). In Fig. 7(d) there is a small discrepancy between measurements and simulations at large $N_b$ which caused by the lack of sufficient number of averaged measurements.

5.1.2 Bias in diffusion estimation

The bias in the average diffusion coefficient is determined by comparing the averaged $D$ to a ground truth. We average 1100 and 10000 values of $D$ in measurements and simulations, respectively. With this large set of averages we can assume that the random errors are negligible and only the systematic errors remain. Figure 8(a,b) shows the measured diffusion coefficients as a function of SNR for both $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ with different averaging schemes. Figure 8(c,d) show the same results but obtained using simulations. In this case the input $D_0$ is also displayed. Both in measurements and simulations we can observe that mixing different correlation functions slightly reduces the bias. However, the lowest bias is obtained for $D_{\bar {g}}$ when averaging the correlation functions. This is more evident for $g_2(z,\tau _m)$. The differences for the three techniques are minimal for $\Re \big (g_1(z,\tau _m)\big )$. Overall, the bias is much smaller and less affected by SNR when using $\Re \big (g_1(z,\tau _m)\big )$ compared to $g_2(z,\tau _m)$.

 

Fig. 7. Diffusion coefficients as a function of SNR and averaging scheme determined using (a,b) measurements and (c,d) simulations.

Download Full Size | PDF

The trends from both measurements and simulations in Fig. 8(a,b) and Fig. 8(c,d) are similar, however for the measurements we lack a-priori knowledge of the ground truth $D_0$. Therefore, it is not possible to compare the bias directly between measurements and simulations. This can be rectified on a relative basis by not comparing the absolute deviation from $D_0$, but the relative bias with respect to the most accurate method (correlation averaging). Figure 9(a,b) displays the relative bias when averaging the fitted diffusion coefficients with respect to averaging $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ first and then fitting the diffusion coefficient. The obtained measurements and simulations are in perfect agreement except for $g_2(z,\tau _m)$ at very low SNR values. This emphasizes the reliability of our simulations and that they can be used for determining the absolute bias.

 

Fig. 9. Measured (points) and simulated (lines) relative bias in the diffusion coefficient as a function of SNR for (a) $\Re \big (g_1(z,\tau _m)\big )$ and (b) $g_2(z,\tau _m)$.

Download Full Size | PDF

The bias of the mean diffusion coefficient from Fig. 8(c,d) is determined by averaging over a large number of simulations. Figure 10(a,b) show simulated systematic errors as a function of the number of averaged measurements $N_b$ for $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$. For $\Re \big (g_1(z,\tau _m)\big )$ the bias is very low and almost independent of $N_b$, as expected. The bias is more significant for $g_2(z,\tau _m)$. Figure 10(b) shows the estimation bias due to mean subtraction modeled using Eq. (22). This is the main source of bias for $g_2(z,\tau _m)$ which is absent in $\Re \big (g_1(z,\tau _m)\big )$. The estimation bias also contains a random component which averages out with increasing number of measurements. As a result, the bias in $g_2(z,\tau )$ initially decreases with $N_b$ until reaching the minimum and then remains constant. Overall, the bias is minimal when averaging the autocorrelation functions. The second best technique is mixing the autocorrelation functions before the curve fitting and averaging the fitted diffusion coefficients. The worst approach is to the fit model parameters to the standard autocorrelation functions and then average them.

 

Fig. 10. Simulated bias of the fitted diffusion coefficient as a function of the number of used correlation functions at $\text {SNR}=100$ for (a) $\Re \big (g_1(z,\tau _m)\big )$ and (b) $g_2(z,\tau _m)$.

Download Full Size | PDF

5.2 Non-dilute DLS-OCT flow measurement

Precision analysis of flow was performed for DLS-OCT $\big |g_1(z,\tau _m)\big |$ and $g_2(z,\tau _m)$ measurements of a laminar transverse flow corrected for diffusion. Figure 11(a,b) show measured and simulated variances at different flow speeds in the channel, good agreement between the two is observed. The variances are higher with lower flow speeds, but do not depend much on the SNR. The highest SNR value in the channel was above 100, and the lowest was 6. This is sufficiently high to neglect the SNR dependence of $\sigma _v$ for plotting purposes and display it only as a function of the flow speed. Velocity profiles, shown in Fig. 12(a), were obtained by averaging the measured correlation functions and were subsequently used as ground-truth input for simulations.

 

Fig. 11. Measured and simulated variances for flow (with a fixed and calibrated diffusion coefficient) of (a) $\big |g_1(z,\tau _m)\big |$ and (b) $g_2(z,\tau _m)$.

Download Full Size | PDF

 

Fig. 12. Non-dilute flow measurements. (a) Measured flow profiles. (b) Simulated bias in the fitted velocity. Measured and simulated $\sigma _v$ using (c) $\big |g_1(z,\tau _m)\big |$ and (d) $g_2(z,\tau _m)$ and the corresponding CRLB.

Download Full Size | PDF

Figure 12(c,d) show measured and simulated standard deviation in the fitted flow speed as a function of velocity. Results are obtained using both the standard method with correlated errors and the mixing technique with uncorrelated errors are given. Similar to diffusion, mixing decreases the spread in the obtained velocity. For comparison, the Cramer-Rao lower bounds for $g_2(z,\tau _m)$ and $ |g_1(z, \tau ) |$, calculated using Eq. (15) and (21), are also displayed. The theoretical Cramer-Rao bound for $ |g_1(z, \tau ) |$ is lower than for $g_2(z,\tau _m)$. For both methods there is a good agreement between measurements and simulations. The obtained $\sigma _v$ from measurements and simulations is significantly improved when mixing removes the error correlations and matches well with the Cramer-Rao bound, especially for $g_2(z,\tau _m)$. In general, the standard deviation in the flow speed decreases with lower velocities. However, it does not decrease to zero and has a certain minimum value. Below the threshold velocity $\sigma _v$ starts to increase because of diffusion. This is only clearly visible when the errors are not correlated.

Figure 12(b) displays the simulated bias for different averaging schemes both for $ |g_1(z, \tau ) |$ and $g_2(z,\tau _m)$. The bias was calculated with respect to the simulation input which included both depth-dependent velocity and SNR profiles. The simulated bias is lowest when averaging the complex $g_1(z,\tau _m)$, taking the absolute value and then fitting $v(z)$. This is denoted by the black curve. The second lowest bias is obtained when averaging $g_2(z,\tau _m)$ and then fitting $v(z)$ or when mixing $g_2(z,\tau _m)$, fitting and averaging $v(z)$. In this case the errors are identical and given by the red curve. This approach is only marginally less accurate than the former method. The third best method relies on using the standard $g_2(z,\tau _m)$, fitting and then averaging $v(z)$. This approach is less accurate than the second best method only at low flow speeds. The largest bias is obtained when using $ |g_1(z, \tau ) |$ to fit $v(z)$ and then average the obtained velocities. In this case the sample population is significantly skewed and the mixing method has virtually no effect.

5.3 Dilute DLS-OCT flow measurement

Similar to non-dilute, we analyse precision of dilute DLS-OCT flow measurements where the effect of diffusion is absent in $g_2(z,\tau _m)$. Also for this case we obtained the variance of the correlation function for different flow speeds. The results look similar to the non-dilute case, except there is no overshoot at small $\tau _m$. Figure 13(a,b) show dilute flow measurements using the second-order normalized autocovariance function incorporating the number fluctuations. Velocity profiles are shown in Fig. 13(a) which were again obtained by averaging the measured correlation functions. The lowest SNR value in the channel was again around 6, and the highest was above 100. In Fig. 13 we have again neglected the SNR dependence for the same reasons as mentioned in Sec. 5.2 and plot $\sigma _v$ only as a function of the velocity.

 

Fig. 13. Dilute flow measurements. (a) Measured flow profiles using $g_2(z,\tau _m)$. (b) Measured $\sigma _v$ using $g_2(z,\tau _m)$ and the corresponding CRLB.

Download Full Size | PDF

Figure 13(b) shows the measured standard deviation in the fitted flow speed as a function of velocity. Results using both the standard $g_2(z,\tau _m)$ with correlated errors and the mixed $g_2(z,\tau _m)$ with uncorrelated errors are given. Measurements using $\big |g_1(z,\tau _m)\big |$ are not shown because they are too noisy and have a significantly larger $\sigma _v$. For comparison, the Cramer-Rao lower bounds for $g_2(z,\tau _m)$, calculated with Eq. (17), and $ |g_1(z, \tau _m) |$, calculated with Eq. (21), are also displayed. Here we used the measured variances of $\sigma ^2_{g_2}(z,\tau _m)$ and $\sigma ^2_{|g_1|}(z,\tau _m)$ for calculating the Cramer-Rao lower bounds. The Cramer-Rao bound for $g_2(z,\tau _m)$ is considerable lower than for $ |g_1(z, \tau _m) |$ and agrees very well with our measurements. When using $g_2(z,\tau _m)$ the velocity standard deviation decreases, with decreasing velocity, to zero. However, this is not the case with $\big |g_1(z,\tau _m)\big |$. As the Cramer-Rao lower bound shows, $\sigma _v$ increases dramatically at lower flower speeds when using $ |g_1(z, \tau _m) |$. This happens because $\big |g_1(z,\tau _m)\big |$ in Eq. (3) is ultimately limited by particle diffusion and is independent of number fluctuations. The second-order normalized autocovariance function in Eq. (5), on the other hand, contains the number fluctuation term that is independent $D$. This increases the velocity precision for this method at extremely low flow speeds. For number fluctuation DLS-OCT no bias analysis is performed because of boundary condition limitations and increasing computational complexity of our simulations.

6. Discussion

In this work we review autocorrelation averaging and mixing techniques for reducing random and systematic errors. In the absence of an analytical model for error correlations, only using the mixing method suggested by us, it is possible to quantify the precision in DLS-OCT and verify whether the CRLB is reached. In diffusion measurements using the average real part of $g_1(z,\tau _m)$ results in the highest accuracy and precision. In flow measurements it is more reliable and convenient to use the average $g_2(z,\tau _m)$ because it is not affected by the sampling distribution bias, does not depend on phase and can be implemented in phase-unstable systems. It is important to note that mixing and/or averaging of the autocorrelation functions can only be performed using the autocorrelation functions with the same SNR (and identical optical or sample properties). For single scattering diffusion measurements where $D$ is constant as a function of depth (SNR), we can also average the measurements with different signal-to-noise ratios. In this case the measured $g_1(z,\tau _m)$ and $g_2(z,\tau _m)$ first have to be noise-corrected [13] before any further processing.

6.1 DLS-OCT diffusion estimation

Random error correlations in $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ are the main reason for $\sigma _D$ not reaching the theoretical Cramer-Rao bound. Once these correlations are removed by mixing, the standard deviation in the fitted diffusion coefficient is reduced significantly and reaches the CRLB. The improvement is bigger for larger SNR values since $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ are lower in the presence of more noise and hence the random errors are less correlated. The improvement with mixing is less for $g_2(z,\tau _m)$ because it decays faster and therefore correlations play less of a role in the precision. We can conclude that faster decorrelations and smaller autocovariance magnitudes, i.e., more noisy measurements, are less affected by mixing. We have not derived analytical expressions for $\rho _e$, but our observations suggest that it depends on the autocorrelation decay rate and also experiences a delta function noise decorrelation at $\tau _m=0$.

Even though the improvement in $\sigma _D$ due to mixing is significant, practical applications of the mixing technique for increasing the measurement precision are limited. In order to reduce some of the error correlations we need at least two measurements. Mixing just 2 measurements already provides us with a significant improvement in $\sigma _D$. However, the precision of the average diffusion coefficient when using multiple measurements is given by the standard error of the mean, $\sigma _{\text {average}\:D}$, which is identical for all averaging techniques with or without mixing the correlation functions. So, even though $\sigma _D$ is significantly lower when mixing the autocorrelation functions, it does not decrease as fast with averaging because the mixed autocorrelation functions are not statistically independent anymore. As a result, the fitted diffusion coefficients are also interdependent and $\sigma _{\text {average}\:D}\neq \frac {\sigma _{D}}{\sqrt {N_b}}$. When $\Re \big (g_1(z,\tau _m)\big )$ or $g_2(z,\tau _m)$ are not mixed, $\sigma _{\text {average}\:D}$ is proportional to $N_b^{-1/2}$ because every fitted $D$ is independent. With mixing we basically move the error correlations from $\Re \big (g_1(z,\tau _m)\big )$ and $g_2(z,\tau _m)$ to the fitted $D$. Therefore, the actual measurement precision is identical for all averaging methods and cannot be improved any further.

We noticed that systematic errors are much larger when using $g_2(z,\tau _m)$ compared to $\Re \big (g_1(z,\tau _m)\big )$ and are dominated by the estimation bias. Averaging $N_b$ measurements slightly reduces the bias in $g_2(z,\tau _m)$ but only to a certain limit, beyond which the errors remain constant. This is probably caused by variability and randomness in the estimation bias which averages out and diminishes with increasing $N_b$. We have mentioned that in phase-stable systems it is always preferable to utilize $\Re (g_1(z,\tau _m) )$ instead of $ |g_1(z, \tau ) |$. Using $ |g_1(z, \tau ) |$ adversely affects the bias because the autocorrelation coefficients are always positive. This leads to a higher sampling distribution bias which is otherwise negligible. Overall, random errors are larger than systematic errors unless we average thousands of measurements.

There is an asymptotic dependence of precision and bias on $\text {SNR}$. Therefore, the benefits of using extremely high signal-to-noise ratios are rather limited. As Fig. 6(b) and Fig. 9(b) show, there is a mismatch between simulations and measurements at very low SNR when using $g_2(z,\tau _m)$. This is caused by the specific noise model used in our simulations [33] that is based on the assumption that the signal intensity is much larger than the noise intensity. However, this assumption is not valid at low SNR values. As a result, this noise model overestimates $A_2(z)$ at low SNR values, leading to the lower $\sigma _D$.

6.2 DLS-OCT flow estimation

Both in dilute and non-dilute flow measurements, the standard deviation (spread) in the fitted velocity is significantly reduced when the error correlations are removed. In non-dilute flows, $\sigma _v$ decreases with decreasing velocity and then starts to increase at very low flow speeds. A minimum in the error occurs because the diffusive limit is reached. At this stage, any further reduction in the velocity does not affect autocorrelation functions or their variances. According to Eq. (15), the velocity is explicitly included in the denominator of $\sigma ^2_{\text {CRLB},\,v}$. As a result, $\sigma _v$ starts to increase when the diffusive limit is reached. This behaviour is overshadowed by error correlations and is only visible when employing correlation function mixing.

For non-dilute flow the Cramer-Rao bound for $g_2(z,\tau )$ is slightly lower than the simulated $\sigma _v$ at larger velocities, which itself is marginally lower than the measured $\sigma _v$. The difference between simulations and the Cramer-Rao bound is probably because of the Siegert approximation. In this case, the approximate number of particles in the scattering volume is around 49 both in simulations and measurements. This could be sufficiently low to violate the large number of particle assumption and have a small effect on $g_2(z,\tau _m)$. Furthermore, the offset between measurements and simulations can be caused by galvo or pump instabilities. The Cramer-Rao bound of $ |g_1(z, \tau ) |$ shows that on paper it is possible to achieve higher precision in non-dilute flow measurements when using $g_1(z, \tau )$ instead of $g_2(z, \tau )$. However, in practice this is more problematic. The mismatch among measurements, simulations and the Cramer-Rao bound at all velocities is considerably larger for $ |g_1(z, \tau ) |$ compared to $g_2(z, \tau )$. This is largely caused by the bias of the sampling distribution and the limited phase stability.

Theoretically, the lowest bias in non-dilute flow measurements can also be achieved when using $ |g_1(z, \tau ) |$ by averaging the complex autocorrelation functions first, then taking the absolute value, and then fitting $v(z)$. However, getting rid of the sampling distribution bias when taking the absolute value requires sufficient averaging of the complex $g_1(z,\tau _m)$. In our work we averaged 1000 autocorrelation functions which in real-time flow measurements is impossible. Insufficient averaging will result in oscillations with $\big |g_1(z,\tau _m)\big |>0$. In this work we have not investigated the dependence of the sampling distribution bias on the number of measurements. However, it is clear that with fewer measurements it is preferred to use $g_2(z,\tau )$ as it does not suffer from the population skewness.

Similar to diffusion, random velocity errors in flow measurements are also generally larger than systematic errors. However, at very low speeds both errors become comparable. Despite some advantages of $\big |g_1(z,\tau _m)\big |$, for practical reasons it is more convenient to use $g_2(z,\tau _m)$ in non-dilute flows. In dilute low-speed flows it is always preferable to use $g_2(z,\tau _m)$ because it is not limited by the particle diffusion.

7. Conclusion

We have investigated precision and bias of diffusion coefficient and flow speed measurements using DLS-OCT based on the first and second-order normalized autocovariance functions. We found that errors in the autocovariance functions are strongly correlated. This significantly reduces the precision and bias of fit parameters and prevents us from reaching the Cramer-Rao lower bound. We demonstrated that mixing different autocovariance functions at every time delay before the curve fitting procedure reduces the standard deviation in the fitted parameters and reaches the Cramer-Rao lower bound. This proves the validity of our DLS-OCT models. When using the mixing technique the correlations are transferred to the fitted parameters and the standard error of the mean remains unchanged. We conclude that the precision in DLS-OCT is identical for all averaging techniques and ultimately limited by correlations between the random variables. The bias is lowest when averaging the measured normalized autocovariance functions before fitting the model parameters.