*in vivo*.

## 1.

## Introduction

Retinal metabolic rate of oxygen (${\mathrm{rMRO}}_{2}$) is a parameter critical for the fundamental investigation and clinical diagnosis of several blinding diseases such as age-related macular degeneration, glaucoma and, most importantly, diabetic retinopathy (DR).^{1}2.^{–}^{3} Previous studies suggested that retinal metabolic abnormalities play an important role in the development of DR^{2}^{,}^{4}^{,}^{5} and that ${\mathrm{rMRO}}_{2}$ is a potential biomarker for early diagnosis of DR before irreversible damage occurs.^{6}7.^{–}^{8} In order to calculate ${\mathrm{rMRO}}_{2}$, we need to measure oxygen saturation (${\mathrm{sO}}_{2}$) and the flow rate in major retinal blood vessels, both of which can be effectively quantified simultaneously by visible-light optical coherence tomography (vis-OCT).^{9}10.^{–}^{11} As a single imaging modality for ${\mathrm{rMRO}}_{2}$ measurement, vis-OCT detects ${\mathrm{sO}}_{2}$ by retrieving the absorption spectrum of blood from spectroscopic analysis and measures the flow rate using the established Doppler OCT method.^{12}13.^{–}^{14} Compared with conventional OCT using near-infrared (NIR) light, vis-OCT provides reliable ${\mathrm{sO}}_{2}$ measurement and higher resolution, but has reduced measurement range of flow velocity.^{13}

This drawback of vis-OCT in flow measurement originates from the underlying mechanism of Doppler OCT, which relates phase shift between sequentially acquired OCT A-lines to the axial displacement of scatterers. The axial displacement is converted to flow velocity based on the A-line rate of system and Doppler angle, which is the angle between the flow direction and light probing axis.^{15}16.17.^{–}^{18} Since phase variation is periodic, the actual detectable Doppler phase shift always lies between $-\pi $ and $\pi $. Phase wrapping occurs when flow induced axial displacement exceeds the center wavelength of the probing beam.^{15}^{,}^{17} Therefore, using shorter wavelengths, vis-OCT is more sensitive to low flow velocity, but more susceptible to phase wrapping in measuring high flow velocity.

There are a few potential ways to extend the measurable velocity range of vis-OCT. First, phase unwrapping can be applied to increase the range beyond $-\pi $ and $\pi $. However, although various unwrapping techniques have been investigated, it remains challenging to remove artifacts induced by phase wrapping in clinical settings.^{19}^{,}^{20} Second, the phase shift induced by a certain flow velocity can be reduced by decreasing the time interval between two sequential observations, which is the inverse of A-line acquisition rate. This can be done by decreasing the spectrometer integration time of spectral-domain OCT (SD-OCT) or by switching to a fast swept-source system. However, a reduced spectrometer integration time leads to a lower signal-to-noise ratio (SNR) and a swept-source for visible light is not yet available. Third, the imaging operator can intentionally position the OCT probing beam and retinal blood vessel for a larger Doppler angle to decrease the projected axial velocity. However, this procedure can significantly increase the operation complexity and may be impractical for the region around the optic nerve head.

We seek to extend the flow velocity measurement range for vis-OCT without using sophisticated phase unwrapping. In most existing Doppler OCT algorithms, the center wavelength or wavenumber of the broadband light source is used to translate the detected phase shift to the axial displacement of scatters.^{15}^{,}^{21} However, if we split the broad spectrum into a number of narrower subbands for OCT reconstruction, a series of Doppler phase shifts can be obtained and these phase shifts follow a linear relationship with the center wavenumbers of the subbands.^{12}^{,}^{22}^{,}^{23} The slope of this linear function, instead of the actual phase shift values, can be used to retrieve the axial flow velocity of the scatterers, which is independent of phase wrapping.^{15}^{,}^{21} Hence, spectroscopic Doppler analysis may extend the maximum measurement limit beyond $-\pi $ and $\pi $, which is especially beneficial in vis-OCT.

In this paper, we report a spectroscopic Doppler analysis method to quantify flow velocity in spectral-domain vis-OCT. We performed a series of short-time Fourier transforms (STFT) of the original broadband OCT interferogram and obtained Doppler phase shifts from different subbands. We applied linear regression on the wavenumber-dependent phase shift values and retrieved the axial flow velocity from the fitted slope. We validated our algorithm using a vis-OCT system with a 100-nm bandwidth at a center wavelength of 560 nm. Phantom experiments showed that both conventional Doppler and the spectroscopic Doppler analysis methods measured flow velocities accurately when the velocity was low; but spectroscopic Doppler analysis increased the measurable range of the flow velocity. We also demonstrated the spectroscopic Doppler analysis method in measuring *in vivo* retinal blood flow in rodents.

## 2.

## Method of Spectroscopic Doppler Analysis

## 2.1.

### Principle of Doppler Optical Coherence Tomography

SD-OCT measures flow velocity by detecting the phase difference between two sequential A-lines.^{21} If a single moving scatterer appears at depths of ${z}_{0}$ and ${z}_{0}+\delta z$ in the first and second A-lines, we can describe the corresponding interferograms $i(k)$ of the two measurements as

## (2)

$${i}_{1}(k)=\xi s(k)\sqrt{{R}_{R}{R}_{S}}\text{\hspace{0.17em}}\mathrm{cos}[2nk({z}_{0}+\delta z)],$$## (4)

$${I}_{1}(z)=\xi \sqrt{{R}_{R}{R}_{S}}[S(z-2n{z}_{0}){e}^{j2n{k}_{\mathrm{equ}}\delta z}+S(z+2n{z}_{0}){e}^{-j2n{k}_{\mathrm{equ}}\delta z}],$$For a light source with a symmetric spectrum, center wavenumbers ${k}_{c}$ and ${k}_{\mathrm{equ}}$ are interchangeable in Eq. (6). However, in real cases when the source spectrum is asymmetric, using ${k}_{c}$ directly leads to bias in flow velocity measurement.

We use $T$ to represent the acquisition time interval between two sequential A-lines, which is usually determined by the spectrometer integration time in SD-OCT or the wavelength sweeping cycle in a swept-source OCT (SS-OCT). We can rearrange Eq. (5) and obtain the axially projected flow velocity ${v}_{\parallel}$ by taking the ratio of the axial displacement $\delta z$ and $T$

Phase wrapping happens when the axial displacement induced absolute phase shift exceeds $\pi $, which results in ambiguous axial flow velocity measurement. The phase wrapping limited maximum absolute axial velocity (PLV) is given by^{15}^{,}^{18}^{,}^{21}

Another factor that limits the flow detection is the washout effect, which describes the signal loss due to averaging of motion-varied interferometric fringes over the integration time of the optical detector. For SD-OCT, the washout effect starts at around PLV and will cause a Doppler phase signal unrecoverable beyond twice PLV.^{19}

## 2.2.

### Principle of Spectroscopic Doppler Optical Coherence Tomography

STFT has been used in spectroscopic analysis of OCT signals to extract wavenumber-dependent absorption and scattering properties in tissue,^{24}25.^{–}^{26} which can be described as

Therefore, we can use spectroscopic analysis and linear fitting to obtain $a$ and, therefore, ${v}_{\parallel}$.

We performed numerical simulation to better demonstrate how Doppler phase shift varies with respect to center wavenumber of the OCT subband. Though OCT reconstruction and Doppler analysis are achieved in wavenumber domain, researchers often refer to OCT bandwidth in wavelength. Therefore, we provided both wavenumber and wavelength values in the following description. The simulated OCT source had a square-shaped spectrum with $(0.19\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (100 nm) bandwidth and $(1.10\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (570 nm) center wavelength. The Gaussian windows used in STFT had $(0.02\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (9 to 12 nm, depending on the wavelength) standard deviation, and the center wavenumbers ranged from $(1.04\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (604 nm) to $(1.16\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (542 nm) at $0.03\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ interval, as shown in Fig. 1(a). In Sec. 2.1, we use ${z}_{0}$ and ${z}_{0}+\delta z$ to represent the depth locations of a moving scatter in the two sequential OCT A-lines and $\delta z$ is the axial displacement. In the simulation, we simulated seven flow velocities whose axial displacements between two sequential A-lines were $0,\delta {z}_{0},2\delta {z}_{0},\dots ,6\delta {z}_{0}$, where $\delta {z}_{0}$ was 35.7 nm. According to Eq. (5), $6\delta {z}_{0}$ would lead to a $2\pi $ Doppler phase shift assuming that the center wavelength was 570 nm and that the refractive index was 1.33. Figure 1(b) shows the subband-dependent Doppler phase shifts under different axial displacements, where the colored dots indicate the measured phase shifts and the dashed black lines show the fitted linear relationship. Each color represents one of the OCT subband in Fig. 1(a). One would note that 0 and $6\delta {z}_{0}$ displacements would both yield zero flow velocity if measuring the Doppler phase shift with the whole spectrum due to phase wrapping. However, the slope of subband-dependent phase shift array is uniquely related to displacement values, unaffected by phase wrapping. Another interesting phenomenon occurs in the case of $3\delta {z}_{0}$, where a $2\pi $ shift exists in the B3 subband. This is because the displacement induced phase shifts are so close to $\pi $ that phase wrapping occurs among different subbands. We call this phenomenon spectroscopic phase shift discontinuity, which needs to be corrected before linearly fitting different phase shift values with subband center wavenumbers. Although this phase shift discontinuity can also happen when flow induced phase shifts are around $3\pi $, $5\pi $, and etc., washout effect limits the measurable axial flow velocity within twice PLV, corresponding to $2\pi $ Doppler phase shift. Therefore, we only need to correct the spectroscopic phase shift discontinuity that happens around $\pi $.

## 2.3.

### Flow Velocity Calculation in Spectroscopic Doppler Optical Coherence Tomography

The major application of Doppler OCT is to measure blood flow velocity by averaging Doppler phase shift values of all pixels within a blood vessel cross section. In cases where the flow velocity exceeds PLV, phase unwrapping needs to be implemented in conventional Doppler OCT to obtain correct the measurement result. If the flow velocity is constant throughout the blood vessel, phase unwrapping can be done simply by adding or subtracting $2\pi $ to all the pixels. However, the actual flow profile within the vessel follows an approximate three-dimensional paraboloid. When only the flow velocity in the center of a blood vessel exceeds PLV, it is challenging to find the boundary between the wrapped and nonwrapped areas. The state-of-the-art two-dimensional algorithms to correct this type of phase wrapping are usually iterative processes and often generate biased offset error in noisy datasets.^{27}28.^{–}^{29}

In spectroscopic Doppler OCT, we can retrieve flow velocity based on the linear relationship between phase shifts and center wavenumbers of OCT subbands. Phase wrapping only affects the intercept of the fitted line without changing its slope. Since we extract the flow velocity from the slope of the fitted line, phase unwrapping is not required. But we do have to correct the spectroscopic phase shift discontinuity, which is demonstrated in Fig. 1(b). We summarize the major steps of spectroscopic Doppler analysis in Fig. 2.

Step 1: STFT

We perform STFT on the broadband OCT B-scan and obtain a set of subband B-scans. We calculate phase shift values in each subband image in the same way as conventional Doppler OCT processing. We segment the region within blood vessel from one subband phase shift image and apply the obtained boundary to all subbands.

Step 2: Correction for phase shift discontinuity

In this step, we aim to calculate the average phase shifts within the blood vessel for each subband. If there is no spectroscopic phase shift discontinuity, the average phase shift ${\overline{\mathrm{\Delta}\phi}}_{D}(\tau )$ can be calculated as follows:

where $\tau $ is the center wavenumber of subband used in STFT, $N$ is the number of pixels inside a blood vessel used for spatial averaging, and $(x,z)$ is the coordinate of the pixel in the B-scan image.## (12)

$${\overline{\mathrm{\Delta}\phi}}_{D}(\tau )=\frac{\sum \mathrm{\Delta}{\phi}_{D}[\tau ,(x,z)]}{N},$$However, if the phase shift values of some pixels are close to $\pi $, spectroscopic phase shift discontinuity may occur. To preserve the linear relationship, the average phase shift should be calculated as follows instead

where a correction term is added to compensate for the discontinuity. To calculate the correction term, we select a phase shift image reconstructed on an arbitrary subband ${\tau}_{0}$ as the reference and subtract the reference image from phase shift images of all the subbands as## (13)

$${\overline{\mathrm{\Delta}\phi}}_{D}(\tau )=\frac{\sum \mathrm{\Delta}{\phi}_{D}[\tau ,(x,z)]}{N}+\frac{\sum \mathrm{\Delta}{\phi}_{\mathrm{corr}}[\tau ,(x,z)]}{N},$$## (14)

$$\psi [\tau ,(x,z)]=|\mathrm{\Delta}{\phi}_{D}[\tau ,(x,z)]-\mathrm{\Delta}{\phi}_{D}[{\tau}_{0},(x,z)]|.$$$\psi [\tau ,(x,z)]$ is expected to be small if no discontinuity occurs. Therefore, we compensate the discontinuity by adding or subtracting $2\pi $ when $\psi [\tau ,(x,z)]$ exceeds a threshold value

where $\mathrm{sgn}(\xb7)$ is the sign function. We mentioned in Sec. 2.2 that the spectroscopic phase shift discontinuity only happens when the flow induced phase shift is around $\pi $ and that the maximum measurable phase shift should be smaller than $2\pi $ due to the washout effect. Therefore, we define the threshold value th as## (15)

$$\mathrm{\Delta}{\phi}_{\mathrm{corr}}[\tau ,(x,z)]=\{\begin{array}{cc}2\pi \xb7\mathrm{sgn}\{\mathrm{\Delta}{\phi}_{D}[{\tau}_{0},(x,z)]\},& \psi [\tau ,(x,z)]>\mathrm{th},\\ 0,& \psi [\tau ,(x,z)]\le \mathrm{th},\end{array}$$Step 3: Linear fitting

We do linear regression to the subband-dependent average phase shift values using least-squares fit and extract the slope of the fitted linear relationship.

Then, the axial flow velocity is calculated by Eq. (11).

## 3.

## Experimental System

We demonstrated flow velocity measurement in vis-OCT using spectroscopic Doppler analysis. The system schematic is shown in Fig. 3(a). We used a broadband supercontinuum laser (Superk EXTREME, NKT Photonics) to generate a 90-nm bandwidth visible light source centered at 565 nm, which is also used in clinical settings.^{11}^{,}^{30} The source beam is coupled into a $2\times 2$ fiber coupler (Nufern 460-HP, 50:50 splitting ratio, GouldFiber Optics) which delivers the light to a sample arm and a reference arm. In the sample arm, the illumination light is collimated and scanned by a set of galvanometer mirrors (6210H, Cambridge Technology). In phantom experiments, we focused the illumination beam using an achromatic doublet (Edmund Optics) with a focal length of 35 mm [Fig. 3(b)]. In *in vivo* experiments, we used a pair of achromatic doublets with focal lengths of 75 and 15 mm, respectively, to relay the illumination beam onto the cornea [Fig. 3(c)]. In the reference arm, after collimation, the beam propagates through a BK7 glass block for dispersion compensation and is reflected by a mirror. The backscattered light from the sample arm interferes with the back-reflected light from the reference arm in the fiber coupler and the interferogram is collected by a home-built spectrometer, consisting of a transmission grating ($1800\text{\hspace{0.17em}\hspace{0.17em}}\text{line}/\mathrm{mm}$, Wasatch Photonics), an SLR lens (85 mm, $\mathrm{f}/1.4$, Samyang), and a line-scan camera (spL4096-70 km, Basler). Though there are 4096 linear pixels in the camera, we only use the center 2048 for spectrum acquisition.

## 4.

## Experimental Results

## 4.1.

### Phantom Validation

We imaged a flow phantom made from 0.5% intralipid (Sigma, I141-100ML) in a plastic tube with a $125\text{-}\mu \mathrm{m}$ inner diameter. The refractive index of intralipid is 1.34.^{31} Flow rate was controlled by a motorized syringe pump (Model A99-EM, Razel Scientific Instruments) within a speed ranging from 0 to $0.25\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$. The Doppler angle of the tube was 83.5 deg. The OCT A-line rate was 20 kHz. We acquired 2-mm wide B-scans crossing the plastic tube, which consisted of 2048 A-lines. In each measurement, we performed 16 repeated B-scans for averaging. Note that in the phantom experiment, we used flow rate instead of flow velocity because flow rate, containing both velocity and cross-sectional area of flow, is of great interest in hemodynamic and metabolic studies. Using flow rate is also more convenient since the true flow rate values can be obtained from the syringe pump without the need to accurately measure the diameter of plastic tube. Therefore, we converted the flow velocity provided by Doppler OCT to flow rate in the following phantom study.

We measured flow induced phase shift using both conventional and spectroscopic Doppler analyses. In conventional Doppler analysis, we obtained the averaged phase shift inside the tube both with and without phase unwrapping. Before phase unwrapping, we further filtered the averaged phase shift images with a $3\times 15$ ($\text{axial}\times \text{lateral}$) pixels Doppler filter following the method reported in Ref. 22, since phase unwrapping is sensitive to noise. We selected an area close to the outer boundary of the flow area on the phase shift B-scan as a reference, where no phase wrapping occurs due to its low flow velocity. Then, we searched throughout the flow area and labeled the phase wrapped pixels if their values changed more than $\pi $ compared with adjacent pixels. Finally, we added or subtracted $2\pi $ to the labeled pixels and recovered the unwrapped phase shift values.^{27} In spectroscopic Doppler analysis, we generated 16 Gaussian windows for STFT, whose center wavenumbers were evenly distributed from $(1.06\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (593 nm) to $(1.16\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (542 nm). The windows had a standard deviation of $(0.02\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (10 nm at 565-nm wavelength). We applied STFT, corrected for possible spectroscopic phase discontinuity, and calculated subband-dependent phase shifts.

Figure 4(a) shows phase shift B-scans under two flow rates before and after phase unwrapping. There is no phase wrapping in flow 1 ($0.021\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$) and obvious phase wrapping in flow 2 ($0.218\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$). Therefore, before phase unwrapping, the averaged phase shift values throughout the flow area are similar between flow 1 and flow 2, which can be shown in Fig. 4(b) by the vertical positions of linear functions (data 1 and data 2). However, there is a clear difference between the slopes of the two linear functions. After phase unwrapping, the difference in flow rates between flow 1 and flow 2 can be better appreciated in Fig. 4(b) from both the vertical positions and the slopes of the linear functions (data 3 and data 4). The slopes of linear functions for data 2 and data 4 are the same, indicating that spectroscopic Doppler analysis can retrieve flow velocity without phase unwrapping.

To further compare conventional and spectroscopic Doppler analyses, we used both methods to measure seven different flow rates in phantom without phase unwrapping, as shown in Fig. 5. To increase flow measurement accuracy, we applied spatial averaging across the entire flow cross section in both methods, which provided a fair comparison. The flow rates ranged from 0 to $0.24\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$ with a step size of $0.04\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$. Phase wrapping occurred in flow rates greater than $0.12\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$. Figure 5(a) shows subband-dependent phase shifts at different flow rates. In order to avoid crossing among different phase shift lines and allow better illustration of their slopes, we manually offset these lines according to the preset flow rates in Fig. 5(a). Therefore, the intercepts of different lines do not represent real calculated values. Figure 5(b) compares measured flow rates using both methods. Since no phase wrapping exists in the four lower flow rates, both methods can accurately retrieve the preset values. For the three higher rates, conventional Doppler OCT fails to provide accurate measurement results due to phase wrapping. However, spectroscopic Doppler OCT is still able to quantify 0.18 and $0.21\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$. In the case of $0.24\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{L}/\mathrm{s}$, a severe washout effect occurs in the center of the flow area, which decreases the SNR. Therefore, the subband-dependent phase shift array is unstable, and its fitted slope value fails to provide an accurate estimate for the flow rate.

## 4.2.

### In Vivo Retinal Blood Flow Measurements

We used both conventional and spectroscopic Doppler analyses to compare measured retinal blood flow velocity without phase unwrapping in rats (Long Evans, 300 g, Charles River) *in vivo*. The details of animal preparation and imaging protocol have been described elsewhere.^{32} All experimental procedures complied with the Association for Research in Vision and Ophthalmology statement for the use of animals in ophthalmic and vision research. The laboratory animal study protocol was approved by the Institutional Animal Care and Use Committee at Northwestern University.

In brief, vis-OCT measurements were conducted in rat eyes using a circular scanning protocol.^{33}^{,}^{34} The diameter of the scanning circle was $\sim 0.5\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{mm}$ on the rat retina. Each circular scan contained 4096 A-lines. The vis-OCT probing power was 0.8 mW on the cornea. We first used conventional Doppler analysis, where we calculated phase shift images using the whole broadband spectrum. Then, we randomly selected 24 vessels whose phase shifts were roughly distributed between 0.2 and 1 rad. We applied spectroscopic analysis to these selected vessels, where we used the same set of Gaussian windows as in the phantom experiment, but adjusted the standard deviation to $0.03\times {10}^{5}\text{\hspace{0.17em}\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ (around 15 nm at 565-nm wavelength). We applied STFT on the raw OCT interferograms and obtained a series of phase shift images. We manually segmented the flow areas within retinal blood vessels from subband phase shift images. We obtained axial flow velocities, converted them to phase shift values, and compared them with the results from the conventional Doppler method. There are two reasons why we did not provide flow rates as in the phantom study. First, in the phantom study, since the true flow rate is continuously controllable, it is convenient to compare between the preset flow rate values and the measured values. Second, in the animal experiment, the flow rates across different vessels may be similar. However, variations in the Doppler angles among different blood vessels give rise to distinct axial flow velocities. Therefore, by using axial flow velocity without calculating the flow rate, we can compare the performance of conventional and spectroscopic Doppler methods at a wide range of phase shift values.

Figure 6 compares the rat retinal axial flow velocities represented by Doppler phase shifts using conventional and spectroscopic Doppler analyses without phase unwrapping. The scatterplot shows a good correlation between phase shifts measured by two methods when the mean Doppler phase shift is smaller than 1 rad. However, there is a larger variation as compared with the phantom results shown in Fig. 5 due to a reduced SNR. We can also observe in Fig. 6 that the five highlighted data points all lie below the diagonal line, suggesting lower flow velocity estimation by the conventional Doppler analysis. This is consistent with the phantom experiment that indicates that the conventional Doppler OCT method tends to underestimate flow velocity when phase wrapping occurs.

## 5.

## Discussion and Conclusion

We developed a spectroscopic Doppler analysis for flow velocity measurement. Instead of calculating Doppler phase shift using broadband OCT interferogram as in conventional Doppler OCT, we used STFT to reconstruct a series of complex OCT B-scans in different subbands, whose Doppler phase shifts were then calculated separately. Instead of averaging the phase shift values for improved SNR as proposed previously by Tan et al.,^{35} we fit these phase shift values against the center wavenumbers of the subbands using a linear function whose slope was related to the axial flow velocity. We demonstrated that spectroscopic Doppler analysis was unaffected by phase wrapping since the absolute shift values were not directly used in velocity calculation. Therefore, spectroscopic Doppler analysis does not rely on phase unwrapping to measure a vascular flow whose center velocity exceeds PLV. Using subbands with reduced bandwidth will decrease the axial resolution in the Doppler phase shift images, which affects measurement accuracy if the flow area cannot be clearly resolved at the reduced axial resolution.

Compared with conventional Doppler OCT and its phase unwrapping techniques, spectroscopic Doppler analysis may be a good alternative for flow velocity measurement in vis-OCT. Working on shorter wavelengths, vis-OCT has much improved axial resolution, but is more susceptible to phase wrapping. Due to higher resolution, the phase shift images of vis-OCT retain sufficient information even when only reconstructed on a spectral subband. Due to a much lower PLV, vis-OCT requires more careful phase unwrapping procedures when measuring flow velocity than using a conventional Doppler OCT method. Though phase unwrapping can provide satisfactory results when the image SNR is guaranteed, as shown in our phantom experiment, it is not trivial, and is even biased if the image quality is moderate, especially in clinical settings. In addition, the unwrapping errors may accumulate during the iterative process adopted in most unwrapping algorithms.^{27}28.^{–}^{29} Therefore, spectroscopic analysis is especially helpful in vis-OCT, though it can also be applied to conventional NIR OCT as well. However, we expect that spectroscopic analysis may be less valuable for OCT systems working beyond $1\text{\hspace{0.17em}\hspace{0.17em}}\mu \mathrm{m}$ because (1) these OCT systems have a much higher PLV and can depend on a conventional Doppler method and (2) the phase shift values do not vary much with the center wavenumber of different subbands and the fitted slope value may be more vulnerable to noise.

The total number and bandwidth of the spectral subbands affect the flow measurement accuracy. A limited number of subbands leads to fewer points for linear regression, which makes the fitted slope susceptible to noise. Too many subbands, on the other hand, requires more intensive computation with a marginal benefit. Therefore, we selected 16 subbands to perform spectroscopic analysis, which was a tradeoff between accuracy and computation efforts. We further confirmed our selection using numerical simulation. Selecting proper bandwidth for subbands is also important. A large bandwidth decreases the spectral fitting range and a small bandwidth leads not only to reduced resolution, but also lower SNR.^{36} In our experiment, we empirically used $(0.04\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ for phantom experiments and $(0.06\times {10}^{5})\text{\hspace{0.17em}}\mathrm{rad}/\mathrm{cm}$ for animal experiments. (Note that the subband width is twice the standard deviation of the Gaussian window used in STFT.) We used a larger bandwidth in animal experiments to limit the phase noise,^{36} since the SNR of *in vivo* data is naturally lower than that of the phantom data. A more thorough investigation will be necessary to provide comprehensive evaluation of the effects of subband selection for spectroscopic Doppler analysis.

Though spectroscopic Doppler analysis is robust to phase wrapping, it is still limited by the fringe washout effect in spectrometer-based SD-OCT. The subband-dependent phase shift becomes noisy when significant washout occurs and the slope extracted from linear fitting fails to reveal the flow velocity. SS-OCT is almost unaffected by the fringe washout effect since it increases the washout threshold by a factor that equals the total number of spectral samples in each A-line.^{19} Therefore, although not tested here, our spectroscopic Doppler analysis can potentially extend the velocity measurement range in SS-OCTs.

We notice that in the current study, sufficient SNR is critical to obtain reliable flow velocity using spectroscopic Doppler analysis. To guarantee sufficient SNR, we applied spatial averaging across the entire flow cross section in both phantom and animal studies. Since the flow velocity is translated from the slope of the fitted straight line, accurate measurement requires consistency of phase shift values calculated from spectral subbands, which themselves could be more susceptible to influence from noise due to reduced spectral bandwidth. Therefore, our future work will include comprehensively investigating the validity of spectroscopic Doppler analysis under different SNR conditions.

In summary, we present a spectroscopic Doppler analysis for flow velocity measurement in vis-OCT and tested it in both phantom and *in vivo* experiments. This method uses spectroscopic information of each phase-shifted image pixel to unambiguously retrieve flow velocity instead of depending solely upon the pixel spatial relationship to correct the phase wrapping as in the conventional Doppler OCT method. Though we demonstrated our method using spectrometer-based OCT only, we expect that SS-OCT will benefit from this method as well.

## Disclosures

Hao F. Zhang, and Wenzhong Liu have financial interests in Optician Health, which did not support this work.

## Acknowledgments

We acknowledge the generous financial support from the National Institutes of Health (NIH) Grant Nos. DP3DK108248, R24EY022883, R01EY026078, and R21EY027502.