1. Introduction

Optical Coherence Tomography (OCT) is an imaging technique that appeared about 15 years ago [1]. A few years later, in 1997 [2], as an important new development, the first spatially resolved Doppler OCT images of liquid flow in a phantom have been obtained. Since then, a number of experiments both ex-vivo and in-vivo have taken place. This technique can detect the flow velocity by measuring the flow-induced Doppler shifts during each A-scan. A great deal of effort has been spent on searching for accurate and efficient methods for processing the experimental data. In addition to the early Doppler OCT experiments which used short-time Fourier transform (STFT), there already exist a number of alternative implementations that have been developed for that purpose. For instance, the phase-resolved Doppler OCT technique deduces the structural image and velocity map of a sample from the amplitude and phase information contained in the real and imaginary parts of the complex signal [3–5]. In the time-domain OCT, the recent phase-resolved techniques rely on the Hilbert transform to extract the phase information from the OCT interferograms [4,6]. A hardware method based on autocorrelation of the OCT signal with variable delay time was proposed [7]. Doppler implementations have also been successfully tried in Spectral Domain (SD) [8].

For structural OCT, the main advantages of SDOCT compared to TDOCT lie in a faster image acquisition rate and a higher signal-to-noise ratio (SNR) [9]. In terms of global superiority, for Doppler OCT, we would like to argue that the situation is not as clear as generally assumed because other considerations such as the achievable Doppler-velocity span (dynamic range) and measurement accuracy must be taken into account. Usually, high accuracy is sacrificed for a high scan rate and justly so in many situations. On the other hand, for applications where precise Doppler-speed metrology is essential, the high scan rate imperative can be relaxed. From a survey of the literature, we found that the achieved speed sensitivity for time-domain systems was of the order of 7 µm/s [6] whereas for spectral systems, it was more in the range of 50 µm/s [10] but no doubt, the situation evolves constantly. Similarly, for the maximum measurable speed, we found 560 mm/s for TDOCT [6] and 191 mm/s for SDOCT [11]. In terms of A-scan rates, Doppler SDOCT systems achieve 25 kHz [12] and Doppler TDOCT systems, 8 kHz [6]. For a number of applications such as microfluidics and ex-vivo experiments, the scan rate is less important than accuracy and the 8 to 10 kHz range can be quite comfortable, way above the values used in recent publications which are in the sub 10 Hz range [13–16]. The previous figures, while representative of the current situation, do not set the limits of the best achievable performance, of course. For instance, we estimate that the time-domain zero-crossing method is much less limited in speed sensitivity (at low speed) than what has been quoted so far because as the frequency gets lower and lower, the number of samples per period at fixed sampling rate gets higher and hence gives a better sensitivity: 1 µm/s is certainly achievable. At the high frequency end, the limit is set by the sampling hardware which nowadays is easily in the range of 100 MHz at 14 bit resolution. One can envisage a Doppler velocity of a few tens of meters per second. In summary, we observe that the Doppler velocity dynamic range in present day SDOCT systems is given to be in the vicinity of 45 to 50 dB but we find that in TDOCT systems, it can achieve 80 dB or more especially at slow scan rates.

Another important issue, actually the main concern for this paper, is the accuracy and ease in Doppler velocity determination. In SDOCT, the only available method to calculate the Doppler velocity is to use phase resolved methods from the inverse Fourier Transform of the spectral signal [11,12]. Taking the imaginary part of the temporal signal gives the phase information for each A–scan. The flow speed is determined by calculating the phase difference between two successive A-scans [10]. This method suffers from an aliasing phenomenon caused by 2π ambiguities in the relation between the flow and the phase that limits the maximum determinable Doppler shift. This maximum velocity corresponds to a phase change of 2π within a time between two successive A-scans. For typical SDOCT systems, the axial scanning rate is around 20 kHz, which brings typical maximum aliased speed values of several tens of mm.s^-1. Velocities higher than this range need to be obtained with phase unwrapping methods [11], but their use appears to be limited either to simple or laminar flow profiles. SDOCT based on CCD detection is also limited by another issue which is fringe washout during the exposure time for A-scan acquisition [18]. In addition to possibly more accurate values and wider dynamic range in the case of time domain methods, these phase ambiguities and indirect phase extraction methods are strong motivations for our decision to consider time-domain OCT as the preferred approach for accurate Doppler shift measurement.

In addition, we feel that the claim of better SNR performance for spectral-domain systems as opposed to time-domain systems needs to be revisited in the case of Doppler OCT because in the latter case, the signal ought to be the Doppler frequency shift and the noise, ought to be the FM noise, and not so much the AM (amplitude modulation) SNR considered previously [9,19], although both are related since the phase noise variance is inversely proportional to the square root of the SNR. Also, some special features associated with our zero-crossing method proposed in this paper directly affect the SNR as discussed in Section 4.4.

Obviously, all the proposed approaches aim at obtaining rapidly, precisely and as much as possible effortlessly, the instantaneous frequency of the time-varying signal present in the Doppler data. As we have shown in a previous paper [20], the Wigner method turned out to be the most precise data processing method in Doppler TDOCT so far, at least in the laminar flow example we have studied. We expect this will also apply to non laminar flow. The global deviation from the theoretical value for laminar flow, averaged over one A-scan, is respectively 22.08 % for the autocorrelation method, 10.32 % for STFT, 2.97 % for the Hilbert method and 1.56 % for the pseudo-Wigner method. In this paper, we show that, as far as precision is concerned, the zero-crossing (ZC) method outperforms the Wigner method in general, and hence all other existing signal processing methods in time-domain Doppler OCT. The main reason is that, because in Doppler OCT we obtain the spatially (or temporally) resolved frequency shifts from strongly localized scatterers (to first order), much of the drawbacks observed in laser velocimetry applications of the zero-crossing method are strongly mitigated. Zero-crossing methods and algorithms have been applied in various fields such as signal and image processing [21], optics [22], biomedical engineering [23] and fluid mechanics [24]. Indeed, zero-crossing detectors were often used for Doppler flow velocity measurements in both laser Doppler velocimetry [25] and ultrasound Doppler velocimetry [26].

In the case of ultrasound Doppler velocimetry, the hardware detectors provide a sharp pulse that matches very closely the zero-crossing voltage of the primary signal. The frequency response is then obtained by counting the number of pulses over a certain period of time (typically several tens of ms). Initially, the wide dynamic range and relative simplicity of hardware implementation made these detectors very attractive for applications such as blood flow measurements [27]. This technique can be shown to be adequate for simple periodic signals. Nonetheless, in these applications, the received signal is the superposition of composite Doppler signals induced by the presence of many scatterers with different velocities present in an extended volume of the sample. Under the assumption that there is no phase relation between the different scatterers, it has been shown that the output frequency from the zero-crossing detector is the root mean square frequency of the signal [28]. Thus, in the case of wide and complex frequency spectra (for instance in the regions of turbulent flow), the zero-crossing detector only gives a rough approximation of the mean frequency [29, 30]. Furthermore the various noise sources cannot be separated out and may strongly corrupt the results. For these reasons, software methods based on the computation of the complete power spectrum gradually replaced the ZC detectors and other more precise approaches such as the Fast Fourier Transform [31] were preferred.

In the case of laser Doppler velocimetry, software and hardware implementations of zero-crossing counting were also developed. Like the ultrasound Doppler velocimetry counterpart, the signal results from a superposition of many Doppler signals originating from scatterers of different velocities found in an extended sample volume. In this case, the zero-crossing method is based on the measurement of the time elapsed between each successive zero-crossing point. A number of researchers have applied a method which is based on the statistical properties of the zero crossing intervals to obtain the complete Doppler frequency spectrum [32]. It appears that this method works only for low SNRs and low depths of field [33].

The aim of this paper is to propose a simple software zero-crossing detection method for time-domain OCT Doppler implementation. Its performance in terms of accuracy is then compared to other techniques. As a benchmark, we will pay special attention to the pseudo-Wigner method since this method was found to be the most accurate so far [20]. Nevertheless other methods (such as the Hilbert transform method) will be considered as well because of their widespread use. For our investigation, we take 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 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.

As we show in this paper, the zero-crossing algorithm applied to the data processing of the Doppler OCT signal provides the best estimate so far of the Doppler frequency with an acceptable spatial resolution in the ≪textbook≫ case of the laminar flow in a rectilinear tube which we used as a phantom. This is due to the characteristics of the OCT signal since it is time resolved. In OCT, an interference between the reference beam and the backscattered sample beam occurs only when their optical path difference are matched within the coherence length of the source, which is usually on the micrometer scale (temporal gating). Thus, the frequency analysis of an OCT signal within a small time interval relates to the Doppler flow of the scatterers that are located within the corresponding minute volume. In the case of laminar flows, the local OCT interference signal is almost a single periodic waveform since the local Doppler frequency spectrum is so narrow. We believe that the relative simplicity of the spectral content of such signals is beneficial to the use of the zero-crossing approach in Doppler OCT. The situation for non-laminar flow remains to be investigated. In our view, one of the main reasons the zero-crossing method gives better results is that it is essentially local giving the instantaneous frequency accurately at a given time whereas the other methods are non-local since they involve integrals or mathematical transforms (Fourier, Wigner, and Hilbert) which need to probe the data over an extended range.

2. Principles

2.1 Principles of Doppler OCT

The principles of Doppler OCT have been well explained in the literature [2–6]. We provide here a brief description. The flow velocity measurement is based on the principle that the Doppler frequency shift Δf_D in the interference fringe intensity modulation of an OCT A-scan results from the motion of the scatterers within the 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 the following equation:

where V_ODT is the velocity at a specific point, λ₀ 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.

2.2 The Wigner method

As demonstrated in our previous paper [20], the Wigner method is up until now the most precise data processing approach for velocity tracking, especially for high velocities. Some details about the Wigner and Pseudo-Wigner distributions are given in the Appendix A. From now on, we mean the Pseudo-Wigner method when we refer to the Wigner method.

The Wigner method provides the frequency spectrum of the signal at each sampling point. The broad bandwidth of each spectrum comes from many factors, the most important of them being the time-frequency trade-off introduced by the windowing operation of the Wigner distribution. We use a Hanning window when passing from the Wigner to the Pseudo-Wigner distribution as explained in [20]. Other spectral broadening factors come from 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. Based on our analysis [20], a better estimate of the frequencies in the signal can be achieved if we retain only the peak frequencies of the spectra at each sampling point. This way, the previously obtained 3D time-frequency distribution becomes a time-frequency or position-velocity 2D plot which maps the velocity distribution within the sample.

2.3 The zero-crossing method

As mentioned before, the zero-crossing method works well for Doppler OCT in our case because the local OCT interferogram is almost a singly periodic waveform. The 2D time-frequency plot can therefore be obtained directly by calculating the half-period Δt=(t_z)_N+1-(t_z)_N between two consecutive zero-crossing points (t_z)_N and (t_z)_N+1 as shown in Fig. 1. We assume that within two such points, the instantaneous frequency of the signal F can be taken to be constant and equal to the inverse of twice the half-period Δt:

The zero-crossing points t_z are calculated by finding all the pairs of adjacent signal values S₁ and S₂ that have opposite signs. If the sampling instants are t₁ and t₂, then the zero-crossing time t_z is determined by a simple linear triangulation formula as illustrated in Fig. 1:

where τ is the sampling rate.

In contrast to the Wigner method, the time resolution of which is determined by the sampling rate only, this zero-crossing method has two factors that may reduce its time-frequency resolution, as inferred from Fig. 1. The first one is the uncertainty δF that arises from the sampling of the signal. Because of the phase and frequency variations inherent to the OCT signal, the actual zero-crossing point may be located anywhere between the points S₁ and S₂, even though there is a strong probability that it lies near the point determined by Eq. (4). The error on the half-period is:

The relative error on the frequency is then:

The other factor that reduces the time-frequency resolution is the spatial resolution δx in the axial direction of the A-scan. As the frequency can only be calculated at zero-crossing points of the signal along the axial direction, δx is determined by the distance between two consecutive zero-crossing points:

where V_scan is the speed of the axial scan, and n is the index of refraction of the sample. In order to avoid confusion with the spatial resolution of the structural OCT image, δx is called the effective sampling parameter of the frequency. From Eq. (3), we then obtain:

 

Fig. 1. A sample signal used to calculate the zero-crossing points, red cross, and the defined time durations.

Download Full Size | PDF

Thus, for a fixed sampling rate (i.e. for fixed τ) and fixed axial scan speed, Eqs (5) and (7) reveal a tradeoff between effective sampling and relative error on frequency that depends only on the flow velocity. In other words, when the frequency of the signal under measurement is very low (i.e. for a long wavelength optical source and low flow velocities), the large separation between two consecutive zero-crossing points leads to a less precise localization on the spatial velocity map. On the other hand, there will be more sampling points per half-period leading to a better velocity precision expressed by δF/F. In the case of high frequency signals (i.e. for a short wavelength optical source and high flow velocities), the velocity map will display a higher effective sampling of the velocity while the velocity precision will be less due to a lower number of sampling points per half-period.

To illustrate this point, the effective sampling parameter and the frequency precision were calculated in the case of an ad-hoc A-scan signal containing different Doppler frequency components from 2 to 10 kHz obtained at a fixed sampling rate of 100 kHz and for an axial scan speed of V_scan=1.5 mm.s^-1. The results are shown in Fig. 2. We see that a combination of limited frequency precision and limited effective sampling is obtained for intermediate frequencies (around 5 kHz). In that case, a typical effective sampling parameter of 0.2 µm is combined with a frequency precision of 20 %. At first sight, a further increase of the measured Doppler frequency would decrease its relative precision (higher) δF/F. Nevertheless, when this precision becomes too poor, (i.e. the number of data points per half-period is too low) the measurement of the half-period is too coarse. In that case, the frequency measurement is done by counting several zero-crossing points over a time range that is large enough to collect a sufficient number of data points. As a consequence, the frequency precision becomes better at the expense of a poorer effective sampling. Anyway, in our experiment, we are far from this case because the data sampling rate is one order of magnitude higher than the maximal average measured frequency.

 

Fig. 2. Effective sampling parameter of map velocity and precision of calculated frequency as a function of Doppler frequency F.

Download Full Size | PDF

In a nutshell, the zero-crossing method processes data by dividing the signal into variable length time windows, each of which is bordered by consecutive zero-crossing points. The duration of each time window is determined and taken to be one half the local period. The frequency of this time window is then calculated as the inverse of the period. After that, the window is shifted to the next time segment in order to generate another local frequency. Apart from accuracy considerations discussed previously, this method leads to a computational time saving of 26 compared to the Wigner method. There are two reasons for that. First, no integral or convolution needs to be calculated and the arithmetic operations involved within each time step remain simple. Also, the calculations are performed only on a subset of all the data points, namely at the positions where the signal amplitudes are equal to zero, but this feature leads to a lower-density spatial grid for the velocity map, as explained previously.

3. Experimental setup

The proposed zero-crossing method has been tested with the experimental conditions described in a previous paper [20]. The source used is a superluminescent diode characterized by a central wavelength of λ₀=1545 nm and a bandwidth of Δλ=32 nm. In that experiment, the sample is a 1.4 mm inner-diameter glass capillary tube with a 1% intralipid in water solution flowing inside. 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 obtain 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 [16, 34]:

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

4.1 Comparison with Wigner method

Let us select a single A-scan taken along the center of the tube to illustrate the procedure. Figure 3(a) shows the original signal at a sampling rate of 100 kHz. We first checked that the signal was symmetrical with reference to the zero amplitude axis. This important condition is ensured by the balanced detectors used for measuring the output signal from the OCT system. By dividing the interferogram by its local envelope, we obtained the normalized signal shown in Fig. 3(b). We have verified that the normalization procedure has a negligible effect on the position of the zero crossing points which determine the frequency distribution.

 

Fig. 3. (a). The original interferogram signal along the center of the tube and (b) the normalized data.

Download Full Size | PDF

First, we submit this signal to the Wigner method and a color-coded 3D time-frequency plot, as shown in Fig. 4, is produced as a first step. The frequency axis is limited by the Nyquist frequency which is half the sampling rate of 100 kHz. A local, broad-bandwidth spectrum containing all the frequency components of the signal at each sampling point is obtained. In our previous paper [20], we have concluded that the peak frequencies of the spectrum at each sampling point give the best estimate of the local frequency of the signal at that moment. We therefore obtain a 2D time-frequency plot by choosing the maximum frequency of the frequency distribution at each temporal sample point, as illustrated in Fig. 5 (black dots). Since the received Doppler frequency information is quite noisy, further filtering is necessary. We use a statistical maximum-likelihood method to process the data. The essence of the procedure is to generate a histogram of the frequencies found within a given (short) time span. The data outside a FWHM (full width at half maximum) of the histogram for that time span are considered to be noise and are discarded since we assume that at least for a laminar flow across the tube, the frequency variation should be smooth and only dependent on the speed [35]. Hence, the complete A-scan data are divided into several time segments, each of which contains 500 sampling points. Within one time segment, the histogram was built in 500 Hz steps ranging from 0 to 50 kHz. The yellow trace in Fig. 5 corresponds to the data remaining after the filtering procedure.

 

Fig. 4. Color-coded 3D time-frequency distribution of the signal along the center of the tube using the Wigner method.

Download Full Size | PDF

 

Fig. 5. Two-dimensional frequency distribution of the signal using the Wigner method. Black dots indicate the maximum frequency values of the spectra at each measured point while the yellow trace corresponds to the filtered frequency distribution.

Download Full Size | PDF

The next step implements the zero-crossing method. The approach described in Paragraph 2.3 is applied to the A-scan. Fig. 6 provides an example of zero-crossing points (indicated by stars) calculated over a small part of the normalized signal shown in Fig. 3. We zoom in on 94 points of data sampled at a rate of 100 kHz leading to a 10 µs time interval between two points. We now assess the performance of the zero-crossing method by comparing its results with those provided by the Wigner method. In Fig. 7, we display the calculated instantaneous frequency of the signal shown in Fig. 6 using both the Wigner and the zero-crossing methods. The blue points indicate the instantaneous peak frequencies of the local spectra calculated using the Wigner method. The stars in Fig. 7 correspond to the local frequencies obtained by the zero-crossing method. Each resulting local frequency is positioned in the middle of the time interval between two adjacent zero-crossing points.

 

Fig. 6. Small segment of the detected signal at a 100 kHz sampling rate. The stars indicate the zero-crossing points.

Download Full Size | PDF

 

Fig. 7. Comparison of Wigner and zero-crossing time-frequency distributions of the sample signal of Fig. 6. The points correspond to the Wigner method using the maximum values of the local spectra. The stars come from the zero-crossing method.

Download Full Size | PDF

As we can see in Fig. 7, the results given by the Wigner and the zero-crossing methods match rather well except at the ends of the plot because the Wigner distribution is known to suffer from artifacts brought by cross-terms especially when the frequencies change abruptly. Since we show truncated data on this plot, the observed discrepancy at both ends of the plot results from the onset and ending of the computation which can be assimilated to an abrupt frequency change. The Wigner method has the distinct advantage of a higher effective sampling locally. In contrast, the zero-crossing method produces only 16 data points and does not need a Fourier transform at each data point. As said before, the processing time is extremely short. Both methods are capable of tracing locally the instantaneous frequency variation with an acceptable spatial resolution when the full data set is considered. The zero-crossing method is sparser but follows the instantaneous frequency of a signal arguably as well as the Wigner method and is potentially better for abrupt velocity changes. Note that the Wigner method also requires the intermediate processing steps of getting the local frequency distribution and finding the local maximum of that distribution.

Now, we apply the zero-crossing method to the full normalized A-scan data set. The 2D frequency distribution of the signal can be obtained directly by this method as shown in Fig. 8 (black dots). The same noise rejection procedure as the one described previously for the Wigner method is applied next to the signal. The yellow trace in Fig. 8 shows the sampling points remaining after the procedure.

 

Fig. 8. Two-dimensional frequency distribution of the signal by the zero-crossing method. The black dots indicate the obtained frequencies and the yellow trace shows the resulting frequency distribution after noise rejection.

Download Full Size | PDF

4.2 Other methods

For a more extended performance comparison, another widespread method was also used for processing the same data of the A-scan. As mentioned in our previous paper [20], the Wigner method was found to be the most precise up until then, and the next best method appeared to be the phase-resolved Hilbert distribution method. We have therefore chosen to compare these two processing methods with the zero-crossing method using a common platform based only on the signal processing aspects to the exclusion of other features in order to discover how well they compared with the theoretical result, at least in the present capillary tube experiment.

4.3 Methodology and results

Our analysis starts with the normalized interferogram data presented in Fig. 3. The time-frequency two-dimensional (2-D) scatter plots are obtained for the Wigner and zero-crossing methods and are then noise-filtered using the procedure described previously. In the case of the Hilbert phase-resolved method, the phase and frequency are obtained from the quadrature components. The resulting data points are finally fitted to a parabola based on Eq. (9) using the two parameters V_c and d, which means that the velocity at the wall of the capillary tube is left floating. The results are plotted in Fig. 9. The flow velocity information is finally deduced from the frequency shift via Eq. (2). The uncertainty shown on the plots comes from the error on the determination of the position of the interface between the fluid and the tube wall.

 

Fig. 9. Computed flow profiles using a quadratic fit.

Download Full Size | PDF

The maximum values of the curves give the central peak velocity V_c. In Fig. 9, 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. (9).

 

Fig. 10. Comparison of the deviations of the three methods from the theoretical results.

Download Full Size | PDF

From our investigation, these three methods give results in the same range and match the theoretical curve more or less accurately. At the zoomed scale on Fig. 9, the zero-crossing method appears to fit the theoretical curve almost perfectly. A more detailed comparison is necessary and has been illustrated in Fig. 10. The deviations of these methods from the theoretical result are shown. It appears that the zero-crossing method is the one that offers the best precision overall except in a very limited region close to the tube wall. The Wigner method appears to be the next best one. It is less precise in the higher velocity region though, but it is slightly better for low velocities. Finally, the phase-resolved Hilbert method appears to be less accurate.

The comparison of the data processing methods was then extended to a full 2D-OCT profile. The gray-scale background of Fig. 11 shows the OCT image of a section of the capillary glass tube with intralipid flow described in the experimental setup.

 

Fig. 11. OCT image of the tube profile. The superimposed contour plots correspond to the frequencies calculated using experimental values and the Pseudo Wigner method (black solid curves). The white curves indicate the calculated frequencies using the zero-crossing method. The white curve near the center of the tube is the contour at 9 kHz given only by zero-cross method.

Download Full Size | PDF

 

Fig. 12. 3D velocity distribution of the cross section of the tube calculated from experimental values using the zero-crossing method.

Download Full Size | PDF

By performing the described data processing procedure on each A-scan from the cross sectional image, we can build a velocity distribution within this image. In Fig. 11, the contour frequency-shift plots are 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. The black solid curves indicate the frequency values of 2, 4, 6, 8, and 8.5 kHz calculated by the Wigner method. The white curves correspond to the frequency values of 2, 4, 6, 8, 8.5, and 9 kHz obtained with the zero-crossing method. This figure shows that, in the region near the glass wall where the velocity is relatively slow, the Wigner method produces results similar to the zero-crossing method. On the other hand, in the region near the center, where the velocity is relatively higher, the zero-crossing method can provide more precise results. This confirms the observations deduced from Fig. 10. Hence, the zero-crossing method is the one that gives the most precise velocity distribution estimation within the tube. For a better visibility, a 3D velocity distribution plot of the cross section of the tube obtained with the zero-crossing method is given in Fig. 12.

4.4 Susceptibility to noise

In Fig. 9, we showed that the zero crossing method gives results more precise than those with the pseudo-Wigner method at least in the higher Doppler-velocity regime. The question arises if this is also the case for lower signal-to-noise ratios. In Fig. 13(a), we display a portion of the normalized interferogram for an intralipid solution at rest in which there are lower SNR regions (in the vicinity of point number 36400) resulting from speckle effects or lower backscattering. The raw signal appears in black and the normalized signal, in blue. The “instantaneous” frequency is extracted using the different methods. Figure 13(b) shows the results prior to any filtering. The theoretical value for the frequency is indicated with a dashed line and corresponds to the mirror sweep. One immediately notices that the Hilbert method has an overall wider spread in frequencies (FM noise) than the Wigner and the zero-crossing methods. Actually, the zero-crossing method tracks the theoretical value very closely except at the points where there are rapid fluctuations in the vicinity of the zero-crossing points. In Fig. 14, we show the resulting statistics for the solution at rest (in red) and the solution in motion (in blue) corresponding to 3000 points in the center of the capillary, where the frequency is assumed to be nearly constant. Note that as the Doppler shift increases, the distribution widens, which is expected [35]. The distribution for the zero-crossing method is close to that of the pseudo-Wigner method and even narrower, except for a spill-over of frequencies at high values. For the fluid at rest, it is clear that the Hilbert method has a wider distribution. We then use the histograms of Fig. 14 to implement our maximum likelihood method whereby we reject all frequencies outside the FWHM (full width at half maximum) value of the distribution. This amounts to using an adaptive bandpass filter which follows the velocity distribution curve. By this filtering procedure, the spill-over effect of the zero-crossing method is strongly reduced. The rejection procedure is also helped by the fact that we normalize the signal, going from the black to the blue trace in Fig. 13(a). Indeed, upon rescaling, noisy transitions acquire higher slopes which translate into high frequencies which are filtered out more easily. In addition, the noise signal between the zero-crossing points never shows up in our procedure which means that in calculating the noise power to obtain the appropriate SNR value, one should integrate the noise power only over a sparse subset which puts in question the relevance of previous spectral-domain versus time-domain comparisons in our case.

 

Fig. 13. (a). OCT raw signal sample (black curve) and normalized signal sample (blue curve) used for frequency calculation. (b). Local frequencies calculated with different methods.

Download Full Size | PDF

 

Fig. 14. Normalized frequency distributions obtained with the three methods with fluid at rest (red) and with fluid in motion (blue).

Download Full Size | PDF

In order to get some feeling for the SNR sensitivity, we did several experiments changing the source power using a neutral density filter. In Fig. 15, we present the results for the 3000 points corresponding to the center of the tube, where the frequency is assumed to be nearly constant. The figure shows the FM SNR f_max/Δf (where f_max and Δf are the maximum frequency and the FWHM width of the frequency distributions of Fig. 14, respectively) as a function of incident optical power for the fluid at rest (plot (a)) and with the fluid in motion (plot (b)). It is clear that the zero crossing method almost always gives the best FM SNR, especially for high input optical power (i.e. high AM SNR). For lower optical power (i.e. low AM SNR), the zero crossing advantage compared to Wigner is not significant anymore. The FM SNR of the Hilbert method is almost always the worst one. This result shows that for acceptable AM SNRs, the zero crossing method tends to be less sensitive to noise than the Wigner and Hilbert methods. This observation holds whether there is flow inside the tube or not. The fluctuations in Fig. 15 are the results of individual scans which may arguably convey more information than the averages.

 

Fig. 15. FM SNR f_max/Δf at the center of the capillary tube calculated by the three methods (a) without and (b) with flow motion inside the tube.

Download Full Size | PDF

A similar comparison was done for data points close to the capillary wall and the results are shown in Fig. 16. Again the FM SNR turns out to favor the zero-crossing method most of the time in this low velocity regime for the experimental conditions we used, being followed closely by the pseudo-Wigner method, but the difference is not overwhelming. The Hilbert method looks inferior.

 

Fig. 16. FM SNR f_max/Δf at the edge of the capillary tube calculated by the three methods (a) without and (b) with flow motion inside the tube.

Download Full Size | PDF

 

Fig. 17. The frequency spread Δf at the edge of the capillary tube (a) and in the middle of the tube (b) as a function of the source power. Black line has a slope of -0.5. Dark brown dotted line represents the linear fit of the zero-crossing method results and has a slope of -0.29.

Download Full Size | PDF

Figure 17 shows the frequency spread at the edge of the tube and at its middle. This is to be compared with the usual $\frac{1}{\sqrt{\mathrm{SNR}}}$ dependency which is indicated by the black line. The dark brown dotted line is a linear fit of the results given by zero-crossing method. Instead of a -0.5 slope, we found a slope of -0.29. These figures are indicative of the fact that our method does not seem to follow precisely the straightforward phase noise analysis often found in the literature.

As a final point of comparison, we have measured the execution time to process the data for an A-scan, including the statistical maximum likelihood method. A word of caution must be given though. This is not the best achievable time because we used a possibly outdated machine running Matlab (Intel Pentium IV, 2.6 GHz processor, 1GB RAM). The numbers should be taken as a relative measure only. We obtained 48.9 seconds for the pseudo-Wigner method and 1.876 second for the zero-crossing method which is a gain factor of 26 in speed.

5. Conclusion

In conclusion, we presented a novel, simple and less time-consuming processing method based on the zero-crossing points for Doppler OCT. The performance of this method was compared with two other signal processing methods for the determination of the Doppler frequency shifts or velocities in an intralipid solution undergoing laminar flow in a rectilinear cylindrical tube. We found, based on experimental data, that among the techniques we investigated, the zero-crossing method is more precise globally than the Wigner method, especially in the higher velocity regime. A comparison with the commonly used phase-resolved Hilbert distribution method shows noticeable improvements in accuracy. The zero-crossing method should therefore be preferred for precision work. To the best of our knowledge, this is the first time that such a method is proposed for Doppler OCT. The next step will be to implement it in hardware thus offering a novel real-time data processing method in Doppler OCT systems.