High speed, wide velocity dynamic range Doppler optical coherence tomography (Part I): System design, signal processing, and performance

Improvements in real-time Doppler optical coherence tomography (DOCT), acquiring up to 32 frames per second at 250 × 512 pixels per image, are reported using signal processing techniques commonly employed in Doppler ultrasound imaging. The ability to measure a wide range of flow velocities, ranging from less than 20 μm/s to more than 10 cm/s, is demonstrated using an 1.3 μm DOCT system with flow phantoms in steady and pulsatile flow conditions. Based on full implementation of a coherent demodulator, four different modes of flow visualization are demonstrated: color Doppler, velocity variance, Doppler spectrum, and power Doppler. The performance of the former two, which are computationally suitable for real-time imaging, are analyzed in detail under various signal-to-noise and frame-rate conditions. The results serve as a guideline for choosing appropriate imaging parameters for detecting in vivo blood flow. 2003 Optical Society of America OCIS codes: (110.4500) Optical Coherence Tomography; (170.3880) Medical and biological imaging; (170.1650) Coherence imaging; (100.6950) Tomographic image processing References 1. D.Huang, E.A.Swanson, C.P.Lin, J.S.Schuman,W.G.Stinson, W.Chang, M.R.Hee, T.Flotte, K. Gregory, C.A.Puliafito, and J.G.Fujimoto, “Optical coherence tomography,” Science 254, 1178-81 (1991). 2. Z.Ding, Y.Zhao, H.Ren, J.S. Nelson, Z.Chen, “Real-time phase-resolved optical coherence tomography and optical Doppler tomography,” Opt. Express 10, 236-45 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-5-236. 3. V.Westphal, S.Yazdanfar, A.M.Rollins, and J.A.Izatt, “Real-time, high velocity-resolution color Doppler optical coherence tomography,” Opt. Lett. 27(1), 34-6 (2002). 4. A.M. Rollins, M.D. Kulkarni, S. Yazdanfar, R. Ung-arunyawee, and J.A. Izatt, “In vivo video rate optical coherence tomography,” Opt. Express 3, 220-29 (1998), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-3-6-219. 5. M. Laubscher, M. Ducros, B. Karamata, T. Lasser and R. Salathe, “Video-rate three-dimensional optical coherence tomography,” Opt. Express 10, 429-35 (2002), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-10-9-429. 6. S. Yazdanfar, M.D. Kulkarni, and J.A. Izatt, “High resolution imaging of in vivo cardiac dynamics using color Doppler optical coherence tomography,” Opt. Express 1, 424-31 (1997), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-1-13-424. 7. J.A.Jensen, Estimation of blood velocities using ultrasound (Cambridge, 1996). 8. G.J. Tearney, B.E. Bouma, and J.G. Fujimoto, “High-speed phaseand group-delay scanning with a grating-based phase control delay line,” Opt. Lett. 22, 1811-3 (1997). 9. V.X.D.Yang, M.L.Gordon, A.Mok, Y.Zhao, Z.Chen, R.S.C.Cobbold, B.C.Wilson, and I.A.Vitkin, “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation,” Opt. Commun. 208, 209-214 (2002). 10. Y.Zhao, Z.Chen, C.Saxer, Q.Shen, S.Xiang, J.F.de Boer, J.S. Nelson, “Doppler standard deviation imaging for clinical monitoring of in vivo human skin blood flow,” Opt. Lett. 25, 1358-60 (2000). 11. J.F. de Boer, C.E. Saxer, J.S. Nelson, “Stable carrier generation and phase-resolved digital data processing in optical coherence tomography,” Applied Opt. 40, 5787-90 (2001). 12. S.I. Fox, Human Physiology (W.C. Brown, Dubuque, Iowa, 1987). 13. B.E. Bouma, G.J. Tearney, Handbook of Optical Coherence Tomography (Marcel Dekker, 2002). (C) 2003 OSA 7 April 2003 / Vol. 11, No. 7 / OPTICS EXPRESS 794 #2105 $15.00 US Received February 04, 2003; Revised March 31, 2003


Introduction
Optical coherence tomography (OCT) [1] and its Doppler variants, optical Doppler tomography (ODT) [2] or color Doppler optical coherence tomography (CDOCT) [3], are analogous to ultrasound imaging where light is used instead of sound.High speed microstructural OCT imaging has been demonstrated at various pixel-resolution and framerate combinations, using a rapid scanning optical delay line (RSOD) [4] to achieve 250 × 125 pixels at 32 frames per second (fps) or a smart-pixel detector array [5] to acquire 58 × 58 × 58 pixels at 25 fps.To image microvasculature blood flow, the velocity sensitivity can be increased by decreasing the frame rate.Most demonstrated systems [2,3] use averaging and/or correlation algorithms to obtain the velocity information, thus limiting the effective number of lines per second available for color Doppler imaging to a fraction of that in Bmode, depending on the number of axial scans used for averaging and/or correlation (typically 4 ~ 16).The high spatial resolution (1 ~ 20 µm) of OCT and high velocity sensitivity (~ 50 µm/s) of ODT allow in vivo imaging of microvasculature, which can be useful for studying angiogenesis, embryo cardiac dynamics in developmental biology [3,6], and vascular response to treatments that are beyond the reach of current clinical ultrasound systems.
There are considerable differences between ultrasound and OCT (and their Doppler extensions), such as the underlying physical interactions, the different wavelength involved, and the much higher velocity of light than sound.Nevertheless, from the signal-engineering point of view, there are many similarities.Since ultrasound signal processing is a mature field with many validated techniques for real-time imaging [7], it can guide the design of highspeed OCT/ODT hardware and software processors.Here, we present a comprehensive approach to the design and implementation of a hardware processor (an analog coherent demodulator) and software processors that perform signal amplitude and phase estimation and motion artifact rejection.This system allows blood flow visualization in a number of modes, including color Doppler, velocity variance, Doppler spectrum, and power Doppler.We will refer to the general technique as Doppler optical coherence tomography (DOCT), and describe the signal-processing algorithms for each mode in detail.The performance of the DOCT system will be characterized by a number of stationary and flow phantoms.This information will allow appropriate trade-off decisions for detecting blood flow in different in vivo imaging applications.

DOCT system
A schematic of the DOCT system is shown in Fig. 1(a).A broadband light source (JDS, Ottawa, ON, Canada) with a polarized output power of 5 mW at center wavelength λ 0 = 1.3 µm and bandwidth ∆λ = 63 nm is coupled into a single-mode fiber optic interferometer.A rapid scanning optical delay line and phase modulator (JDS, Ottawa, ON, Canada) form the reference arm.The light is focused to a spot size of ~ 25 µm in the sample arm, which can be scanned in three dimensions via high-speed linear translators (Newport, Irvine, CA) and/or a galvanometer (Cambridge Technology, Cambridge, MA).Since the RSOD and the phase modulator are polarization sensitive, multiple polarization controllers (OZ Optics, Ottawa, ON, Canada) are used to optimize the signal transmission in the reference arm and to match polarization states between the reference and sample arms.An optical circulator (O-Net, Boston, MA) provides balanced optical detection to reject the high DC background and increase the signal-to-noise ratio (SNR).A wide bandwidth (~ 10 MHz), custom-made, transimpedance amplifier detects the AC interference signal.This signal, high-pass filtered at 30 kHz, is passed through a depth-gain-compensation amplifier, analogous to the time-gaincompensation amplifier employed in ultrasound, to increase the gain as a function of imaging depth.A coherent demodulator yields both the amplitude and phase of the AC signal, as explained later.Care must be taken to minimize the frequency response differences between the in-phase and quadrature channels of the analog hardware.Residual differences in DC bias and AC gain can be balanced digitally after the signals are analog-to-digital (A/D) converted.These steps are critical to reduce background noise in both structural and flow imaging and to ensure accurate estimation of flow velocity.
One of the most important challenges in real-time DOCT imaging with high spatial/velocity resolution is achieving high axial scan frequency, analogous to pulse repetition frequency (PRF) in ultrasound.Commercial clinical ultrasound imaging systems (transducer center frequency 3 ~ 8 MHz) typically have PRFs in the range of 4 ~ 10 kHz for array transducers, whereas high frequency (20 ~ 100 MHz) systems may exceed 20 kHz for single element transducers [7].For a given image size (number of pixels), the B-mode imaging frame rate is proportional to the PRF.In the case of flow imaging, signal averaging between axial scan lines is often implemented, and higher PRFs with lower frame rates typically result in better flow images.
Fourier-domain rapid scanning optical delay lines have been used in a number of highspeed OCT systems, offering axial scan frequencies from ~ 1 kHz (galvanometer [8]) to 4 kHz (resonant scanner [4]), corresponding to frame rates of 2 ~ 8 fps for 500 line images, if only the forward scans are used to form the images.The frame rate can be doubled if both the forward and reverse scans are recorded.To achieve real-time (e.g., 32 fps) imaging with high resolution (e.g., 500 lines per image), the axial scan frequency must be increased to at least 8 kHz (Table 1).Only resonant scanners can achieve such high frequencies.The resultant OCT signal bandwidth varies sinusoidally as [4]: where D is the axial scan depth, f a is the axial scan frequency, ∆λ = 63 nm is the light source (with a Gaussian spectrum) bandwidth, and λ 0 = 1.3 µm is the source center wavelength.We have constructed two RSODs operating at 8.05 kHz (D = 1.3 mm) and 12.95 kHz (D = 1.0 mm), resulting in maximum signal bandwidths of 2.4 and 3.0 MHz, respectively, assuming a Gaussian source spectrum (Fig. 2).The RSODs are adjusted to provide group delay only by centering the light spectrum on the scanning mirror.A phase modulator operating at 4.3 MHz provides the phase delay using a saw-tooth waveform triggered by the RSOD (Fig. 1).This arrangement ensures a stable carrier frequency (f C = 4.3 MHz) [11], which is essential for hardware coherent demodulation, and allows the use of small mirrors with low inertia in the RSOD to achieve high resonant frequencies.Compared to sinusoidal waveform, the saw-tooth phase modulation also suppresses the sidebands, which contribute to noise in the OCT signal.The full implementation of the coherent demodulator involves mixing the detected OCT signal with not only , yielding the in-phase (I) and quadrature (Q) components, as shown in Fig. 1(b).Since the I and Q signals are shifted to the base-band, low-pass-filters (LPF) instead of band-pass-filters (BPF) are used.The LPF used in our system has a bandwidth of approximately 1.55 MHz, and greater than 70 dB stop-band rejection.Such a LPF is easier to construct than a BPF with equivalent performance, and is required to avoid aliasing before A/D conversion.The required base-band filtering bandwidth is 1.2 and 1.5 MHz, for the 8.05 and 12.95 kHz RSOD, respectively (Fig. 2).(2) To increase data throughput from the A/D card to computer memory, an A/D conversion clock pulse train triggered by the RSOD is used, with data acquisition occurring during the central ~ 80% of the depth scan.The acquired signals can be expressed as 2D arrays: I m,n and Q m,n with m and n denoting the indices in the depth and lateral directions, respectively.These signals are transferred by direct memory access without computer CPU intervention, reducing the computational overhead.

Digital signal processing
In structural imaging mode, the intensity of the OCT signal is proportional to the local reflectivity at a particular pixel.Spatial averaging in the depth and lateral directions are performed to improve the SNR.The value of a given pixel in the structural image is calculated as: Due to the large dynamic range involved, the logarithm of the OCT signal is displayed in the structural image.Simultaneous to the structural OCT imaging, blood flow information can be presented in one of several modes, as the following.

Color-Doppler mode
The mean velocity at any pixel can be evaluated by the Kasai autocorrelation equation [7,9]: ( ) where f D is the Doppler frequency shift, n t is the tissue index of refraction (~ 1.4), θ is the Doppler angle, and Y, X will be used for estimating the variance of velocity estimate.The arctangent function evaluates the mean phase shift between axial scans (Fig. 3) and is calculated in all four quadrants with outputs ranging from -π to +π.The calculation of mean phase change from Eq. ( 5) is numerically efficient, especially in fixed-point, because it relies mostly on summations and multiplications, with only one each of division and arctangent operations using floating-point representation to maintain accuracy (Table 2).

# of calculations (single pixel) # of calculations* (on average per pixel)
* Estimated over the entire image with the total number of pixels >> MN The velocity estimation is performed within an N M × window (Fig. 4(a)), where M is the depth window length (often called the gate length in ultrasound literature) and N is the number of axial (depth) scans used in the autocorrelation (often called the ensemble length).Additional computation efficiency can be achieved if the N M × windows for adjacent pixels overlap when the flow image is processed, because the calculations for the numerator and denominator in Eq. ( 5) within the overlapping portions (αΜ and βN) can be shared between pixels.The accuracy of the velocity estimation improves via larger M and N, i.e., with more spatial averaging, as shown in Fig. 4(c).For a given axial scan frequency (f a ) and for a given number of axial-scans displayed per frame (K X ), the frame rate is: when only inward or outward axial scans are used.Clearly, a higher frame rate can be obtained by increasing β, at the expense of increased blurring in the lateral direction due to overlapping of sequential axial scans.In this work, the values of αΜ and βN were determined as follows: Because the A/D conversion rate is more than four times higher than the LPF bandwidth (i.e., over-sampled), Μ is set to 4 and α to 0.75, which does not result in significant axial blurring.The number of pixels per axial scan is set to K Y = 512.For the 8.05 kHz RSOD, β and N values are set such that the resultant horizontal resolution is 500 pixels for all frame rates except 32 fps (Table 3).At high frame rates (16 and 32 fps), high β and low N values are chosen to maintain the horizontal resolution.At lower frame rates (2 -8 fps), low β and high N values are used to reduce lateral blurring and improve velocity resolution. .Although aliasing can be phase unwrapped in principle, it is difficult in practice in the presence of noise.In addition, because velocity estimation accuracy deteriorates with increasing flow velocity (Fig. 10(b)), there is a physical upper limit for phase unwrapping techniques.The non-aliased velocity range is which is 3.9 mm/s for f a = 8.05 kHz and is 6.3 mm/s for 12.95 kHz.The velocity resolution is determined by the phase stability of the interferometer and electronics [3,11]: where φ ∆ is the phase error of the system measured using a stationary reflector as the target.At N = 8 and f a = 8.05 kHz, φ ∆ is ~ 0.06 radians (Fig. 4(c)), corresponding to a minimum detectable velocity of ± 17 µm/s.The velocity dynamic range (VDR), is 40.5 dB for the 8.05 kHz system.This can be improved by increasing M or N at the expense of spatial blurring or frame rate.For the 8.05 kHz system, the VDR can be as high as 63.9 dB with a minimum detectable velocity of ± 2 µm/s when N = 128 (Fig. 4(c)).
During in vivo imaging, the biological sample can exhibit bulk motion, i.e., the entire sample can move relative to the DOCT system and create motion artifacts.We have shown previously [9] that these artifacts can be removed by calculating the velocity histogram along an axial line.The histogram of the axial line with K z pixels is defined as: where ( ) is the number of pixels at velocity level  ( ) (11) Equation (11) ensures that the histogram has sufficient velocity resolution.The histogram h i reaches a maximum at the bulk tissue motion velocity: , where (12) Note that two conditions must be met: (a) the sample acceleration must small, i.e., the change in bulk motion velocity is smaller than the velocity resolution during one axial scan, and (b) the blood flow must occupy only a small portion of the axial scan depth.The second condition depends on the velocity distribution and the amount of bulk tissue in the depth scan.For example, for plug flow within bulk tissue that occupies the entire scan depth, the histogram method breaks down when blood flow occupies more than 50% of the axial scan.For flow conditions with a broader velocity distribution (e.g., laminar), this constraint is more relaxed.
The corrected flow velocity, accounting for aliasing due to the block tissue motion, is: The velocities estimated by Eqs. ( 4) and ( 5) and corrected for bulk tissue motion with Eqs. ( 10) -( 13) are color-coded (see Fig. 4(b)) and form the color Doppler display mode of the DOCT system.
Even with high axial scan frequency, the non-aliased velocity range (Eq.( 7)) is insufficient to measure high-speed blood flow in arterioles and larger vessels [12].In these cases, it is beneficial to use a display mode with a wider velocity measurement range.One approach is to compute the velocity standard deviation or variance [10,13], which can vary monotonically with increasing flow velocity (see below).If * S is the complex conjugate of S , the normalized velocity variance is [7]: Equation ( 14) can be evaluated efficiently, since X , Y and 2 S have all been calculated for the color Doppler and structural display modes.Hence, obtaining the velocity variance display only slightly increases the overall computation.Velocity variance mapping eliminates aliasing and can greatly increase the velocity measurement range, but the blood flow orientation is lost.

Power Doppler mode
Another approach to velocity mapping is power Doppler, in which the area under the Doppler spectrum (excluding the DC component due to stationary tissue) is calculated [7]: ( ) where I' and Q' are high pass filtered from I and Q using an infinite-impulse-response (IIR) filter to remove signal from bulk tissue.The power Doppler signal is related to the volume of moving blood within the imaging volume and results in the loss of blood velocity and directionality information [7].Since in vivo bulk tissue motion can place the stationary tissue OCT signal outside the IIR filter (especially with the high velocity sensitivity associated with DOCT), V btm must be calculated from the Kasai estimator and velocity histogram.Then the IIR filter must be dynamically generated, different for each axial scan line, essentially creating a band-pass filter centered at V btm to remove the stationary tissue signal and provide an artifact-free power Doppler image.In OCT imaging, this approach increases the computation complexity considerably, in contrast to ultrasound where power Doppler requires less computation than color Doppler.Hence for this work, power Doppler OCT can only be calculated during post processing rather than in real-time.

Doppler spectrum mode
A clinically important display mode for Doppler ultrasound is the so-called sonogram or spectral display [7], essentially a joint time-frequency analysis of the flow signal using a transform method such as short-time Fast Fourier Transform (ST-FFT).The spectral display is usually calculated at a particular location within a blood vessel, often identified by color Doppler imaging, to illustrate the velocity (or Doppler frequency) distribution as a function of time.The Doppler spectrum of that location at time nT a (T a is the period for one axial scan, equivalent to PRF 1 ) is then calculated by FFT over a window length N FFT , which is selected using similar criteria as h B : where ( ) is the complex OCT signal sequence from time (n−N FFT )T a to n T a , M is the range gate size, and W is a window function.This display mode is especially useful for imaging pulsatile flow.Associated with the spectral display, an audio output of the Doppler frequency distribution is often presented to the physician in Doppler ultrasound [7].Audio presentation can also be accomplished in DOCT (Fig. 5).The audio signals are digital-toanalog converted at 8 kHz (matching the axial scan frequency of ~ 8.05 kHz) with 16-bit resolution in stereo, and are played along with the spectral display (Figs. 14 and 15 below).

Software
Detailed software processing algorithms were first implemented in LabVIEW (National Instrument, Austin, TX), operating in post-processing to validate the code.The LabVIEW program can be used for real-time imaging when frame rate is less than 2 fps.Subsequently, a multi-threaded C++ program utilizing DirectX display and sound subroutines (Microsoft, Redmond, WA) was developed for a personal computer with dual Xeon (Intel, CA) processors.The data processing and image display time was minimized by code optimization and balancing the computation load on each processor, which allowed real-time signal acquisition, processing, and display of the structural and flow images.
In summary, we have constructed a high speed DOCT system that can perform axial scans 2 ~ 8 times faster than previously reported, using a full implementation of a coherent demodulator.The digital signal processing technique is flexible so that the high axial scan frequency can provide high frame rate, high velocity sensitivity, or an optimized compromise between them.The combination of hardware and software processing strategy for velocity estimation allowed efficient computation for real-time operations.

Performance results and discussion
A key parameter describing the overall speed of an imaging system is the pixel-wise throughput rate.Implementing mature ultrasound signal processor designs, the current DOCT system achieved 4.1×10 6 pixels/sec at f a = 8.05 kHz and 6.6×10 6 pixels/sec at f a = 12.95 kHz, when digitizing only the forward scans.This system is comparable to the throughput rate of the smart detector array [5], but it is still compatible with catheter-based OCT designs to ensure high frame rate.At the present time, the software processing speed limits the efficiency to utilize all of the digitized data.With two 1.7 GHz Xeon processors, the software can process 3.1×10 6 pixels/sec (~75%).The efficiency can approach 100% when using two 2.4 GHz processors, resulting in few dropped frames.
Compared to the simplified auto-correlation method using a demodulating logarithmic amplifier with limiter output [3], the phase noise of the full implementation of coherent demodulator is considerably lower, allowing detection of slower blood flow.In addition, the range for velocity detection is doubled from ± f a /4 to ± f a /2.The decreased noise and the increased velocity detection range both contribute to a wider VDR.Furthermore, the demodulating logarithmic amplifier and its limiter output are also more sensitive to noise [3], limiting the SNR for both structural and flow imaging.
In the following sections we present example images and examine the performance of the DOCT system using stationary, steady flow, and pulsatile flow phantoms.The results presented here are obtained using the 8.05 kHz RSOD (imaging with the 12.95 kHz system will be reported separately).

Background noise in color Doppler and velocity variance modes
Even with no blood flow present, some background noise in the color Doppler and velocity variance images exists, arising from phase instability (Fig. 4), and physical factors such as Brownian motion in solutions.In addition, regions with low OCT signal SNR produce noisier velocity estimates (Fig. 7).With lateral scanning, especially at high frame rates, mechanical vibration of the sample arm fiber optics and inherent speckle modulation drastically increase the background noise.To evaluate this, we imaged a stationary target (0.5% Intralipid solution), and calculated the root-mean-square (RMS) values of the background Doppler frequency and normalized velocity variance at various frame rates and SNR levels (Figs. 6  and 7).The noise background in the Doppler frequency shift reaches a minimum, determined by the phase stability, when no lateral scanning is involved (i.e., 0 fps, or equivalent to M-mode ultrasound).At low frame rates (2 fps) and high SNRs (> 10 for Doppler frequency, and > 100 for velocity variance), the measured Doppler shift frequency background noise levels approached the minimum level.Slow (~ 100 µm/s), non-pulsatile blood flow can be imaged at 0 -2 fps, yielding a VDR ~ 40 dB.At intermediate frame rates (4 -8 fps), the background noise levels increased due to vibration and speckle modulation, decreasing the VDR to 30 ~ 32 dB.These conditions are probably adequate to image steady or low pulsatility (< 2 Hz) blood flow at 0.1 ~ 2 mm/s.Dynamic cardiac imaging demands high temporal resolution and fast frame rates (16 -32 fps); however, the much-reduced VDR (15 ~ 24 dB) means that only high speed flow can be reliably detected.Because speckle modulation is related to the lateral scanning speed, so that for a given frame rate a reduction in the lateral field of view can reduce the background noise.Finally, the velocity variance required much higher SNR than Doppler frequency to achieve the same background level, as expected since the variance estimation is affected more by the presence of OCT signal noise.In summary, Fig. 7 establishes the minimum detectable flow velocities at different frame rate and SNR conditions.

Steady flow velocity calibration in color Doppler and velocity variance modes
The accuracy of flow velocity measured by DOCT was verified using a calibrated infusion pump on a flow phantom.The flow phantom was in a plastic tube with inner diameter 0.75 mm through which 0.5% Intralipid solution was pumped (Harvard Apparatus, Holliston, MA).The flow was laminar, and a parabolic profile was assumed (and seen in color Doppler mode) such that the peak flow velocity could be calculated from the volumetric flow rate, set by the infusion pump.Imaging was performed at 0 fps (i.e., M-mode).The Doppler angle was 79.5° (81.4° after correction for refraction, assuming an refraction index of 1.33).The flow rate was adjusted from 10 µL/min to 3 mL/min, corresponding to peak velocities ranging from 0.7 mm/s to 22.6 cm/s.After accounting for Doppler angle, this corresponded to The accuracy of the color Doppler mode was assessed by comparing the peak velocity set by the flow pump and the measured velocity (Fig. 10).The results suggested that color Doppler mode DOCT could accurately detect flow velocities within the non-aliased velocity measurement range, which was about ± 1.9 mm/s.The increased scattering of data points with increasing velocity is due to the greater variance at higher flow velocities (Fig. 11), and is not due to the error in the infusion pump set velocity (< 0.5%).Measured velocity [mm/s] In color-Doppler mode, the non-aliased velocity measurement range is given by Eq. ( 7), and at 81.4° the maximum measurable velocity is ~ 1.3 cm/s, shown by the dashed line in Fig. 11.Using velocity variance as a measure of flow speed, DOCT could estimate flow speeds well beyond the Doppler aliasing range, in this case above 10 cm/s.The increase in flow speed detection range was more than 7-fold compared to color Doppler estimation at 81.4°.Because the aliasing range was Doppler angle dependent, this increase was also angle dependent and ranged from 2 to 50-fold.The inverted Gaussian fit in Fig. 11 can serve as a calibration curve for obtaining speed information from measured velocity variances.The gain in increased measurement range comes at the cost of loss in directionality and reduced accuracy in velocity estimation.Since the variance measurement does not distinguish flow direction, only the speed, not velocity, can be obtained.Due to the strong SNR dependence, the speed estimation by variance loses reliability when the SNR in the structural image is reduced (Fig. 7(b)).In addition, at high flow speeds (e.g., > 10 cm/s), the error in the speed estimation by variance can be as high as 100%.Hence, there is an upper limit in velocity (and speed) estimation by phase-shift measurements between axial scans, regardless of the method (color Doppler or velocity variance).This limit can be overcome using information along the axial scan, similar to ST-FFT methods [6] which estimate velocity from Doppler frequency shift of a single axial-scan.
One key objective in the DOCT development was to increase the frame rate while preserving the velocity detection sensitivity, such that dynamic flow profiles can be visualized.In the next section, we demonstrate DOCT imaging of flow phantoms under pulsatile conditions.

Pulsatile flow imaging by color Doppler, velocity variance, and Doppler spectrum modes
The discontinuity in stepping motor motion of the infusion pump (Fig. 8) was exploited to provide pulsatility in the flow phantom.Using larger step sizes at low repetition rate, reproducible pulsatile flow (peak velocities 1 -2 mm/s) was obtained (Fig. 12).Cross-sectional videos of the pulsatile flow phantom were acquired at various frame rates (Fig. 13).Once a region of interest in the flow (e.g., the center of the flow phantom) was identified by the color-Doppler mode, the pulsatility could be characterized in detail using the Doppler spectrum mode, as shown in Fig. 15.For comparison, we include the result from steady flow conditions in Fig. 14, with the flow rate set at 150 µL/min or 1.7 mm/s peak velocity at the center of the phantom.
The Doppler spectrum, indicating the velocity distribution at the center of the flow phantom, is plotted as a function of time in Fig. 14(a).The mean Doppler frequency was approximately 3.6 kHz, as expected for 1.7 mm/s flow velocity.The vertical spread of this 3.6 kHz signal is wider than the expected measurement error at this velocity (see Fig. 10(a)).This was due to the fast pulsatility of the flow as explained and confirmed by the vertical streaks in Fig. 14(b).The audio output, which can be calculated in real-time, does correspond to this artefact, since it is not a monotone.It should be pointed out that the Doppler spectrum contains more information than the color-Doppler mode, since the latter only displays the mean velocity, not the entire velocity distribution.This will be important in investigating high shear-rate flow conditions when large velocity gradients exist within the focal volume.
When pulsatility was deliberately introduced, Fig. 15 clearly demonstrated the velocity distribution variation over time for the ~9 Hz pulsatile flow.The vertical spread at any instance in time in Fig. 15(a) is smaller than that of Fig. 14(a), since the pulsatility was slow enough to be resolved.The audio output confirms this observation, since it sounds like a clear tone with varying frequency.We believe the combination of structural, color-Doppler, and Doppler spectrum modes can offer a valuable system for clinical use that would provide every step in microvascular imaging, including identification of the vessel, locating of a region of interest for temporal analysis, and calculating not only of the mean velocity, but also of the velocity distribution, including the pulsatile profile.

Comparison of power-Doppler and velocity variance imaging modes
A typical set of images measuring steady flow in a rectangular flow phantom, are shown in Fig. 16.The color-Doppler image contained several rings due to aliasing phase wrap-around.In contrast, these rings were completely removed in the power-Doppler image (although the directional flow information was also lost).To compare power Doppler and velocity variance imaging modes, a similar set of images including the velocity variance is shown in Fig. 17.While both the velocity variance and power-Doppler images remove the aliasing effect (both losing the directional information), there is one important difference exists between them.The variance image is brighter in the center, corresponding to higher velocity flow (Fig. 17(c)).This gradient is not present in the power-Doppler image (Fig. 16(c)), since it does not measure flow velocity but provide a "flow" versus "no-flow" determination by thresholding.
The power-Doppler image in Fig. 16 was processed using a high-pass IIR filter centered at zero velocity, since no bulk tissue motion was present.The IIR filter was identical for each axial scan, and provided a sharp cut-off at ~ 500 Hz.As explained earlier, if motion artifacts are expected, a band-pass IIR filter (probably different for each axial scan) must be used, which will significantly reduce the processing speed.

Conclusions
In conclusion, we have demonstrated a DOCT system capable of high speed imaging with wide velocity dynamic range using a full implementation of hardware coherent demodulation.
A number of software processors for blood flow visualization, including color-Doppler, velocity variance, power-Doppler, and Doppler spectrum (including audio output) were demonstrated.Real-time imaging, at 16 fps with 500 × 512 pixels per frame and 32 fps with 250 × 512 pixels per frame, allowed visualization of 9 Hz pulsatile flow without motion artifacts.The wide velocity dynamic range of the system, covering from 2 µm/s to 10 cm/s, allowed detection of bi-directional flow up to ± 1.9 mm/s along the optical axis in color-Doppler mode.The minimum detectable velocity was ± 17 µm/s when 8 axial-scans were averaged, and ± 2 µm/s when 128 axial scans were used.When used in velocity variance mode, the upper limit for flow speed detection was more than 10 cm/s, significantly higher than previous reports using phase-shift methods [2,3].

Fig. 2 .
Fig. 2. Signal and filter frequency response during coherent demodulation.Solid line: Calculated frequency response for 8.05 kHz RSOD when scanning 1.3 mm of depth.Dashed line: Frequency response of 12.95 kHz RSOD when scanning 1.0 mm of depth.After mixing with the 4.3 MHz carrier frequency, the signals are down-shifted to the base-band, and upshifted to 8.6 MHz.Blue line: Low pass filter response, approximating a matched filter for the base-band signals, and completely removing the up-shifted signal around 8.6 MHz.The I and Q signals are simultaneously digitized at ~ 10 MHz (for f a = 8.05 kHz) or ~ 16 MHz (for f a = 12.95 kHz) with 12-bit accuracy.The complex base-band OCT signal can be then expressed as jQ I S + = .(2) To increase data throughput from the A/D card to computer memory, an A/D conversion clock pulse train triggered by the RSOD is used, with data acquisition occurring during the central ~ 80% of the depth scan.The acquired signals can be expressed as 2D arrays: I m,n and Q m,n with m and n denoting the indices in the depth and lateral directions, respectively.These signals are transferred by direct memory access without computer CPU intervention, reducing the computational overhead.

Fig. 3 .
Fig. 3. (a) The mean phase shift can be calculated by evaluating the phase angle of the individual OCT signals and then computing the phase differences between axial scans.(b) The result is equivalent to computing the phase angle of the mean (〈X〉,〈Y〉) vector.The computation in (b) is numerically less complex than in (a), and so (b) is more suitable for realtime processing.

Fig. 4 .
Fig. 4. (a) Spatial averaging masks for improving the SNR in structural and flow images.Shaded areas indicate regions with shared computation between pixels.(b) Calculated axial velocity using two-quadrant arctan (red) and four-quadrant arctan (blue) functions, showing the aliasing effect and the color map used in color Doppler OCT.(c) Measured phase noise processed from different ensemble lengths (N = 2 to 128, and M = 1).Each data set (dots) represents the normalized distribution of phase noise measured from 20,000 pixels using a stationary diffuse reflectance target, with a Gaussian fit (line).
h v ∆ is the velocity step size.The choice of h B is based on two considerations: (a) h B should be a power of 2 to allow binary tree sorting for computing the histogram efficiently, and (b) the velocity step size should match the expected velocity resolution: (C) 2003 OSA 7 April 2003 / Vol.11, No. 7 / OPTICS EXPRESS 801

Fig. 5 .
Fig. 5. Schematic for generating audio output to accompany the Doppler spectrum display.The directional information is encoded in the stereo's left and right channels.H: digital Hilbert transform performing 90° phase shift.DAC: digital-to-analog convertor.

Fig. 7 .
Fig.7.SNR varied with depth; therefore, the image was divided into 5 horizontal regions of interest (ROI) to investigate velocity noise as a function of SNR.Each data point is the mean value within a ROI containing 25,000 pixels.(a) Background noise levels in the color Doppler mode at 2 -32 fps and the corresponding structural image SNR conditions.(b) Background noise levels in the velocity variance mode at various frame rates and corresponding structural image SNR conditions.Note: 0 fps indicates no lateral scanning, i.e., analogous to M-mode ultrasound.
-3.4 cm/s.Two examples are shown in Figs. 8 and 9, illustrating low flow rate (not aliased) and high flow rate (aliased) conditions, respectively.

Fig. 8 .
Fig. 8. (a) Structural, (b) color Doppler, and (c) normalized velocity variance images of 0.5% Intralipid flow phantom, acquired at 0 fps.At 0 fps, the x-axis represents time, not distance.Flow rate = 60 µL/min, corresponding to 0.68 mm/s peak velocity after accounting for Doppler angle.Vertical streaks in (b) are due to the discrete stepping motor motion on the infusion pump, prominent at low flow rates.Scale bar = 0.5 sec.

Fig. 9 .
Fig. 9. (a) Structural, (b) color Doppler, and (c) normalized velocity variance images of 0.5% Intralipid flow phantom, acquired at 0 fps.At 0 fps, the x-axis represents time, not distance.Flow rate = 1.5 mL/min, corresponding to 1.7 cm/s peak velocity after accounting for Doppler angle.Notice the aliasing pattern in (b), and the disappearance of the vertical streaks.Scale bar = 0.5 sec.

Fig. 10 .Fig. 11 .
Fig. 10.Comparison of measured peak velocity in color Doppler mode (measuring the velocity vector along the optical axis) and the expected flow velocity, as set on the infusion pump.Each data point is the mean of 25,000 pixels at the center of the phantom.(a) Low flow rate conditions without aliasing.Solid line: Unity slope.Error bars: Standard deviation within the 25,000 pixels.(b) With flow velocities greater than 1.9 mm/s, aliasing occurs.Solid line: Expected relationship, following Fig. 4(b).Notice the increased scattering of data points with higher velocity, which forms the basis of velocity variance imaging.

Fig. 12 .
Fig. 12.(a) Structural, (b) color Doppler, and (c) normalized velocity variance images of a 0.5% Intralipid flow phantom, acquired at 0 fps.The x-axis is time and flow goes from right to left.The flow rate was adjusted to be pulsatile at ~ 9 Hz, with the peak velocity just below the aliasing velocity limit.Scale bar = 250 ms.

Fig. 13 .
Fig. 13.(77, 172, 359, 737 and 1423 kB, respectively) Left to right: 2, 4, 8, 16, and 32 fps videos of the pulsatile flow phantom, each containing structural (top), color Doppler (mid), and velocity variance (bottom) images.Each video is 2 seconds long.The diameter of the tube is ~ 0.75 mm.Notice the peak velocity drifting laterally in the 2 and 4 fps videos, which is an artifact of low frame rates.The slow pulsatility (~ 1 Hz) observed in the 8 fps video is an aliasing artifact, since the actual pulse rate is ~ 9 Hz.At 16 fps, the pulsatility can be visualized properly.At 32 fps, although the temporal resolution is improved, the reduced velocity sensitivity degrades the image quality markedly.

Fig. 14 .
Fig. 14.Click to hear the associated audio output [33 kB].(a) Doppler spectrum and (b) color Doppler images of the flow phantom driven at a steady flow rate of 150 µL/min, corresponding to 1.7 mm/s peak velocity after Doppler angle correction.The Doppler spectrum is obtained at the center of the phantom, along the white line shown in (b).

Fig. 15 .
Fig. 15.Click to hear the associated audio output [33 kB].(a) Doppler spectrum and (b) color Doppler images of the flow phantom driven under pulsatile conditions.The ~ 9 Hz pulses had peak flow rates ~ 2 mm/s.Notice the small reflection observed in (a), likely due to bias and gain mismatch in the I and Q channels.

Fig. 16 .
Fig. 16.(a) Structural, (b) color-Doppler, and (c) normalized power Doppler OCT images of a glass channel flow phantom of 0.25% Intralipid with peak velocity ~ 6 mm/s.Note the aliasing effect in (b) due to phase wrap-around, and the loss of flow directionality in (c).

Table 1 .
Requirements of axial scan frequency and analog-to-digital sampling rates for 500 × 500 pixels per frame.

Table 2 .
Computation complexity for velocity estimation for a single pixel

Table 3 .
Ensemble length and overlap values for various frame rates.