Quantitative technique for robust and noise-tolerant speed measurements based on speckle decorrelation in optical coherence tomography.

Intensity-based techniques in optical coherence tomography (OCT), such as those based on speckle decorrelation, have attracted great interest for biomedical and industrial applications requiring speed or flow information. In this work we present a rigorous analysis of the effects of noise on speckle decorrelation, demonstrate that these effects frustrate accurate speed quantitation, and propose new techniques that achieve quantitative and repeatable measurements. First, we derive the effect of background noise on the speckle autocorrelation function, finding two detrimental effects of noise. We propose a new autocorrelation function that is immune to the main effect of background noise and permits quantitative measurements at high and moderate signal-to-noise ratios. At the same time, this autocorrelation function is able to provide motion contrast information that accurately identifies areas with movement, similar to speckle variance techniques. In order to extend the SNR range, we quantify and model the second effect of background noise on the autocorrelation function through a calibration. By obtaining an explicit expression for the decorrelation time as a function of speed and diffusion, we show how to use our autocorrelation function and noise calibration to measure a flowing liquid. We obtain accurate results, which are validated by Doppler OCT, and demonstrate a very high dynamic range (> 600 mm/s) compared to that of Doppler OCT (±25 mm/s). We also derive the behavior for low flows, and show that there is an inherent non-linearity in speed measurements in the presence of diffusion due to statistical fluctuations of speckle. Our technique allows quantitative and robust measurements of speeds using OCT, and this work delimits precisely the conditions in which it is accurate.


Introduction
Optical coherence tomography (OCT) is a technique that provides high-resolution imaging with moderate penetration depths. It is commonly used in biomedical and clinical applications due to its combination of optical sectioning, axial resolution and high-speed imaging [1,2]. OCT techniques have also been applied to optical inspection for non-destructive industrial applications [3,4,5,6]. Signal and processing methods, such as polarization-sensitive OCT [7], spectroscopic OCT [8], OCT-based elastrography [9] and vibrography [10], have additionally been developed to provide functional and molecular contrast.
A valuable piece of information that can be acquired with OCT techniques corresponds to motion or deformation of sample. In OCT there have been two approaches traditionally applied for motion sensing: Doppler-based and speckle-based. As the tomogram provided by OCT measurements contains information about the reflected complex amplitude, Doppler OCT utilizes the phase of the reflected light wave as a function of time to derive information related to motion of the sample [11]. Doppler has been the technique of choice when quantitative measurements are important [12], but the stringent requirements on phase-stable detection makes its deployment in fast frequency-domain OCT systems difficult [13]. Additionally, robust Doppler methods are sensitive only to motion along the axis of the optical beam [11]. In comparison, speckle techniques rely on the relation between the statistics of the intensity signal and the movement of scatterers within the sample. Speckle techniques have primarily been applied to qualitative measurements, for example to characterize the morphology of vascular networks in biomedical applications [14,15,16]. Few studies have addressed quantitation from speckle decorrelation measurements [17,18,19], and we are unaware of studies validating quantitative speckle analysis. A recent report acknowledged this lack of validation and discussed the work needed to bring a particular implementation of speckle decorrelation, split-spectrum angiography, into the realm of quantitative flow measurements [20]. Despite their advances, however, the authors recognized that more insight was needed in order to realize this goal. Dynamic light scattering techniques have been used in OCT to extract speed, axial velocity and diffusion information from the complex amplitude of the OCT signal [21,22,23,24]. Considering the literature on ultrasound speckle decorrelation for flow speed measurements in medical applications, it has been found that velocity estimations are strongly dependent on the noise of the data [25,26], a point that has not been addressed in speckle OCT systematically. Interestingly, we are aware of only a recent study that approaches the issue of a thorough validation of ultrasound speckle decorrelation measurements [27].
This work aims to overcome the current limitations of speckle-decorrelation speed measurements. First we analyze the effect of background noise on the intensity autocorrelation function, and define a new autocorrelation function that is immune to the primary effect of additive noise. This autocorrelation function also provides flow contrast information, which is lacking in previous speckle-decorrelation approaches. We then extend the signal-to-noise ratio (SNR) range in which speckle decorrelation measurements work by performing a calibration of the secondary effect of background noise using a solid phantom. In order to show the validity of this technique beyond solid samples, we derive the explicit relation between decorrelation time, flow speed and diffusion, and make measurements in a flow phantom setup and compare the results with Doppler OCT. In order to have Doppler OCT measurements with high accuracy to be used as references, we implemented a posteriori correction of phase [13]. Finally, we explain the effect of statistical noise on the measurement of low flow speeds, and how the non-linear behavior of speckle decorrelation is inherent to the statistical nature of the technique.

Materials and methods
The experimental part of this work used an OCT system with a wavelength-swept laser source having a 134 nm sweep range, centered at 1285 nm and a repetition rate of 102 kHz (one Aline every ∆t = 0.0098 ms). The signal was digitized at 180 MS/s, recording 1600 samples per sweep. The system utilized an acousto-optical frequency-shifter for depth-degeneracy removal [28] and the tomograms were reconstructed using zero-padding to have a reconstruction size of 1024 pixels in the axial direction. The system has a polarization-diverse receiver (described in [29]), and the signal from each channel was added to calculate the tomogram intensity I = |E x | 2 + |E y | 2 . We define the tomogram amplitude as |E| = √ I. Solid samples consisted of white rubber and liquid samples of intralipid at 0.5% volume in water. The waist radius of the lateral point spread function (PSF) of the optical probe was approximately w n = 15-18 µm in water.
We use a coordinate system defined by assigning the z axis to the light propagation direction, as shown in Fig. 1. The solid sample was placed on a motorized translation stage that moved the sample at a known speed in the x direction.
Liquid phantom measurements were performed using the geometry shown in the inset of Fig. 1. The flow setup consisted of a peristaltic pump, a closed 3.2 mm-diameter tubing circuit and a pulse dampener in order to have steady, laminar flow in the area of the OCT probe, which was used to take measurements from within the tube. The probe emitted light at an angle of 82.5 • with respect to the flow direction, which allowed for a small Doppler signal to be acquired from the liquid. The Doppler signal was analyzed using the Kasai autocorrelation [11] in the same correlation window that was used for decorrelation measurements.
To obtain a clear Doppler signal from noisy phase measurements, the OCT system was setup in order to have two calibration reflection signals to compensate for the phase-unstable acquisition, in a similar way to the technique described in [13]. The first reflection signal, close to zero-depth, was used to compensate for the phase offset introduced by the acousto-optical frequency-shifter used for depth-degeneracy removal [28]. The second was used to compensate for the phase slope introduced by the phase jitter of the data acquisition board.

Time autocorrelation of OCT signals
The theory of dynamic light scattering applied to the OCT signal has been studied thoroughly previously [21], and it is well known that the time autocorrelation of the complex OCT signal z x y rubber phantom translation stage intralipid optical fiber connected to OFDI system z x y θ Fig. 1. Measurement setup. The xyz coordinate system was defined by the light propagation direction z and its associated perpendicular plane xy (y towards the reader). The solid rubber phantom traveled in the x direction on a motorized linear stage. Inset: configuration for measuring signals from flowing intralipid. The beam was reflected almost perpendicular to the fiber after internal reflection (a rotation around y). The angle between the beam and the flow was θ = 82.5 • .
has the following form when M m (M s ) is the fraction of moving (static) scatterers, and all moving scatterers inside the resolution volume follow the same behavior Here the first exponential is related to the Doppler signal, the second exponential is related to the decorrelation of the amplitude signal due to the diffusion of the scatterers, and the last two exponentials are related to the decorrelation of the amplitude of the signal due to movement of the scatterers in the axial z and transverse xy directions. n is the refraction index, k 0 the central wave number of the spectrum, D the diffusion constant of the scatterers, w z the axial 1/e size of the PSF of the system, and w t the transversal 1/e size of the PSF of the system. M m can be considered as a figure of motion contrast, because areas where the OCT signal has a significant value of M m correspond to areas with motion, such as vascular networks. We do not include the term for entering/exiting scatterers into the resolution volume in [21], because the PSF of the system is Gaussian and by definition it does not have hard edges that particles can enter or exit. In Eq.
(1) and in the following equations, the integral is replaced by a summation in the discrete sampling case. The second-order correlation function deals with the intensity I(t) = |E(t)| 2 and is defined as This function is related to the first-order autocorrelation through the Siegert relationship when the probability distribution of E is Gaussian, as in the case of speckle. If we assume M m = 1 for simplicity, Eq. (3) becomes It is advantageous to consider ensemble averages when estimating g (2) (τ) experimentally. The normalized second-order correlation function makes use of the concept of observation ensemble to make a more robust estimation. It is defined as [31] g (2) (τ) = I(t)I(t + τ)dt where I is the intensity, and · · · denotes ensemble average. The ensemble average allows to take into account multiple correlation windows in the calculation, reducing the statistical fluctuations on speckle size and its effect on the calculated decorrelation time. Although definitions Eqs.
(1) and (3) are commonly used in dynamic light scattering [32,21], in the ultrasound and OCT speckle decorrelation literature the Pearson correlation function g (P) has been traditionally used instead, either between sequential A-lines [17] or at a single depth as a function of time [33,34,19]. In the case of the autocorrelation of the intensity, it is defined as whereĪ denotes the average intensity and σ I the standard deviation. g (P) is defined between -1 (for anticorrelated) and 1 (for correlated signals), having a value of 0 for decorrelated signals. For this reason, comparison between the Pearson autocorrelation and g (2) requires subtraction of 1 from g (2) , because decorrelated signals produce g (2) = 1. The (I(t) −Ī)/σ I terms are equalizing variations in the intensity, regardless of their magnitude, increasing the effect of noise. For example, static or low speed areas with high, near constant intensity will exhibit strong decorrelation because the fluctuations of the noise determine the value of σ I . In opposition to Eq. (6), the definition in Eq. (5) carries information about the speckle contrast. In particular, g (2) (0) = 2 for fully developed, unity contrast speckle. For this reason, we define a motion contrast function C as C can be used to discriminate between static and dynamic areas of the tomogram, and is a good indicator of the presence of flow, similar to the conventional speckle variance technique [14]. Experimentally, static areas with strong signal and low noise will produce a small value in C. It is possible to use the definition in Eq.(5), replacing I(t) by A(t) = |E(t)|. This modified autocorrelationg (2) (τ) behaves asg (2) in our particular case, where γ takes into account the different intercept of this function.

Background noise effects in the time autocorrelation of OCT signals
Now consider the autocorrelation of a signal that is affected by additive noise. We can express the total detected signal, up to a multiplicative constant, as where S is the original signalS normalized by its standard deviation σ S , and similarly for the noise N. The relative weight of each if given by 0 ≤ w ≤ 1, and serves as a measure of the signal-to-noise ratio (SNR), which in dB equals We define the cross-correlation of functions f and g as R f g . The autocorrelation of I using Eq. (3) is In the case of OCT, we will consider N(t) as background noise. Effects of other types of noise (such as that produced by multiple scattering) will not be considered here. It is important to note that, rigorously, noise is added to the complex OCT signal however the terms that arise are essentially the same as those in Eq. (11), with duplicity to account for the real and imaginary parts of each contribution. Because this does not change the behavior that we will describe below, we will focus on Eq.(11) for simplicity.
In the present case, the noise will have no correlation with the signal, which implies R SN = R NS = 1 as mentioned above. If we additionally assume that the frequency spectrum of the noise n( f ) is white, we have where δ denotes the impulse function. For this reason, background white noise will affect the very first point of the autocorrelation τ = 0, but τ > 0 will contain w 2 R SS . Because w 2 is a constant, it will not modify the shape of the autocorrelation, which enables, in principle, to omit the τ = 0 point in the autocorrelation function to regain R SS . This is what we consider the primary effect of background noise. If the Pearson autocorrelation is used instead, we have The denominator arises because of the change in the standard deviation of the combined signal, and the third term vanishes because R (P) Experimentally R II is being estimated by a limited, discrete sampling in time I(t) → I t with total number of samples n s . For this reason, although N(t) → N t will be generated by white noise, its discrete autocorrelation function, being estimated from a finite number of discrete time points, will not be a perfect impulse function. We can expect that R NN will be given by with | f τ (τ)| < 1. f τ (τ) will behave as white noise itself and oscillate around zero, but the weight of its mean square value will be directly related to the number of samples n s used to estimate the autocorrelation function. Numerical simulations suggest that the variance of f τ behaves as σ 2 f τ = 1 n s . When using the Pearson definition, R (P) NN = R NN − 1. f τ will directly affect the estimation of the autocorrelation R II with weight (1 − w) 2 , and will start dominating the autocorrelation function for small values of w and n s . The autocorrelation function becomes with a similar definition when using the Pearson autocorrelation. The primary effect of background noise is seen as an additional drop ∆g in R II between the first τ = 0 and second time points τ = ∆t with height If the SNR of the signal is known, ∆g can be used to compensate for the primary effect of noise without losing the first point of the autocorrelation function.. Due to the random nature of f τ , the autocorrelation function R II when noise dominates will consist of oscillations with near-unity magnitude in terms of the normalized autocorrelation function. When calculating the decorrelation time defined as the crossing of the autocorrelation function below a given threshold, these oscillations will make very probable a crossing in the first few time difference samples. For this reason, once noise starts to dominate the autocorrelation function, we can expect for the correlation time to drop, reaching a point where the correlation time will be very short irrespective of the dynamics of the scatterers in the sample. This is what we refer to as the secondary effect of background noise.
As mentioned before, to compensate the primary effect of background noise we can define a new discrete autocorrelation function as where ∆t is the time interval between A-lines, and the new motion contrast value that replaces Eq. (7) is With these definitions, the main contribution from R NN is suppressed by omitting the value at τ = 0, and the subtraction by median{g (2) (τ > n s /2)} will allow a renormalization of the autocorrelation function between 0 and 1, by taking as 0 the median of f τ when the pure signal autocorrelation R SS is negligible. This assumes that n s /2 is larger than any time scale of the dynamics of the sample. Although this enables the calculation of the decorrelation time for medium to high SNRs, there is still a significant range of low SNRs that will be affected by the secondary noise effect. The determination of the behavior of the secondary effect, and how to account for it, will be considered in Sec. 4.2.

Relation between speckle decorrelation time and speed
Assuming an isotropic voxel (w z = w t ) for simplicity, and representing the total speed as v = v 2 z + v 2 t , we can find using Eq.(4) the correlation time τ c when the autocorrelation function g (2) (τ) [g (2) (τ)] reaches the threshold 1 + g c [1 + √ g c ] for moving and diffusing scatterers as with the diffusion constant K D = 2n 2 k 2 o D andĝ c = −1/ ln(g c ), which is positive definite because g c < 1. The different definition of the threshold when usingg (2) (τ) from Eq. (8) makes the rest of the derivation the same for both cases. For autocorrelations defined between −1 and 1, 1 is subtracted from the threshold. Now, if the diffusion term is small compared to the speed term so that v 2 /w 2 t K D , we arrive at the familiar result which is also in agreement with full-field techniques [35,36]. Note that for most common thresholds 1 ĝ c 3.
In a practical implementation of speckle decorrelation measurements in OCT, it is important to know that OCT exhibits two speckle sizes that come from different origins. The speckle size in the axial direction is related to the axial resolution of the system. When A-lines are acquired as a function of time, the intensity at each depth will evolve in time. We can attribute a speckle size along a second axis to this temporal evolution, which will exhibit a decorrelation time given by Eq. (20), which is linked to the dynamics of the scatterers reflecting light at that depth. Figure 2 shows a clear example of this relation, where we show an OCT tomogram of a tube filled with flowing liquid (inset of Fig. 1) at two different flow rates. The tomograms are so-called M-mode measurements, which consist of the measurement of A-lines as a function of time while maintaining the OCT probe position constant. For this reason the horizontal axis in Fig. 2 corresponds to time and the vertical axis to depth. The beam enters at the top, where the reflections from the collimating lens and other surfaces can be seen. The tube material that encloses the scattering liquid can be seen in the area of weak scattering at the bottom of the tomograms. It is easily seen that the speckle-decorrelation time (and therefore the speckle size in the horizontal direction) matches the expected relationship with scatterer speed. Note how the horizontal speckle size gets larger near the interfaces due to the no-slip condition of the flow.
For simplicity, in the rest of this work we will consider exclusively motion in the x direction, although diffusion still exists in three dimensions. In the case of movement of a liquid, Eq. (20) can be considered as the inverse decorrelation time as a function T of all the contributions to the decorrelation of the signal where k absorbs the parameters dependent on the beam characteristics and the chosen threshold, and might be dependent on z if the beam is not well collimated in the area of interest. K BM contains the contribution from diffusion and the chosen threshold as well. It is possible have more complex scenarios in which different processes apart from scatterer translational motion affect τ and can be included inside T , especially in the case of anisotropic scatterers. For example, it is known that red blood cells (RBC) can exhibit motions such as tumbling, tank-treading, swinging and deformation under certain circumstances [37], which would add additional terms to T , and which also potentially depend on other factors such as temperature and RBC concentration.

Experimental results
The setup of Fig. 1 was arranged using a solid rubber phantom placed on a motorized linear stage for precise speed control. The OCT probe was placed near the surface of the sample. Because the optics of the OCT probe were designed to be used while immersed in water, we expect to have a small defocus which modifies w t in Eq. (20), and therefore k in Eq. (22). The sample was moved in front of the beam at different speeds ranging from 40 to 600 mm/s, in steps of 40 mm/s. Because the sample was a solid we have K BM = 0. Except otherwise noted, the correlation window used was n s = 64, and the ensemble consisted of 32 correlation windows in time (2048 A-lines total) and 8 px in depth. This means that each group of 2048 samples in time and 8 px high (roughly 32 µm) produce one point in the inverse speckle decorrelation-time profile, at a rate of approximately 50 depth decorrelation profiles per second. We found that using the OCT signal amplitude |E| as described in Eq. (8) produced decorrelation profiles that were significantly more homogeneous than using the tomogram intensity. This seems to be linked to the effect that the square has on outliers of the amplitude in the statistical fluctuations of speckle intensity. We speculate that these outliers might be produced by extraneous reflective particles in the liquid. Because of γ, this modifies the range of values for C, and all autocorrelations had a maximum C around 1.15, the reason for which the threshold value chosen was √ g c = 0.025 for the traditional autocorrelation function, and √ g c = 0.5 for the renormalized version. A linear interpolation was used to calculate the τ c time-difference crossing of the threshold at a sub-∆t accuracy.  Fig. 3 we show only the top 80 px in order to simplify visualization. The OCT probe presented an artifact due to multiple reflections in its interfaces, which is clearly visible in the decorrelation time profile as two lines of high τ −1 near the top of the sample. Because these reflection artifacts varied considerably depending on the particular probe, we avoided these areas in the subsequent analysis.

Primary effect of background noise
The sample had an irregular surface and discontinuities that appear as shallow lines with low τ −1 . The deep discontinuities correspond to the end of the travel of the translation stage, a point where the sample came to a stop before gaining speed in the opposite direction. The acceleration was chosen in order to go from a complete stop to nominal speed in less than 30 ms (approximately 3 time pixels in the τ −1 profile for a complete deceleration-acceleration cycle). Ignoring these discontinuities, it is easy to see that when using Eq. (5) there is a strong gradient with depth, even though the sample was a solid, moving at a well-defined speed. This gradient is more pronounced as the threshold √ g c increased. Note, however, that the current threshold is too low for crossing the autocorrelation function for points near the top of the sample, where the SNR is very high: although the sample starts at a depth of 20 px, measurable decorrelation is only found deeper than 24 px. When using our new autocorrelation function Eq. (18), the profile is significantly more homogeneous, which confirms our derivation on the primary effect of noise on the autocorrelation function. Also, the area with high SNR shows a (a) In order to appreciate more clearly the effect of noise in τ −1 , the top row in Fig. 4 shows the traditional autocorrelation function at different depths and at different sample speeds, and the bottom row shows the new autocorrelation function of Eq. (18). The depth of 20 px corresponds to a point just after light enters the sample, where there is a significantly higher SNR level than in the rest of the depths shown. The dashed line represents the threshold 1 + √ g c that was used to calculate the correlation time profile of Fig. 3. The autocorrelations show clearly the excessive drop in correlation between τ = 0 and gτ = ∆t produced by the primary effect of background noise, described by Eq. (17). This dip increases with depth as the signal decays as expected. In contrast, the new autocorrelation functions at the bottom row exhibit higher similarity at a given speed, including the autocorrelation at pixel depth 20. The secondary effect of background noise is also visible, as an oscillation for low g (2) and large depth, as well more spread in the correlation time at low speeds. This indicates that addressing the first effect of noise enables the use of speckle decorrelation in a limited range of SNRs. In the next section we will expand the SNR range by addressing the secondary effect of noise.

Secondary effect of background noise
To understand the contribution from f τ and make a calibration of the proportionality constants in Eq.(22), we first had to define a way to determine the SNR at each correlation window. We estimated the noise floor given by the presence of background noise by measuring a tomogram when the sample arm of the interferometer was blocked. This depth-dependent background sig- nal was then used as the noise level, and the average intensity at each correlation window was used as the signal intensity to calculate the SNR. The scatter plot in Fig. 5(a) shows the valid data points from the τ −1 profiles as a function of SNR, color coded with the sample speed. To asses the performance of the Pearson autocorrelation function in the presence of noise, Fig. 5(b) shows τ −1 as calculated using g (P) in Eq. (6), and modified following Eq. (18). The same correlation window was used, and ensemble averaging was performed by simple averaging among the g (P) calculated for each member of the ensemble. As described previously, the amplification of the noise fluctuations is evident at low speeds, which makes speeds indistinguishable in the dashed box area. Paradoxically, τ −1 estimations using g (P) have smaller dispersion than when using g (2) for a given speed and sufficient SNR, as seen by the more diffuse nature of points in Fig. 5(a) versus Fig. 5(b) at a given speed. Figure 5(a) shows that there is a small increase in τ −1 as the SNR decreases from 40 dB to 16 dB. This was mostly uniform for all sample speeds. In this region the relation between τ −1 and v is linear, with a non-linearity at higher speeds. This is expected, as speeds that decorrelate faster than the time resolution of the system should all give the same τ −1 . Below 16 dB SNR the behavior changes depending on the sample speed. High speeds preserve the slow increase in τ −1 , but the dispersion in τ −1 values increases. Low speeds have a sharp increase in τ −1 until the SNR reaches approximately 10 dB. At this point all the speeds are overlapped, and it is not possible to uniquely identify the speed from the inverse correlation time and the SNR. This happens because at this point f τ completely dominates the autocorrelation function. The limit on the maximum measurable τ −1 is given by the time interval between A-lines ∆t and the elimination of the zero delay autocorrelation value, and is given by τ −1 max ≈ 1/(2∆t) ∼ 50 ms −1 . This can be increased if interpolation is used.
The data in Fig. 5(a) was used to make a parametrization F of the behavior of τ −1 as a function of SNR and calibration speed v x,c for the areas in which F(τ −1 , SNR) is single valued. Selected SNRs usingg (2)   The function F is defined as where we consider the effect of f τ in Eq. (16) as an additive term F SNR (SNR) when calculating the autocorrelation time. The SNR response should be independent of the PSF size, and for this reason we introduce the calibration k constant from Eq. (22) as k c = √ 2ĝ c /w c , where w c is the lateral PSF during calibration. This factor implies that F(τ −1 , SNR) is in inverse time units, an can be considered a noise-corrected inverse correlation time τ −1 corr . In our case we did not know the exact value of w c during calibration, so we used the value for the nominal lateral PSF and replaced the corrected τ −1 by an effective inverse correlation time τ −1 eff . The relation between the noise-corrected correlation time and the effective correlation time is In the case of F τ (τ −1 ), it is clear that the relation between v and τ −1 for a solid should be linear according to Eq. (21), except at low speeds where the correlation window is too small (speeds with correlation times τ > n s ∆t), and at very high speeds (with respect to ∆t, with correlation times τ τ max ). All sample speeds in Fig. 5(a) fulfilled τ < n s ∆t as we are not interested in very low speeds, so the non-linearity at low speeds will not be considered here. For this reason, we proposed F τ (τ −1 ) as given by polynomials of different orders. We found that the simplest polynomial that described well the experimental was given by where the linear term dominates, except for correlation times near the limit of the system τ −1 max where the 5-th order term starts to make significant contribution. In the case of F SNR (SNR), we assumed there was a contribution to τ −1 from noise that increased asymptotically at low SNRs according to The position of the asymptote had a small dependence on τ −1 , as seen in Fig. 5(a).
The rubber material used for this calibration presented significant multiple scattering, which is a source of noise in addition to background noise. It is feasible to expect a change in behavior at low SNRs for materials in which single scattering is predominant. For this reason, similar measurements were carried out using an intralipid solution in water concentrations ranging from 0.5% to 2.0% in volume. In this concentration range, the samples exhibited highly isotropic, single scattering. The solution was placed in a fully filled cylindrical container, and placed on the translation stage. This configuration avoided convection effects due to the acceleration. Very similar results where obtained for all concentrations, and the mean of all parameters were k 0 = 2.4 ms −1 k 1 = 1.28 k 5 = 3.6 · 10 −8 ms −4 α 1 = 5.56 ms −1 α 2 = 15.5 dB α 3 = −3.51 · 10 −2 dB · ms α 4 = 5.55 · 10 −1 , (29) It is clear that behavior at low SNRs changes, but the general behavior given by k 1 is maintained.
We attribute the small difference in k 1 to a difference in the PSF due to geometrical differences: in the solid phantom, the light entered a flat surface which introduces mainly spherical aberration [38], while in the liquid phantom, light entered a curved surface which acts mainly as a lens. Figure 6 shows the solid sample speed as determined by using the model. Figure 6 contains the data for a different set of measurements from those used in the calibration. It is clear that the determined speed is highly homogeneous, and that it is difficult to reliably determine the speed at large depths due to the decreased SNR.

Measurements in presence of diffusion
In order to test the technique in a sample with diffusion, we performed measurements using an intralipid solution in water. To account for the different divergence of the beam in this case, we will take into account the effect of the small correction factor introduced by Eq. (24) in Eq. (22). Considering that using the nominal PSF now we have using Eq. (25) to use our effective correlation time parametrized in the calibration we obtain where for simplicity we definedK BM = k n /k c K BM andk = k n /k c k n , noting that k n /k c = w c /w n . The liquid measurements used an ensemble size in the axial direction of 4 px. The measurements were carried out as shown in the inset of Fig. 1, with an angle between the optical beam and the flow direction. This angle is close enough to the perpendicular in order to expect the same behavior of speckle decorrelation time as a function of scatterer speed as in the SNR calibration, but at the same time allowing for a small Doppler signal to be acquired from the liquid. The power of using F(τ −1 , SNR) in the equations is evident in the case of measuring a liquid: all the complexity regarding the time-resolution of the system, the autocorrelation function used and the effect of noise is encapsulated in it. More complex systems can incorporate additional effects, such as those mentioned before in the case of red blood cells. The only calibration necessary at this step is to find K BM (andk in our particular case, because we do not know the value of k c ). Afterwards, the flow speed can be found easily by inversion of Eq. (31) In order to find the constants of Eq. (32), we performed measurements at different flow rates using the setup shown in the inset of Fig. 1. Each dataset was analyzed using both Doppler and speckle decorrelation. At this step, we included in the analysis the depths that were corrupted by the secondary reflection of the probe, manifested near pixel 30 in Fig. 3. Because this reflection artifact does not change with time, it is possible to account for it by assuming a different value ofk in this area. The net effect is an increase of the value ofk, which makes this area more sensitive to motion. The calibration process to find the constants of Eq. (32) consisted of using the Doppler signal as the absolute velocity reference in a point-per-point fashion, for different flow rates. The Doppler velocity was corrected for the line-of-sight (LOS) angle, θ , as v = v Doppler / cos θ , and phase unwrapping was required at high flow rates due to the reduced dynamic range of Doppler. In order to minimize the statistical variation of the speckle value during the calibration, each data set was time-averaged before plotting the relationship between τ −1 eff and v. We determined the parameter values by performing a least-squares fit to Eq. (32) as a function of depth to allow for changes in the PSF with depth as well as to account for the effect of the reflection artifact. The data and the model fit are shown in Figs. 7(a)-7(b) with the following results (averaged in depth in the area outside the artifact) It is clear that there is an excellent agreement between the speckle decorrelation and the Doppler speed measurements, and the speckle measurements are consistent as a function of time and for different flow rates. Even the area with the reflection artifact shows good agreement after calibration, with only a slight trace of the artifact. It is important to point out that the Doppler signal had to be unwrapped, because of its reduced dynamic range. This is possible in a smooth flow profile such as the one shown here, but unwrapping is not always possible depending on the characteristics of the profile being measured. Figure 8 presents the results when the calibration obtained above was applied to measurements obtained with 20, 40, and 80 mL/min flow rates. The LOS-corrected speed profiles as determined by Doppler processing and after unwrapping are also presented for reference. Note that the scales for the speckle decorrelation and Doppler flow profiles for each flow rate are the same. It is clear that at these flow rates the Doppler signal needs to be unwrapped as there are speeds well beyond the Doppler resolution limit of the system (∼ ±25 mm/s, ∼ ±220 mm/s when corrected for our particular LOS). The results in Fig. 8 demonstrate that the theoretical description in this work describes very precisely the behavior of speckle-decorrelation time as a function of scatterer speed, and confirms our approach for the calibration of the secondary effect of noise and of the PSF as a function of depth in water. Figure 8 illuminates the advantages and disadvantages of speckle decorrelation motion measurements. It clear that speckle measurements are inherently more affected by noise than Doppler measurements due to the statistical nature of the measurement. More importantly, the   dynamic range of speckle decorrelation is considerably larger than for Doppler: the maximum measurable speed is > 600 mm/s, or more than 24-fold greater than possible with LOS Doppler. On the other hand, measuring low speeds is challenging: the slower parts of the flow in Fig. 8 show a deviation between speckle and Doppler, either in the form of oscillations near the OCT probe, or as an increased speed near the opposite wall.
The difficulty in measuring low speeds in speckle decorrelation is an intrinsic property of the technique in the presence of Brownian motion. Consider a measurement of a liquid with no flow, where ideally τ −1 eff = 2K BM . There will be fluctuations in τ −1 eff that come from the estimation of the autocorrelation function from a finite measurement. These fluctuations can be seen as an additional term ξ ( j) in the j-th correlation window, so that The particular distribution of ξ is difficult to determine, because it has contributions from different sources with different statistics. The discrete nature of the autocorrelation function indicates that a Poisson distribution best describes the determination of the correlation time. We will consider, for simplicity, that we can approximate the Poisson distribution with a Gaussian distribution, and that |ξ ( j)| < 2K BM . Also note that the average over j, ξ ( j) j , is zero due to the absence of flow. Using Eq. (32), we obtain a distribution of velocities where we indicate that the real part of the square should be taken, as v x ( j) can only be real. We can see that v x j ∼ |ξ | j /k > 0 (except for ξ ( j) = 0, ∀ j), which explains why even at zero speed, on average speckle decorrelation will always measure some finite speed in a liquid. This is not necessarily the case for dynamic light scattering, because the curve fitting done in this case is able to distinguish, to some degree, between the exponential decay of diffusion and the squared exponential decay of speed. However, even in the case of phase-stable acquisition, this kind of non-linear fitting is very complex and requires ad-hoc fitting algorithms [21] for g (1) . In an intensity-based dynamic light scattering approach, the fitting must be done for g (2) , which is even more complex due to the crossed terms. This makes dynamic light scattering difficult to adapt for real-time measurements.

Conclusions
In spite of the great interest in reliable speed measurements based on decorrelation analysis of the intensity information of OCT data, previous studies have not demonstrated calibration and validation. In this work, we have shown the limitations of using a direct implementation of the traditional autocorrelation function for speed measurements in the presence of noise, and have presented an alternative autocorrelation function that is considerably more robust. Although we focused on speckle decorrelation, the same approach can be used in dynamic light scattering OCT to account for background noise contributions to decorrelation.
We derived a direct expression for the relation between speckle decorrelation time and speed and diffusion and demonstrated a calibration procedure that allowed us to perform quantitative measurements of speed using speckle intensity statistics in a wide range of SNRs. We analyzed the effects of Brownian motion when measuring speeds in liquids, and found an unavoidable non-linearity at very low flow speeds. Our technique demonstrated a high dynamic range that easily surpasses that of Doppler techniques. We validated our results in two different scenarios, using a calibrated linear stage and corrected Doppler measurements. Our experiments delimit more clearly the range of validity of the speckle decorrelation technique. Moreover, our results demonstrate that phase-unstable OCT systems, which are far easier, faster and cheaper to build and calibrate, can now be used to measure speeds quantitatively, repetitively and with a high dynamic range, thereby enabling a host of biomedical applications.