## Abstract

Analyzing the experimental data of the velocity distribution in a fluid flow using Doppler Optical Coherence Tomography (OCT), we compared the Wigner distribution method to the short-time Fourier transform method, the Hilbert-based phase-resolved method and the autocorrelation method. We conclude that the pseudo Wigner-distribution signal processing method is overall more precise than other often-used methods in Doppler OCT for the analysis of cross-sectional velocity distributions.

© 2007 Optical Society of America

## 1. Introduction

Optical Coherence Tomography (OCT) is an imaging technique that appeared about 15 years ago [1]. It can provide in-depth scattering profiles of tissues with a resolution as high as one micron in a non-contact and quasi-noninvasive way [2]. It uses low coherence light interferometry in the near infrared to analyze the light reflected back or scattered back off tissue inhomogeneities. This technique provides high-resolution, real-time, three-dimensional images of the detailed morphology of transparent and semi-transparent biological tissues. More recently [3], as an important development of OCT, the first Doppler OCT images of intralipid solution flow in a channel located 1 mm beneath the surface of a scattering phantom have been obtained. Later that year, Izatt et al [4] reported Doppler OCT images of bidirectional flow in biological tissues where the velocity was color coded to indicate the magnitude and the direction of flow. These early Doppler OCT experiments employed short-time Fourier transform (STFT) to obtain the power spectrum from the temporal interference fringe intensity. The Doppler shift was then calculated with the centroid method. In our view, the STFT spectrogram never gives the «true value» of the instantaneous frequency and group delay even if a good approximation is sometimes obtained. There are a number of alternative implementations that can overcome these limitations. For instance, the phase-resolved Doppler OCT technique deduces the structural and velocity OCT images from the amplitude and phase information contained in the real and imaginary parts of the complex signal [5–7]. These papers and the recent phase-resolved techniques either in the time domain [8,9] or in the spectral domain [10,11], all rely on the Hilbert transform to extract the phase information from the OCT interferograms in one way or another. Correlation processing performs the correlation integral on the sequential in-depth scan signals and determines the complex signal by a phase sensitive demodulation of the photocurrent [9,12,13]. The primary objective sought by these last authors is to achieve real-time and especially *in-vivo* Doppler imaging thus putting the emphasis on hardware implementation. Our purpose is somewhat different. It can be stated as follows: given an interferogram obtained either via the time-domain or frequency domain OCT we seek the most precise time-frequency method for Doppler OCT. For our investigation, we use experimental data and perform the signal processing comparison numerically off-line putting aside hardware considerations or particular implementations in an attempt to provide as much as possible a balanced comparison. As a benchmark, we will pay special attention to STFT but other methods will be considered as well. It is our view that as a time-frequency analysis method, the Wigner distribution can provide more accurate instantaneous frequency and group delay values of the signal at a given time compared to STFT or the Hilbert transform method. With respect to STFT, a similar conclusion was reached in [14] although in a different context. In this paper, we present a Doppler-OCT flow-metrology system based on the Wigner distribution in a setup that yields cross-sectional flow-velocity maps. As an example, the approach is tested by providing the Doppler-OCT velocity distribution of a laminar flow in a rectilinear cylindrical tube using the scattering data obtained from an intralipid solution.

## 2. Principles

In addition to the usual structural information which is extracted from the position and the amplitude of the OCT cross-correlation trace, Doppler OCT also yields precious information on the frequency content of the detected fringe intensity coming from the interference between the reference and target arms. These interference fringes can be observed only when the optical path lengths between the light backscattered from the sample and that from the reference arm are matched within the coherence length of the broad bandwidth light source. The flow velocity measurement is based on the principle that the Doppler frequency shift Δf_{D} in the interference fringe intensity modulation results from the motion of the scatterers within a sample [5]:

where k⃗_{s} and k⃗_{i} are the wave vectors of the incoming and scattered light, respectively, and v⃗ is the velocity of the moving scatterer. Given the angle θ between the probe beam and the direction of motion of the scatterer, the change in frequency is related to the velocity by Eq. (1):

where v_{ODT} is the velocity at a specific point, λ_{0} is the center wavelength of the source, Δf_{D} is the Doppler frequency shift and n̅ is the local index of refraction in the sample.

The simplest and the most commonly used time-frequency analysis is the short-time Fourier transform (STFT). The basic idea of the short-time Fourier transform is to divide the signal into small time segments (also called STFT windows), and apply a Fourier transform to each time segment in order to extract the frequencies that are present in that time interval. The window is then shifted to the next time segment in order to generate another local spectrum. The totality of such spectra provides the power distribution in the signal with respect to time and frequency simultaneously:

where w(t+τ) is a short time analysis window, y(t) is the signal. In STFT, the spectral centroid method is frequently used for each window to obtain the first moment of the power spectrum f_{c} within each time segment.

The final flow velocity map is produced by displaying these first moments at each pixel. In contrast, the Wigner distribution of a signal y(t) is defined as:

Unlike the STFT which extracts the frequency information as a function of time by dividing the signal into many small time segments and Fourier analyzes each one of them, the Wigner distribution at a particular time t is determined by folding the left part of the signal over the right part and Fourier transforming the product. It is known that one of the advantages of the Wigner distribution over STFT is that one does not have to bother with the choice of a window [15]. In particular, the Wigner distribution is always real-valued. It preserves time and frequency shifts and it satisfies the following marginal distribution properties:

where Y(ω) is the frequency energy distribution of the signal. This means that if we integrate the time-frequency energy density along one variable, we obtain the energy density corresponding to the other variable. To improve the signal-to-noise ratio of the Wigner distribution, a function h(τ) peaked around τ = 0, has been introduced to define a pseudo Wigner distribution. This way, we make the Wigner distribution local for the purpose of emphasizing the signal around time t. The definition of the pseudo Wigner (PW) distribution is:

Note that the Hanning window h(τ) represents the same STFT window w(t+τ) as the one used in Eq. (3). In this work, we used the discrete versions of Eqs. (3) and (9) found in [16]. The discrete STFT version can also be found in [4] and the PW version in [17], Eq. (7.10). From now on, we mean the **pseudo Wigner method** when we refer to the Wigner method.

## 3. Experimental setup

To test our technique, we performed flow measurements with a Doppler OCT system depicted in Fig. 1. This experimental setup is based on a fiber-optic Michelson interferometer which uses a superluminescent diode (SLD) as the low-coherence light source. This source is characterized by a central wavelength of λ_{0} = 1545 nm and a bandwidth of Δλ = 32 nm. Light is focused onto the sample with a near-infrared, low-numerical aperture (0.25 N.A., 10X magnification) microscope objective in order to get a long depth of field. The sample and reference mirror scans are provided by motorized stages moving at a typical speed of 1.5 mm/s. The output optical power is measured with a dual balanced detector that cancels the random intensity fluctuations, thus improving signal-to-noise ratio. Polarization controllers are used in the fiber interferometer in order to get the maximum interference fringe visibility between the reference and sample arms. In the sample arm, a 1.4 mm inner-diameter glass capillary tube with a 1% intralipid in water solution flowing inside it was fixed on an adjustable stage. The flow rate was adjusted by a controller and a pump. The angle between the probing beam and the flow direction was 80°. The carrier frequency of the fringe signal was 1.94 kHz since the scan speed of the reference mirror was 1.5 mm/s. OCT images of the tube section were captured to calculate the velocity profile. The mean flow velocity was set at 16 mm/s. As the Reynolds number is much less than the critical number, the flow was always laminar and therefore, the velocity profile in the circular tube was parabolic [18]:

Here V(r) is the velocity at radial position r, V_{c} is the central peak velocity, and d is the inner diameter of the tube.

## 4. Results and discussion

The comparison of the data processing methods was based on a full 2D-OCT profile. The background of Fig. 2 shows the OCT image of a section of the capillary glass tube with intralipid flow described in the experimental setup. As an application of the Wigner distribution method, we calculated the maximum frequencies at each pixel shown in Fig. 2 in order to get a contour velocity map.

The measured flow profiles of each A-scan were fitted with a quadratic fit. Thus, a velocity distribution within the cross sectional image can be built. In Fig. 2, the contour frequency shift map is superimposed to the structural image in order to give a better appreciation of the correlation between the actual sample and the fluid velocity distribution at any position inside the tube.

#### 4.1 Comparison with STFT

Many limitations of time-frequency analysis arise from the compromise between spatial resolution and velocity sensitivity. In a typical approach, the accuracy of the flow velocity calculated with STFT depends on the window size of the Fourier transform for each pixel. Hence, the minimum detectable Doppler intensity has to be calculated from at least one period of the average frequency within the time segment of the window.

The full A-scan Doppler OCT interferometric signal that is indicated by the white line on Fig. 2 has been submitted to both STFT and Wigner distribution methods. By zooming in onto the signal, we can see some details of the signal processing. Fig. 3 shows a small time span which includes 94 points of the data sampled at a rate of 100 kHz leading to a 10 μs time interval between two points. A Hanning window function is used for apodization and the size of the window is chosen so as to contain 17 points which leads to a 0.17 ms time interval. Such an interval guarantees that each window contains at least one frequency cycle.

In Fig. 3, we show data from a single A-scan where the amplitude has been normalized by dividing the interferogram by the local envelope. We have verified that the normalization procedure has a negligible effect on the position of the zero crossing points which determine the frequency distribution.

In STFT, one usually performs the Fourier transform for each window. Although one could display the time-frequency distribution for each window, we show in Fig. 4 a «sliding-window» plot using the above Hanning window characteristics. The same characteristics will be used in the pseudo-Wigner approach.

In comparison, the Wigner distribution is capable of tracing locally this instantaneous frequency variation at each sampling point with a better resolution. The boxed region in Fig. 5 clearly shows the frequency variation with time. An advantage of the Wigner distribution over STFT is that the calculated bandwidth of the spectrum is more precise. Fig. 6 displays two frequency spectra calculated with the Wigner method at two specified time values (points No. 52 and No. 59 from the signal in Fig. 3) and compares them with the frequency distributions calculated by STFT over a window delimited by points 52 and 68.

The spectral bandwidths of the frequency plots using the Wigner distribution approach are significantly narrower than those calculated with STFT. Since one period contains between 10 and 12 sampling points, the frequency of the signal should be between 8.3 and 10 kHz. From the Wigner time-frequency distribution, the calculated frequencies at point 52 and 59 are 8767.1 and 9799.5 Hz respectively if we use the centroid method, while their maximum frequencies are 8708 and 9491.2 Hz, respectively. It seems to make sense that the maximum frequency is close to the frequency calculated by the centroid method for this particular Wigner distribution spectrum. However, while the maximum frequency from the STFT spectrum (9491.2 Hz) lies within the estimated frequency range, the frequency calculated from the centroid method is overestimated (11683 Hz).

Part of the broad bandwidth spectrum corresponds to negative frequency components. This part is not accounted for in the centroid calculation. Moreover, side lobes appear for higher frequency components. This arises as a consequence of the well known spectral leakage in STFT caused by a relatively short window function. While a more extended analysis could provide additional information, one should also account for spurious artifacts caused by distortions in the shape of the signal due to the low sampling rate and to speed fluctuations in the feedback-controlled motorized stage that drives the in-depth scan. Thus, unwanted frequency components appear in the recorded signal. This indicates that the calculated maximum frequency is more likely to give a better estimate of the real frequency than the frequency calculated using centroid method. In Fig. 7 (b), the 55 spectra calculated with the Wigner method at each sampling point shown in Fig. 7 (a) have been plotted. This figure shows that the main problem with the Wigner method might lie as well with these side lobes appearing in the calculated spectra. These side lobes result in a systematic overestimation of the frequency calculated with the centroid method. Moreover, if we retain only the peak frequencies, the side lobes do not even contribute to the final results, leading to a better estimate of the frequencies in the signal. These spectra clearly show that the maximum frequency peaks all lie within a small range and give a reasonable approximation of the actual frequency. Therefore, as just shown, the maximum frequency method is much more adequate than the centroid method and will be used exclusively.

#### 4.2 Other methods

For a more extended performance comparison, various methods were used for processing the data of the A-scan represented by the white line on Fig. 2. First a word concerning the time-domain and spectral-domain OCT methods [10–11]. Both lead eventually to the interferograms which are the starting point for our digital signal processing comparison. The spectral-domain OCT is known to give a better SNR and dynamic range and should be advantageous for weak signals. In our case, this difference is alleviated by the fact that we use a normalized interferograms for the phase determination. With respect to the phase-resolved methods, in most cases found in the literature [6,7], the Hilbert transform is used as a technique for obtaining the phase from the corresponding analytic signal. We will therefore refer to this class of phase-resolved methods as the Hilbert transform method. Another approach also found in the literature [8,13], is the autocorrelation method borrowed from ultrasound Doppler techniques. This method appears to have certain advantages in terms of speed for real-time hardware implementations. Nevertheless, if the velocity range is too large, this method introduces ring features that must be dealt with. Also the delay inherent to this method must be chosen appropriately.

As mentioned previously, we have chosen to compare the various methods on a similar platform singling out the signal processing aspects to the exclusion of other features in order to discover how they compare relatively to the theoretical result, at least in the present capillary tube experiment. Again, in a spirit of increased commonality, we have used the Matlab Time-Frequency Toolbox [16] in all cases, keeping in mind that the comparison bears not so much on the speed of execution of this or that algorithm but on the precision that is eventually achieved. Still, a few comments relative to the latter point will be made later on.

#### 4.3 Methodology and results

Our analysis starts with the normalized interferograms data. The primary objective is to obtain the time-frequency two-dimensional (2-D) scatter plots which we fit to a parabola after noise filtering. In the case of STFT and Wigner, it is done through the three-dimensional frequency distribution plots and only the maximum value of the distribution is retained in order to generate the 2-D plots. In the case of the Hilbert phase-resolved method, the phase and frequency are obtained from the quadrature components whereas for the autocorrelation method, the flow diagram of [13] is reproduced with a Matlab simulation. What is believed to be the best delay was chosen in that particular case. The 2-D time-frequency scatter plots are then processed by the same noise filtering procedure. They are finally fitted to a two parameter parabola where the velocity on the wall of the capillary tube is left floating. The results are plotted on Fig. 8. The uncertainty shown on the plots comes from the error on the determination of the position of the interface between fluid and the tube wall. The curve pertaining to the autocorrelation method was obtained for a delay corresponding to 0.05 msec for data sampled at a rate of 100 kHz which gave us the best result for that method. Going to higher delays introduced the ring structure and worsened the results in spite of phase unwrapping.

The flow velocity of Fig. 8 is deduced from frequency via Eq. (2). The maximum values of the curves give the central peak velocity V_{c}. This parameter was calculated for each method and was then reported in Table 1. In Fig. 8, a theoretical flow profile is also shown. This curve was obtained by repeating several measurements of the mean flow inside the tube with the help of a flowmeter. Assuming a laminar flow model, a theoretical value of V_{c} was then deduced by multiplying the mean velocity by a factor 2. The theoretical flow profile is then plotted by inserting V_{c} into Eq. (10).

Following our investigation, it appears as if the Wigner distribution method is the one that offers the best precision overall. Near the maximum, the Wigner and STFT methods are close but for low velocities, the STFT method overestimates the theoretical result substantially. The phase-resolved Hilbert method gives good results for low velocities but appears to be slightly inferior at higher velocities. Finally the autocorrelation method is definitely not so good even though we phase-unwrapped a secondary ring in our analysis. Apart from the Wigner transform method which is not considered in that paper, our results are consistent with [22].

#### 4.4 Discussion.

In spite of a great deal of research, the Wigner distribution method remains somewhat unfamiliar. We can list a few important features taken from the abundant literature, for instance [15,23]. The Wigner-Ville transform has several advantageous properties for the time-frequency analysis of nonstationary signals. It has a bilinear form and is compatible with linear filtering and modulation. One can get the instantaneous frequency, the spectral density, the instantaneous power and group delay by computing moments. It provides a real-valued representation. The Wigner function gives the clearest picture of the instantaneous frequency and group delay. The pseudo Wigner distribution is a smoothed version of it. Nevertheless while the conditional averages are exactly the instantaneous frequency and group delay, the pseudo Wigner distribution loses this characteristic. This distribution is known to suffer from artifacts in form of cross-terms especially when the frequencies change abruptly. For instance, the observed artifacts of a broadening of the spectra at both ends of Fig. 5 result from the onset and ending of the computation which can be assimilated to an abrupt frequency change. It is often perceived that the discrete Wigner distribution as a time-frequency spectral estimator is computationally intensive. In fact, the Wigner distribution can be done with an FFT network half the size of the one needed to perform a standard STFT [19]. In [17], it is stated that the pseudo Wigner distribution can be computed with a complexity of ½M(log_{2}M+4) where M is the length of the window. This is of the same order of magnitude as STFT. Finally in another paper [20], it was shown that the Wigner transform is 37.5% faster than the Fourier transform while having the similar accuracy and noise sensitivity for signals with high SNR. The Hilbert transform requires 2nlog(n)+n complexity operations where n is the number of samples per A-line as stated in [21]. For fast real-time applications, the Hilbert method with hardware in-phase and quadrature demodulation and Kasai autocorrelation over several A-scans [21] appears to offer the best tradeoff at the present time with a slight disadvantage on precision. A recently introduced method, based on an adaptive notch filter appears to offer fast processing speed and «fine spatial resolution» [24]. It seems though that the step size is a critical parameter which may lead to loss of precision if it is too small or non convergence, if too large. We did not consider this method in the present paper.

## 5. Conclusion

In conclusion, the performance of various signal processing methods in Doppler OCT was compared for the determination of the Doppler frequency shifts or velocity in a flowing intralipid solution. We found, based on experimental data, that among the techniques we investigated, the Wigner method is more precise, especially in the higher velocity regime. Whereas the STFT method is close in the same range, it is substantially worse for low velocities. The next best method appears to be the phase-resolved Hilbert distribution method and, based on our analysis, the worst of all turns out to be the autocorrelation method although it is often used for real-time imaging. To the best of our knowledge, this is the first time that the Wigner time-frequency distribution method is applied to the data processing of Doppler OCT. The results show that the Wigner distribution method offers a better estimate of the flow velocity even though the difference is not overwhelming. It should therefore be preferred for precision work.

## Acknowledgements

The authors acknowledge the financial support of the Canadian Institute for Photonic Innovation and NSERC. Moreover, they thank the anonymous reviewers for their judicious comments.

## References and links

**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–1181 (1991). [CrossRef] [PubMed]

**2. **W. Drexler, U. Morgner, F. X. Kartner, C. Pitris, S. A. Boppart, X. D. Li, E. P. Ippen, and J. G. Fujimoto, “In vivo ultrahigh-resolution optical coherence tomography,” Opt. Lett. **24**, 1221–1223 (1999). [CrossRef]

**3. **Z. P. Chen, T. E. Milner, D. Dave, and J. S. Nelson, “Optical Doppler tomographic imaging of fluid flow velocity in highly scattering media,” Opt. Lett. **22**, 64–66 (1997). [CrossRef] [PubMed]

**4. **J. A. Izatt, M. D. Kulkami, S. Yazdanfar, J. K. Barton, and A. J. Welch, “In vivo bidirectional color Doppler flow imaging of picoliter blood volumes using optical coherence tomograghy,” Opt. Lett. **22**, 1439–1441 (1997). [CrossRef]

**5. **Z. P. Chen, Y. H. Zhao, S. M. Srinivas, J. S. Nelson, N. Prakash, and R. D. Frostig, “Optical Doppler Tomography,” IEEE J. Sel. Top. Quantum Electron. **5**, 1134–1142 (1999). [CrossRef]

**6. **Yonghua Zhao, Zhongping Chen, C. Saxer, Shaohua Xiang, J. F. de Boer, and J. S. Nelson “Phase-resolved optical coherence tomography and optical Doppler tomography for imaging blood flow in human skin with fast scanning speed and high velocity sensitivity,” Opt. Lett. **25**, 114–116 (2000). [CrossRef]

**7. **L. Wang, W. Xu, M. Bachman, G. P. Li, and Z. P. Chen, “Phase-resolved optical Doppler tomography for imaging flow dynamics in microfluidic channels,” Appl. Phys. Lett. **85**, 1855–1857 (2004). [CrossRef]

**8. **D. Morofke, M. C. Kolios, I. A. Vitkin, and V. X. D. Yang “Wide dynamic range detection of bidirectional flow in Doppler optical coherence tomography using a two-dimensional Kasai estimator”, Opt. Lett. **32**, 253–255 (2007). [CrossRef] [PubMed]

**9. **Andrew M. Rollins, Siavash Yazdanfar, Jennifer K. Barton, and Joseph A. Izatt, “Real-time in vivo color Doppler optical coherence tomography”, J. Biomed. Opt. **7**, 123–129 (2002). [CrossRef] [PubMed]

**10. **B. J. Vakoc, S. H. Yun, J. F. de Boer, G. J. Tearney, and B. E. Bouma, “Phase-resolved optical frequency domain imaging”, Opt. Express **13**, 5483–5493 (2005). [CrossRef] [PubMed]

**11. **Brian R. White, Mark C. Pierce, Nader Nassif, Barry Cense, B. Hyle Park, Guillermo J. Tearney, Brett E. Bouma, Teresa C. Chen, and Johannes F. de Boer, “In vivo dynamic human retinal blood flow imaging using ultra-high-speed spectral domain optical Doppler tomography”, Opt. Express **11**, 3490–3497 (2003). [CrossRef] [PubMed]

**12. **S. Yazdanfar, A. M. Rollins, and J. A. Izatt “Ultrahigh velocity resolution imaging of the microcirculation in vivo using color Doppler optical coherence tomography”. Proc. SPIE , **4251**, 156–164 (2001). [CrossRef]

**13. **Andrew M. Rollins, Siavash Yazdanfar, Rujchai Ung-arunyawee, and Joseph A. Izatt “Real time color Doppler optical coherence tomography using an autocorrelation technique”, Proc. SPIE. **3598**, 168–176 (1999). [CrossRef]

**14. **C. Xu, F. Kamalabadi, and S. Boppart, “Comparative performance analysis of time-frequency distributions for spectroscopic optical coherence tomography,” Appl. Opt. **44**, 1813–1822 (2005). [CrossRef] [PubMed]

**15. **Leon Cohen, *Time-frequency analysis*, Prentice Hall, Englewood Cliffs, N.J. (1995).

**16. **F. Auger, P. Flandrin, P. Gonçalvès, and O. Lemoine, “Time-Frequency Toolbox tutorial”, CNRS (France), Rice U. (U.S.A.), http://tftb.nongnu.org/ (2005) and http://gdr-isis.org/tftb/tutorial/tutorial.html.

**17. **T. C. M. Claasen and W. F. G. Mecklenbräuker “The Wigner distribution - a tool for time-frequency signal analysis” Part II: Discrete time signals. Philips J. Res. **35**, 276–300 (1980).

**18. **Wu L., “Simultaneous measurement of flow velocity and Doppler angle by the use of Doppler optical coherence tomography,” Opt. Lasers Eng. **42**, 303–313 (2004). [CrossRef]

**19. **N. Bergmann “New formulation of discrete Wigner-Ville distribution”, Electron. Lett. **27**, 111–112 (1991). [CrossRef]

**20. **T. P. Zielinski, “Wigner transform instantaneous phase estimator”, Eusipco-96, Trieste, PDE.10 (1996)

**21. **Victor X. D. Yang, M. L Gordon, A. Mok, Y. Zhao, Z. Chen, R. S. C. Cobbold, Brian C. Wilson, and I. Alex. Vitkin “Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation”, Opt. Comm. **208**, 209–214 (2002). [CrossRef]

**22. **Daqing Piao, Linda L. Otis, Niloy K. Dutta, and Quing Zhu, “Quantitative assessment of flow velocity-estimation algorithms for optical Doppler tomography imaging”, App. Opt. **41**, 6118–6127 (2002). [CrossRef]

**23. **T. P. Zielinski, “On a software implementation of the Wigner-Ville transform”, Comp. Phys. Comm. **50**, 269–272 (1988). [CrossRef]

**24. **Y. Chen, P. Willett, and Q. Zhu, “Frequency tracking in optical Doppler tomography using an adaptive notch filter”, J. Biomed. Opt. **12**, 014018-1 to 014018-9 (2007). [CrossRef]