1. Introduction

Optical coherence tomography (OCT) is an interferometric noninvasive high resolution imaging modality providing cross-sectional images with micrometer resolution [1]. Because of the high speed and increased signal-to-noise-ratio (SNR), Fourier domain OCT (FD OCT) is today preferred to time domain OCT (TD OCT) to acquire tomographic data [2,3]. There, the structural and velocity information is extracted by a Fourier transform (FT) of the spectral interference fringe signal resulting in a complex-valued depth scan. The absolute value is defined as A-scan representing one vertical line in a cross-sectional morphological image (B-scan). The argument contains the phase information of the light prevalently used for blood flow measurement in high scattering biological samples, which is based on the Doppler principle. There, sample movement in the direction of the incident sample beam results in a frequency shift of the interference fringe signal referred to as Doppler frequency shift. In TD OCT, the axial velocity can be measured by analyzing the fringe frequency or by analyzing the phase shift between sequential A-scans. Since FD OCT spectrally analyzes the interferometric signal, the local Doppler frequency shift at a special depth position can be obtained by the phase change between sequential A-scans known as phase-resolved Doppler OCT (DOCT). Thus, the depth-resolved Doppler information can be determined only from at least two spectral interference signals. This method was originally developed for TD OCT [4] and later adapted in spectral domain OCT (SD OCT) [5,6] – one of two implemented configurations in FD OCT. SD OCT uses a broadband light source and a spectrometer to record spectral interference fringes, whereas optical frequency domain imaging (OFDI), as second technical advancement in FD OCT, generally implies a narrowband wavelength-swept source and a single detector [7]. In OFDI, timing variations of the swept source relative to the acquisition electronics can degrade sensitivity by additional phase noise for the case of insufficient correction of timing induced phase errors. The sensitivity of the Doppler measurement is limited by the phase sensitivity of the OCT system, which is fundamentally affected by additive noise like detector and shot noise [7,8], by speckle noise resulting in a reduced SNR [8–10] and by decorrelation of adjacent signals due to transverse motion [11,12]. For biological samples like blood containing multiple scatterers, the stochastic interference of the phase terms of several scatterers at the detector results in multiplicative noise in form of speckle noise. The resulting inverse dependence of the DOCT precision on the image SNR is commonly known [8,13,14]. Transverse movement between beam and sample either by beam scanning or sample movement causes a decorrelation of the signal. A description of the statistical properties of the decorrelation noise was given in former publications [11,12]. Recently, S. Makita et al. have analyzed the statistics of phase-resolved DOCT showing the different contributions of additive, speckle, and decorrelation noise in single and dual-beam-scan DOCT highlighting the importance of dense oversampling in the case of single beam DOCT [10]. To minimize phase shift noise σ(Δφ) and by this velocity noise, averaging of amplitude-weighted phase shifts in the complex way [11,12,15] is essential. As shown by a former study of our group, an averaging of K > 5 samples achieves a reduced standard deviation of σ(Δφ) ~K^h with h about −0.5 for measurements with σ(Δφ) ≤ 1 rad [12]. Chan et al. [16] compare the Kasai autocorrelation with a maximum likelihood estimator for Doppler OCT. The Kasai autocorrelation [17] in its discrete form is equivalent to DOCT. This model takes only additive and multiplicative noise into account, neglecting the speckle and decorrelation noise, being dominant in many biomedical applications.

Since the phase sensitivity decreases with the reciprocal value of the square root of the SNR of the measured signal [8], OCT applied in the 800 nm wavelength range in combination with an SD OCT configuration is particularly vulnerable to high Doppler noise in blood flow measurements. The reason for this is first a low SNR effected by motion-induced signal blurring [18,19] and second a low backscattering with increasing depth in large blood vessels due to high scattering at erythrocytes in the 800 nm wavelength range. Since both influencing factors generate a strong signal damping, different approaches for quantitative DOCT have been proposed for measurements with low SNR [20–22]. Most promising is the joint spectral and time domain OCT (jSTdOCT) presented by M. Skulmowski et al. [21,23,24], where interferometric OCT signals were analyzed in time and frequency domain by two independent Fourier transformations. By this method, the velocity estimation is determined directly by the Doppler frequency shift without the explicit phase information. It was proposed, that jSTdOCT is less sensitive to phase instabilities present in measurements with low SNR. Although, the results of frequency-based methods and phase-based methods have been compared in several works, to date, the mathematical relation between jSTdOCT and phase-resolved DOCT is neither described in literature nor considered in the research field of Doppler velocity detection by OCT. This work overcomes the lack and introduces first a common model of both established algorithms. Additionally, a modification of jSTdOCT is proposed for effective velocity noise reduction. The enhanced algorithm enhjSTdOCT takes for the first time the influence of noise into account and shows advantages in terms of precision and SNR theoretically and experimentally.

2. Theoretical perspectives

In spectral domain OCT (SD OCT), the spectral interference signals detected by the line scan camera of the spectrometer provide the depth-resolved structural information of the sample in the optical beam. In general, the measured raw interference spectra were corrected for background, spectrally shaped, transformed to equidistant in wave number and finally Fourier transformed to get the backscattering amplitude of the light reflected from the sample as a function of the depth, which is called A-scan. After the Fourier transformation, the depth-encoded complex OCT signal from one single moving sample reflector in SD OCT can be given by:

2.1 Principle of phase-resolved Doppler OCT (DOCT) by using complex A-scans

Although the phase term φ_j(z) of the complex A-scan in Eq. (1) is generally random for biological samples and has no direct relevance, the functional information of flow can be obtained by measuring the phase change relative to a reference signal to obtain structural changes in depth well below the wavelength [4]. For phase-resolved Doppler OCT (DOCT), any depth scan according to Eq. (1) that is detected temporally earlier at the same transverse position of the sample can be used as a reference. Therefore, DOCT requires minimum two adjacent A-scans of identical or transversely overlapping areas in the tissue to measure the axial velocity component v_z toward the incident sample beam. Though, any spatial displacement between the single measurements, due to less oversampling or sample movement in transverse direction, will cause a decorrelation of the phase signal detected [11]. Depending on the amount of overlap of the sample structure, the decorrelation noise degrades the flow sensitivity of the system and results in noisier phase measurements [12].

Because the signal reconstruction in SD OCT provides directly the phase information of the OCT signal, the phase change between two successive A-scans, related to the Doppler frequency f_D(z), is commonly calculated in the spectral domain. In the following, we will consider specifics of DOCT for SD OCT-based systems. Although, the calculation of the phase shift by consecutive absolute phase values can be found in former literature, we chose Eq. (2) to determine the phase difference

2.2 Principle of joint spectral and time domain OCT (jSTdOCT) and its enhancement

Alternatively to DOCT, the flow velocity can be estimated by using joint spectral and time domain OCT (jSTdOCT), which calculates the Doppler frequency shift f_D(z) occurring between consecutive spectral interference fringes [21]. As f_D(z) is related to Δφ(z) (cf. Eqs. (3) and (5)), jSTdOCT is a variation of the conventional DOCT, that determines the temporal frequency shift instead of the phase shift of sequentially detected A-scans. There, jSTdOCT is based on the 2D Fourier transformation (2D FT) of the interference spectra highly temporally and spatially oversampled. For the case of an M-scan, meaning A-scans are detected at identical transverse position of the sample, the spectral fringes depend on wave number k and time t as presented in Fig. 1(a).By using the 2D FT, the wave number domain is transferred directly to the depth z as usual for conventional OCT and the time domain is converted into the Doppler frequency f_D. While the order of the two FTs does not matter, the way of performing the first FT vertically in the direction of k (green vertical line in Fig. 1(a)) and then applying the second FT horizontally over a limited range of adjacent complex A-scans taken at the same (M-scan) (red horizontal line in Fig. 1(b)) or nearly the same location (B-scan) is used in this manuscript and has the advantage that the correlation between jSTdOCT and DOCT can be analyzed because of the same input data Γ_j(z) used for the axial velocity estimation. The result of the second FT is the amplitude of the different Doppler frequency shifts f_D in a certain depth z of the sample (Fig. 1(c)). Since the Doppler frequency shift is not a single frequency peak but a broadened frequency spectrum in the range of [-f_A-scan/2,f_A-scan/2], the value of f_D(z) is determined by the component with the maximum amplitude [21]. Components that obviously belong to depths without signal information were eliminated by a threshold method. Knowing the center wavelength λ₀, the axial flow velocity v_z(z) is then calculated from f_D(z) as shown in Eq. (5).

 

Fig. 1 One way to apply the 2D FT in joint spectral and time domain OCT (jSTdOCT) as used in this manuscript; (a) 2D interferogram consisting of 704 interference spectra of the detected 1% Intralipid flow within a 315 µm glass capillary imaged at a Doppler angle β of 81.7°. (b) Grayscale structural image (M-scan with 704 A-scans) containing the amplitude information of the capillary center with an Intralipid flow rate of 1.3 ml/h. (c) Grayscale image of the Doppler frequency shift f_D distribution in depth z extracted from 704 complex-valued A-scans.

Download Full Size | PDF

The unique measurement range of jSTdOCT and DOCT depends on the time interval T_A-scan = f_A-scan⁻¹ between two A-scans as presented in Eq. (6), where the plus-minus sign indicates the bidirectional flow velocity measurement. For absolute axial sample velocities larger than the maximum v_z, |Δφ| exceeds |π| in DOCT and |f_D| goes beyond the Nyquist limit of |f_A-scan/2| in jSTdOCT.

For further analyses, the velocities v_x in transverse and v_z in axial direction are transformed to dimensionless units v_x’ and v_z’ by v_x’ = v_x/(f_A-scanw₀) and v_z’ = v_z2n/(f_A-scanλ₀), where w₀ is the beam width (FWHM) of the Gaussian sample beam at the focus position [12,26]. By scaling the Doppler frequency to the A-scan frequency, the unique measurement range of the normalized Doppler frequency shift f_D’ = f_D/f_A-scan is in the interval between [-0.5,0.5]. Consequently, f_D’ and v_z’ are equal within the unique measurement range. Because the second FT is applied to the complex analytical result Γ_j(z) of the first FT (cf. Eq. (1)), where the input sequence consists of K samples typically 30 to 40, the output of the second FT has in total K different (positive and negative) frequencies that span the range of f_D’. To quantify the Doppler frequency shift f_D’ from the broadened Doppler frequency spectrum, M. Szkulmowski et al. proposed to determine the frequency with the maximum intensity [21], which defines one velocity value for each depth z, as presented by Eq. (7), where the index j is again the time sample coordinate for the case of an M-scan (cf. Fig. 1(b)).

 

Fig. 2 (a) Normalized amplitude of the Doppler frequency spectrum caused by a single obliquely moving particle passing the sample beam center in the middle of the detection period with K measurements. Simulated data (K = 32) in blue without and in red with zero padding by using 480 elements. (b) Same diagram for two scatterers of equal amplitude separated by 10 time steps with a phase difference of 2/3π. The signal of both scatterers declines at the border of the interval. (c) Doppler frequency spectrum of a single scatterer in the middle of the sample beam at the beginning of the time window with K samples. (d) Linear frequency spectrum of a large number (100) of point scatterers with random phases and amplitudes passing the beam at different points of time. The sample velocities are in all simulations set to v_x’ = 0.1 and v_z’ = 0.1. In all cases the maximum amplitude after zero padding was normalized to 1. The green vertical line shows the set axial velocity, which coincides with the result of DOCT in (a) and (c) and exhibits small deviations to the result of DOCT (b) vx’ = 0.089 and in (d) vx’ = 0.088. Only for (a) and (c), the maximum of f_D’ after zero padding is identical to the set axial velocity component.

Download Full Size | PDF

In practice, a signal of many scatterers with different amplitudes and different phases randomly distributed exists, which will change the simulated curve in Fig. 2(a) in two ways. First, the signal from different particles passing the sample beam at different time points and with different phases will cause an amplitude modulation within the Gaussian peak of the Doppler frequency distribution. In Fig. 2(b), the simulation of two scatterers passing the center of the sample beam ten time steps apart with equal amplitude and a phase difference of 2/3π show a modulated shifted Doppler frequency spectrum. As for the single particle in Fig. 2(a), both scatterers run into and out the sample beam within the window of K measurements. The second change results from the rectangular window in time, truncating the signal of some scatterers, which causes a broadening and an additional offset in the Doppler frequency spectrum. This is shown in Fig. 2(c) at first for a single particle positioned at the center of the sample beam at the beginning of the measurement window. The effect of many particles (N = 100) passing the beam at random time points with random phases and amplitudes is simulated in the data presented in Fig. 2(d). The occurring offset in the image parts (c) and (d) can be reduced by applying a window function to the data prior to the second FT but this broadens the center, which should be avoided.

On closer examination of the previous results, estimating the velocity from the maximum of the Doppler frequency distribution has still the advantage, that it is not biased by noise nor by parts of f_D’ wrapped to the opposite measurement range for velocities exceeding the Nyquist limit but, as seen from Fig. 2, the maximum will in general not coincide with the true axial velocity component leading to inaccurate or noisy velocity profiles. Zero padding may reduce the error and the maximum will be typically within the half width of the Gaussian envelope of the Doppler spectrum. In order to reduce the noise of the axial flow velocities, we propose to measure the center of gravity of the frequency spectrum e.g. the mean frequency as an enhancement to the classic jSTdOCT [28,29]. But as shown for phase data in DOCT [15], the simple averaging of phase differences Δφ itself will cause a bias toward zero velocity for values of Δφ close to ± π. To avoid this effect, we suggest to distribute the data on a circle by multiplying each component of the power spectrum of f_D’ with a factor of exp[i2πm/K] corresponding to the angle in the complex plane, where m is the index of the components of the second FT and K is assumed to be even, and sum over all newly formed complex values (cf. Fig. 3). The phase of the sum, which is defined as Φ, is then transformed back to f_D’ or v_z’ scale, respectively.

 

Fig. 3 Graphic of the Doppler frequency spectrum transferred on the circle (a) and onto the complex plane (b) by drawing each component of the power spectrum as a function of exp[i2πm/K]. The components from the power spectrum (vertical) are the weights for the complex values on the circle. The numerically simulated Doppler frequency spectrum of multiple reflecting particles passing the sample beam, shown in Fig. 2(d), is used. The blue line indicates the angle Φ calculated by Eq. (8).

Download Full Size | PDF

This way of averaging has the advantage, that it does not bias the data if parts of the distribution of f_D’ are above the Nyquist limit. Moreover, the average is not disturbed by additional noise floor as constant values not belonging to flow information multiplied by the complex exponential factor add up to zero. Besides the normalization, this operation is equivalent to the calculation of the Fourier coefficient of the first harmonic by an inverse FT from the power spectrum. We define this as FT⁻¹(P(ω))₁ where P(ω) is the power spectrum of Γ_j(z). Finally, the axial sample velocity is calculated by determining the argument of the first Fourier coefficient, indicated by the index 1, of the inverse FT shown in Eq. (9) and transferring it to the normalized velocity v_z’ as described by the upper part of Eq. (8).

2.3 Common theoretical model of enhanced jSTdOCT and phase-resolved Doppler OCT

On the basis of the enhanced jSTdOCT algorithm (enhjSTdOCT) proposed in the previous section, the power spectrum of FT[Γ_j(z)] in the inner part of Eq. (9) can be rewritten as product of the FT of the complex analytical signal Γ_j(z) and the FT of its conjugate complex form presented in Eq. (10), which is equivalent to Eq. (11).

Consequently, enhjSTdOCT, averaging the frequency spectrum in the non-biasing way described by Eq. (8) instead of taking simply the maximum signal, is equivalent to DOCT if the phases are averaged complex described by Eq. (4). As the algorithms of the enhjSTdOCT and DOCT are mathematically equal, all features attributed to one of the methods will be valid for both. This means for instance that the velocity deviations found in spectral domain OCT (SD OCT) due to the finite detector exposure time [12] will be the same in enhjSTdOCT. Calculations with both methods on simulated and experimental data lead to the same result within the numerical accuracy of the computation.

2.4 Optimal exponent in enhjSTdOCT

The velocity estimation by Eq. (7) for jSTdOCT, taking the argument of the maximum amplitude of FT[Γ_j(z)], is equivalent to the following equation.

 

Fig. 4 (a) Standard deviation of the normalized axial velocity σ(v_z’) calculated as a function of the exponent Ex in Eq. (15) for a transverse velocity component of v_x’ = 0.1. The diagram shows an optimal value for Ex = 1. σ(v_z’) for Ex = 2 is only larger by 17%, while the limiting value for large exponents is higher by a factor of more than 2.5. The position of Ex = 1 and Ex = 2 are marked by vertical black lines. (b) σ(v_z’) as a function of the transverse velocity v_x’ for Ex = 1 (green graph), Ex = 2 (corresp. to DOCT, red graph) and Ex = 10⁷ (corresp. to jSTdOCT, blue graph). Ex = 1 is optimal over the total range of transverse velocities. At high transverse velocities, the results with different exponents converge because all values approach the standard deviation for an isotropic distribution of velocities in the range of ± 0.5 which is about 0.289.

Download Full Size | PDF

From a statistical point of view, averaging data is optimal when the noise in the data is comparable. Unfortunately, the noise in the addends of Eq. (15) is not simply calculated as the noise in the original data Γ_j(z) as a function of j is not white noise. Assuming a white noise in the original data, the noise of the amplitudes of the Fourier coefficients is equal leading to the preliminary conclusion, that simply the amplitudes of the Doppler frequency spectrum and with it an exponent of Ex = 1 in Eq. (15) should be optimal as found for a large range of transverse velocities in the simulation. As the results of the numerical simulation in Fig. 4 are based on an ideal noiseless signal, the question rises if the SNR has an effect on the result. Therefore, σ(v_z’) was calculated as a function of the transverse velocity component v_x’ and the reciprocal of the signal-to-noise ratio of the measured OCT signal (1/SNR). The result is presented in Fig. 5. For the calculation, white Gaussian noise is simply added to Γ_j(z) to simulate lower OCT signals with proportionate noise components. The role of decreased SNR is of particular interest since low SNR occurs especially for blood flow imaging with spectrometer-based OCT systems due to fringe washout at high axial flow velocities [18,19] and also for systems designed for the 800 nm wavelength band because of strong scattering red blood cells and the resulting limited penetration depth of light. In our simulations and measurements, the SNR is defined as ratio of the mean amplitude of the A-scan and the mean amplitude of the noise floor. The axial scaling with 1/SNR is chosen in Fig. 5 because useful parameters are then in the range of zero to one instead of one to infinity for the SNR. The chosen range of the normalized transverse velocity v_x’ from 0.0 to 0.6 is representative for measurements in single-beam DOCT, where it often outperforms dual-beam DOCT [10], whereas values of v_x’ larger than 0.6 are more theoretical. Input data sets with mean amplitude less than 30% above the noise level were discarded. The number of data sets discarded becomes significant for SNR values below 1.5 and reaches 2/3 at an SNR of 1. Without this removal, especially at low transverse velocities, the phase noise increases dramatically. In Fig. 5, the standard deviation σ(v_z’) is presented in the negated way just due to a better visualization of the results (small noise at the top). For each exponent (Ex = 0.5, 1, 2, 4, 8, 16, ∞), σ(v_z’) is shown by a single colored 3D plane described by the legend on the right. As the phase noise depends only smoothly on the exponent, a finer graduation is not needed. Moreover, smaller or larger exponents (i.e. Ex = 32, Ex = 64) reveal no benefit. For better comprehensibility, the parts (a) to (d) of Fig. 5 show different perspectives of the same 3D plot with identical mathematical content.

 

Fig. 5 Standard deviation of the axial velocity σ(v_z’) calculated as a function of the transverse velocity v_x’ and the reciprocal of the signal-to-noise ratio of the measured OCT signal 1/SNR. For a better visualization, σ(v_z’) is negated. Seven specific exponents were used for the calculation of the respective σ(v_z’), each presented by a single colored 3D plane. The result of jSTdOCT with Ex = ∞ corresponds to the red data and the one of DOCT with Ex = 2 to the cyan. Ex = 0.5 (yellow), Ex = 1 (green), Ex = 2 (cyan), Ex = 4 (blue), Ex = 8 (violet), Ex = 16 (magenta) and Ex = ∞ (red). For better comprehensibility, the parts (a) to (d) of Fig. 5 show different perspectives of the same 3D plot with identical mathematical content. Furthermore, the relative noise of DOCT and jSTdOCT compared to the enhjSTdOCT with the optimized exponent is shown in (e) and (f), respectively. The red color signals optimal or nearly optimal performance while other colors show higher phase noise. Note the different scale and color scale in (e) and (f).

Download Full Size | PDF

Starting with the most obvious result in Fig. 5(a), one can identify a quasi new contour plot in the plan view, formed by the single planes of the selected exponents, which shows the optimal exponents for a specific range of transverse sample velocity and signal quality. For the entire measurement range with v_x’ = [0,0.6] and 1/SNR = [0,1], the exponents with the minimal possible σ(v_z’) are increasing with decreasing v_x’ and worsen SNR. For almost purely axial motion and decreasing SNR, exponents greater than two towards infinite give sample velocities with smallest variation. Towards larger v_x’, higher than 0.2, an exponent between 0.5 and 2 shows the lowest dispersion of the calculated velocity values. In this range, the planes of Ex = 0.5 and Ex = 1 nearly touch each other. On the basis of part (b) to (d), σ(v_z’) is generally increasing with higher v_x’, worsen SNR of the backscattering signal and higher exponents for enhjSTdOCT calculation, as expected. Contrary to expectations and to the prior at small SNR, low v_x’ leads to an increase in σ(v_z’) at small exponents due to the low speed of the speckle. The total signal may be hidden in noise leading to arbitrary velocities. In spite of the large number of randomly selected distributions, this might be also the reason for the frayed border between Ex = 4 and Ex = 8 planes. For an ideal signal with 1/SNR = 0, σ(v_z’) maximal amounts to 0.03 for v_x’ = [0.0,0.6]. Staying with this noiseless signal, an exponent larger than two has no advantage over smaller exponents. On the contrary, an exponent converging to infinity (shown by the red plane in Fig. 5), as used in jSTdOCT by taking the argument of the maximum, results in a large noise of the calculated velocity values, especially for measurements with a high transverse velocity component. This plane with Ex → ∞ shows that the maximum intensity signal of the broadened f_D’ gives generally a higher variance of the evaluated velocities except the case that the sample motion is almost axial and the backscattering signal is very low (Fig. 5(c)). For an application with exclusively axial sample movement, Ex → ∞ shows a smaller σ(v_z’) as Ex = 0.5, Ex = 1 and Ex = 2 for the case of a weak signal. Additionally to the absolute noise in Fig. 5(e) and Fig. 5(d) the relative performance of DOCT and jSTdOCT compared to enhjSTdOCT is shown. While at high transverse velocities, DOCT is only worse by up to 40% compared to the optimal exponent (green area on the left of Fig. 5(e)), at nearly axial movement and small SNR, the noise might be nearly four times the noise encountered with the optimal exponent of Ex = 8 (white area beyond the scale of Fig. 5(e)). Note, that the phase noise at almost axial movement is still very small in DOCT. Vice versa jSTdOCT is nearly optimal at purely axial movement (red area at the bottom of Fig. 5(d)), while the phase noise might be worse by a factor of more than three at high transverse velocities (blue area at the top of Fig. 5(d)).

3. Experimental validation

3.1 Data measurement and processing

To validate the theoretical model, the short wavelength band of the dual-band spectral domain OCT system centered at 800 nm with a full spectral width of 200 nm is used with the fiber-coupled customized scanning unit with free-space Michelson interferometer [30]. The light source of the system is a supercontinuum laser light source (SuperK Versa, Koheras A/S, Denmark), which is spectrally shaped to two wavelength ranges of 700-950 nm and 1100-1400 nm. The 800 nm-band used in this study allows a 3D imaging with a high axial resolution better than 3.8 µm in air and a lateral resolution of 7 µm. To guarantee a high backscattering intensity and with it a strong SNR of the sample, the adjustable A-scan rate was set to 12 kHz for the flow measurement. The flow phantom used consists of a syringe pump (Injectomat cp-PS, Fresenius Kabi) delivering a highly scattering and low absorbent 1% Intralipid emulsion through a glass capillary (Paul Marienfeld GmbH&CoKG, Typ 2930203) with 315 µm inner diameter. The set flow rates of the pump ranged from 0.5 to 2.3 ml/h corresponding to mean flow velocities of 1.8 to 8.3 mm/s. The center of the capillary was positioned at the sample arm focus. The Doppler angle β (between flow direction and optical beam) of 81.7° was measured by means of a volume scan. For the quantification of the flow velocities by enhjSTdOCT, M-scans consisting of 704 A-scans were detected at each lateral position of the glass capillary cross-section.

The exemplary results of the measurement at a flow rate of 1.3 ml/h are presented in Fig. 6. The absolute flow velocity was measured by DOCT to be 9.2 mm/s at the capillary center corresponding to a normalized transverse velocity v_x’ of 0.11 and in axial direction to v_z’ of 0.37 in agreement with the theoretically calculated flow velocities. For the noise analysis of the flow velocity, the center M-scan consisting of 704 A-scans was processed by enhjSTdOCT with varying exponent. Figure 6(a) and 6(b) present the classic structural OCT image after the first FT without and with Doppler flow information overlayed (DOCT by Eq. (4)). As for jSTdOCT, the 704 complex A-scans, as result of the first FT, are fed into the second FT (Fig. 1) to calculate the depth-dependent Doppler frequency spectra of the Intralipid flow (Fig. 6(c)). By using the proposed enhjSTdOCT algorithm, the flow velocities can be calculated depending on the exponent set to weight the broadened Doppler frequency spectrum (cp. Eq. (15)). The results of four different exponents (Ex = 1, 2, 8, Ex → ∞) are presented exemplarily in Fig. 6(d) with the result that for this example lower exponents are fortunate. To validate the performance of different exponents for the velocity calculation by means of the standard deviation of the calculated Doppler frequency shift, the flow velocities were additionally calculated from small packages of 32 interference spectra resulting in 22 single Doppler flow profiles when using the entire M-scan of 704 A-scans (Fig. 6(e)). In Fig. 6(f), the standard deviation of the normalized Doppler frequencies σ(f_D’) = σ(v_z’) at each depth position of the capillary center A-scan is shown confirming the strongest noise of the velocity for exponents tending to infinity.

 

Fig. 6 Stepwise processing of the 1% Intralipid flow measurement with an Doppler angle β of 81.7° and an absolute flow velocity of 9.2 mm/s corresponding to v_x’ = 0.11 at the capillary center. There, M-scans consisting of 704 A-scans were acquired with f_A-scan = 12 kHz at the glass capillary cross-section achieving a measured SNR of 14.9 dB at the capillary center. Backscattering information (a) and overlayed Doppler flow information (b) by usual FD OCT and DOCT processing, respectively, is presented. The depth-resolved Doppler frequency spectra (c) are determined by K = 704 complex A-scans applied to the 2. FT of the enhjSTdOCT algorithm. Doppler flow profiles are calculated as a function of the depth position for four different exponents (d). The direction of time (a and b), frequency (c and e) and depth (a, b, c and e) are indicated. To calculate the standard deviation σ(f_D’), 22 input sequences of K = 32 were defined to be Fourier transformed in time (e,f). The exponent resulting in the lowest velocity noise can be read from the black dot in the 3D surface plot (g) being Ex = 1 and Ex = 2, respectively, for v_x’ = 0.11 and 1/SNR = 0.18. The standard deviation σ(f_D’) as a function of the exponent Ex for v_x’ = 0.11 confirm the numerically simulated optimal exponent for the capillary center (h) (red measurement, black theoretical curve). The positions of Ex = 1 and Ex = 2 are marked by vertical lines.

Download Full Size | PDF

With v_x’ of 0.11 at the capillary center and an SNR of 14.9 dB (1/SNR = 0.18), the optimal exponent for this flow measurement can be read from the black point in the 3D plot in Fig. 6(g). There, the minimal noise of the flow velocity is achieved comparably by an exponent of one and two. In order to proof this, the measured standard deviation σ(f_D’) is plotted as a function of the exponent in Fig. 6(h) (red curve) in comparison to the black theoretical curve for v_x’ of 0.11. There, a small offset of σ(f_D’) of 0.004 can be seen, which could be caused by additional intrinsic statistical phase noise since the light of the incident sample beam already passes Intralipid solution before scattering at particles at the capillary center, where σ(f_D’) was measured for the diagram in Fig. 6(h). Other causes for the deviation might be an erroneous estimate of the SNR, v_x’, neglected diffusion or other simplifications in the model. Despite of the small deviation, the measured optimal exponent of Ex = 1 and Ex = 2, respectively, confirm the expectation of the numerical simulation.

3.2 SNR decrease due to sample movement

Since a dependency of the axial velocity noise on the SNR was calculated theoretically, the spectrometer-induced signal power decrease due to sample motion was induced with the flow phantom model to present the effect of different exponents in Eq. (15). Therefore, a measurement series was performed with maximum flow velocities of the Intralipid solution at the capillary center ranging from 3.5 to 16.6 mm/s corresponding to v_x’ of 0.04 to 0.2 (cf. Fig. 7).

 

Fig. 7 Results of the measurement series for presenting the dependency of the optimal exponent Ex on the transverse velocity component and the SNR. Image series show the intensity (a) and flow information (b) of the Intralipid flow at the capillary center A-scan for flow transverse velocities v_x’ ranging from 0.04 to 0.2. Corresponding diagrams showing the Doppler frequency shift f_D’ and the velocity noise σ(f_D’) are presented in (c) and (d). The optimal exponent for the capillary center can be read from image parts (e) and (f).

Download Full Size | PDF

The image series in Fig. 7(a) shows the signal power damping with increasing flow velocity. The interference spectra detected were processed as described above to generate the broadened Doppler frequency shift as a function of the depth z (image series in Fig. 7(b)) and to calculate the normalized axial velocity f_D’ = v_z’ as well as the corresponding standard deviation σ(f_D’) (diagrams in Fig. 7(c) and 7(d)). As seen by the four exemplary measurements, the noise of the flow velocity profiles is strongly increasing with rising exponent and higher signal decrease. For the lowest flow, the first diagram in the series in Fig. 7(d) shows the lowest velocity noise for an exponent of two corresponding to the classic DOCT as proven in subsection 2.3. By means of the three remaining diagrams (Fig. 7(d)), there is no difference between σ(f_D’) achieved with Ex = 1 and Ex = 2 for higher flow velocities. Therefore, the 3D plot in Fig. 7(e), showing the theoretical σ(f_D’) = σ(v_z’) for different exponents, offers the optimal exponent for the corresponding flow situation marked by the four symbols (asterisk, plus sign, point, multiplication sign). By Fig. 7(f), the trend of σ(f_D’) as a function of Ex is presented in detail. There it is confirmed, that for v_x’ = 0.04 the optimal exponent is Ex = 2, whereas for v_x’ = 0.09, the minimal σ(f_D’) is achieved simultaneously with Ex = 1 and Ex = 2 and for higher flows with v_x’ = 0.13 and v_x’ = 0.20, an exponent of Ex = 1 gives the best result.

3.3 Reduced SNR due to additional Gaussian white noise

For the experimental investigation of the impact of the SNR on the axial velocity noise σ(f_D’) independent of the sample velocity, Gaussian white noise was added to the depth-encoded complex OCT signal Γ_j(z) (Eq. (1)) of the flow measurement with a maximum velocity of 7.1 mm/s corresponding to v_x’ of 0.09 and v_z’ of 0.29 at the capillary center. The original measurement is already presented in Fig. 7 with a measured SNR of 15.3 dB at the center of the capillary. For the analysis, the noise in the original data was increased gradually. In this way, the SNR is reduced to 10.0 dB, 7.0 dB, 5.5 dB, 4.5 dB and 3.5 dB as presented in Fig. 8(a).The corresponding broadened Doppler frequency spectra are shown in the image series in Fig. 8(b). The Doppler flow profiles of the original data with 15.3 dB (1/SNR = 0.17) (c) are compared to the one with 3.5 dB (1/SNR = 0.67) (d) for four different exponents (Ex = 1, 2, 8, Ex → ∞). It can be seen, that the noise of the flow profiles is strongly increasing for small exponents Ex of one and two for the measurement with increased noise level (d) in comparison to the original data (c). Additionally, an exponent Ex = 8 results in the flow profile with lowest noise. This visual analysis is confirmed by the trend of σ(f_D’) as a function of the exponent Ex in Fig. 8(e), where the optimal exponent is increasing with decreasing OCT signal from 1/SNR = 0.17 to 1/SNR = 0.67. In a direct comparison with the results of the calculation in Fig. 8(g), it becomes apparent that the optimal exponent is slightly higher than theoretically predicted. Exemplarily, for an SNR of 7.0 dB (1/SNR = 0.44), the optimal exponent experimentally found is Ex = 4, whereas the model proposes Ex = 2. These characteristics continue in regard to the optimal exponent for this flow measurement. Accordingly, the optimal exponent for SNR = 3.5 dB (1/SNR = 0.67) is experimentally Ex = 8 instead of Ex = 4. This divergence might be caused by the measurement itself. The statisticsof our simulation shows, that for the repetition of the model calculation for identical parameters v_x’ and 1/SNR, the optimal exponent varies slightly in spite of many possible distributions of the large number of moving scatterers. Consequently, our model is calculated only for a medium resolution in the 3D diagrams in Fig. 5 to obtain smooth borders between planes of adjacent optimal exponents. A higher resolution would result in even more frayed borders due to the statistics. Since only one distribution of many Intralipid particles is measured in our experiment with 704 A-scans, only one possible constellation is detected, which can indeed deviate from the theoretical plot. Additionally, a small offset in σ(f_D’) is seen again (cp. Fig. 6(h)). In Fig. 8(f), σ(f_D’) is shown as a function of the reciprocal value of the SNR. There, it becomes apparent, that for a backscattering signal with a high SNR, smaller exponents of Ex = 1 and Ex = 2 are again preferable over higher ones. Going to lower signals with strong signal damping, smaller exponents Ex < 4 result in a higher velocity noise in comparison to higher ones. Interestingly, for this measurement with a maximum v_x’ = 0.09 at the capillary center, σ(f_D’) is almost constant for Ex = 8, Ex = 16 and Ex → ∞ over the presented 1/SNR range. The trend confirms the theoretical σ(f_D’) (cp. Fig. 5(c)) despite of a small absolute difference.

 

Fig. 8 OCT intensity images (a) and broadened Doppler frequency shifts (b) of the representative measurement with v_x’ = 0.09 (cp. Fig. 7), where Gaussian white noise is added to Γ_j(z) to generate incrementally lower backscattering signals of the same measurement. Exemplarily, flow profiles are presented for the original data with SNR = 15.3 (c) in comparison to the attenuated one with SNR = 3.5 (d) for four different exponents (Ex = 1, Ex = 2, Ex = 8, Ex → ∞). As seen, (c) offers an optimal exponent of Ex = 1 and (d) of Ex = 8. The diagrams in (e) and (f) present σ(f_D’) as a function of the exponent Ex and the signal quality 1/SNR for exponents ranging from Ex = 0.5 to Ex → ∞ and 1/SNR from 0.17 to 0.67. The theoretical 3D diagram shows the estimated optimal exponents for the flow measurement with different signal amplitudes (marked by specific symbols) due to additive Gaussian white noise in the post processing.

Download Full Size | PDF

4. Conclusion and outlook

In this study, the relation of the enhanced jSTdOCT (enhjSTdOCT) and phase-resolved Doppler OCT (DOCT) was analyzed. The Doppler spectrum for a single scatterer and a large number of point scatterers passing the beam on a ramp was calculated thereby explaining the noise encountered in jSTdOCT. The presented enhanced algorithm of enhjSTdOCT permits a more precise flow velocity measurement in comparison to the conventional jSTdOCT. There, the Doppler frequency shift f_D is not determined by detecting the maximum intensity signal of the broadened f_D, but is measured by calculating the center of gravity via the complex analytical signal as a result of the second FT. By transforming the enhjSTdOCT algorithm to a common mathematical notation with DOCT using identical input data Γ(z), it becomes apparent that both methods are mathematically identical. While DOCT is computationally faster, the benefit of using jSTdOCT is the direct access to the velocity spectrum. In the case that more than one source with different axial velocities at the same depth contribute to the OCT signal, both velocities could be estimated with jSTdOCT while DOCT would only provide an average value with a higher value of the scattering of the individual phase differences.

The notation of jSTdOCT and enhjSTdOCT can be brought to a common notation with only different exponents of infinity in the case of jSTdOCT versus two for enhjSTdOCT leading to the question which exponent is optimal. Numerical simulations, taking the statistics of transversely moving densely packet point scatterers and the limited SNR into account, show a minimal phase error for exponents between 0.5 and 8 with smaller exponents at high SNR and low transverse velocity, and higher exponents at low SNR and nearly axial motion. The proposed method of enhjSTdOCT with an optimal exponent outperforms DOCT up to a factor of nearly four at mostly axial movement and jSTdOCT up to a factor of more than three at high SNR and high transverse velocity within the simulated region of SNR between one and infinity and normalized transverse velocity between 0 and 0.6 (Fig. 5(e), 5(f)).

The experimental comparison of enhjSTdOCT with different exponents affirms these properties in general with slightly higher exponents than numerically calculated. For a practical implementation, the SNR and a raw guess for the axial velocity could be computed in a first step. In most cases, a 3D measurement will give the direction of the flow and thereby the transverse velocity component. With a rough estimate of the transverse velocity and the SNR, Fig. 5 gives the optimal exponent for enhjSTdOCT, which would be applied in a second step for a more precise velocity calculation.

Our model does not take the axial decorrelation into account. In most applications with medium resolution OCT systems, this factor will be much smaller than the transverse decorrelation and may therefore be neglected. If in the case of high axial velocity the axial decorrelation should be considered, the total decorrelation could be accounted for by taking the RMS value of both factors, the transverse v_x’ and the axial v_z2n/(f_A-scan∙δz), with the axial resolution δz of the system instead of the pure normalized transverse velocity v_x’, although the window functions in transverse and axial directions do not have the same form.

Though the data presented were acquired with a spectrometer-based system, the result will be valid for optical frequency domain imaging (OFDI) as well. Actually, the simulated data did not take the fringe washout due to finite exposure time in spectrometer-based systems into account. Known artifacts of velocity measurements by OCT, as fringe washout in SD OCT, and axial image shift in OFDI due to axial movement [19], will not be influenced by the new algorithm. Effects of multiple scattering may change the frequency spectrum and thereby the calculated velocity might depend on the exponent to a certain extend.