1.

Introduction

Optical coherence tomography (OCT) is a high resolution imaging technique.^{1, 2, 3} With more than a decade of development, it has proven to be a useful imaging modality for many biological applications. Doppler OCT (DOCT) or optical Doppler tomography (ODT)^{4, 5, 6, 7} is a functional extension of OCT with great potential to explore subsurface blood flow information. ODT estimates the depth-resolved frequency or phase changes in addition to the backscattering intensity from OCT raw data. The development of suitable estimation algorithms is therefore a key issue. The most straightforward approach is to use the short-time Fourier transform (STFT) to find the spectral centroid in successive time windows.^{4, 5, 6, 7} Extensions of STFT include adaptive centroid⁸ and weighted centroid methods.⁹ The sliding-window filter technique implements a moving digital filter bank in the time domain instead of the Fourier transform to estimate the local spectral maximum.^{9, 10} Our previous studies have shown advantages of a sliding-window filter technique, such as improved noise robustness compared with several other techniques. While these techniques have been found to be effective for ODT signal processing, the use of a finite time window to localize the frequency results in coupled frequency sensitivity and spatial resolution. For example, a longer time window enhances the frequency sensitivity at the expense of reduced spatial resolution. These algorithms are also unable to achieve simultaneous fine velocity sensitivity and fast frame rate, both of which are important in biomedical applications.

Chen ¹¹ and Zhao ¹² have pioneered a phase-resolved technique, which utilizes the analytic continuation of the OCT coherence function to obtain the phase term via the Hilbert transform. This technique can be readily applied to a single A-line, successive A-lines, or A-lines obtained from successive frames. One disadvantage of this technique is the intensive computational cost. For example, a digital Hilbert transform in Matlab involves fast Fourier transform (FFT) and inverse FFT operations. Several groups have implemented the algorithm with additional hardware, such as dual processors, a quadrature detector, or an optical Hilbert transformer to achieve real-time and multifunctional imaging.^{13, 14, 15, 16} Another typical frequency estimation method is the correlation technique, which is well established in ultrasound and can be implemented in hardware for fast processing.¹⁷ It has been successfully applied in OCT to estimate Doppler shifts in real time.^{18, 19} Both phase-resolved and correlation techniques are very sensitive to noise because of the use of adjacent phase changes.

In this work, the Doppler frequency estimation problem is treated from a different approach. Time-domain OCT heterodyned signals can be modeled as data series with a single time-variant central frequency, similar to some signal models in radar, sonar, and communication systems. In these systems, the frequency-tracking filters have been designed to estimate central frequencies based on an adaptive filter theory.²⁰ We introduce the use of an infinite impulse response (IIR) adaptive notch filter (ANF) to track the Doppler frequency.^{21, 22, 23, 24} This filter has a compact computation structure amenable to real-time processing. In the following, we introduce the mathematics of the filter and discuss its performance in a quantitative view by an A-line simulation model. Images of an intralipid conduit and an in-vivo skin sample are also presented. The results using the ANF are compared to those of several established algorithms. Both simulations and sample data demonstrate that the ANF is a very effective frequency estimator for OCT signals.

2.

Theory

A notch filter is a filter that contains one or more deep notches or nulls in its frequency response characteristic. It can remove certain frequency components from the input signal. A pair of complex-conjugate zeros on the unit circle at $z_{1,2}=\exp (\pm j{\omega }_0)$ can create a null at frequency ${\omega }_0$ with the system function given as:

Fig. 1

(a) The frequency response of the IIR notch filter [Eq. 2] with two sets of $b$ and $\alpha$ values. (b) The adaptive notch filter diagram. BPF stands for band-pass filter. Upper branch is the ANF for Doppler frequency estimation. The lower branch is a dynamic bandpass filter for OCT.

3.

Methods

Mathematically, the estimation precision of an estimator can be evaluated via the Cramer-Rao lower bound (CRLB). Yazdanfar ²⁵ have discussed the use of the CRLB for the frequency estimation problem in ODT, and pointed out that the estimated frequency precision could surpass the Fourier limit $f_{\min }=1∕t$ barrier and achieve the lower bound. As a well-developed method, the estimation accuracy of the ANF estimator has been systematically studied in Ref. 22 (stationary frequency case) and Ref. 23 (nonstationary frequency case) in detail using the CRLB. These studies showed that the ANF, as a prediction error method, is a nearly optimal estimator with a mean square error close to the CRLB. Besides the use of established mathematical approaches for estimator evaluation, we also implement a sinusoid-in-noise model to assess quantitatively the estimation performance of the ANF, such as estimation accuracy and processing speed The signal model is given as:

To verify the imaging quality using the ANF technique, ODT image data were acquired from a fiber-based OCT system developed earlier,¹⁰ where two galvanometers were synchronized for depth and lateral scan. The wavelength was $1.3\mu \mathrm{m}$ with an axial resolution of about $16\mu \mathrm{m}$ . The 0.5% intralipid solution was pumped into a circular conduit by a peristaltic pump and the pulsatile fluid flow was imaged. The glass conduit had a $1.0\text{-}\mathrm{mm}$ inner diameter and $1.2\text{-}\mathrm{mm}$ outer diameter. The carrier frequency of the flow images was $60\mathrm{kHz}$ . The incident angle was $80\mathrm{deg}$ and the maximum flow speed was about $50\mathrm{mm}∕\sec$ , corresponding to $14\text{-}\mathrm{kHz}$ maximum Doppler shift. The pulsatile flow data were processed with the seven algorithms. The cross sectional ODT images were also obtained from a female volunteer who had a spot of an aggregated blood vessel on the subepidermal area of her hand. The sample-arm grin lens was placed perpendicularly to the skin surface and laterally scanned.

4.

Results

With the simulation model described in Eqs. 7, 8, the estimated frequency profiles (dotted line) and actual frequency profiles (solid line) are plotted together for the seven algorithms using two SNR levels. In Fig. 2 , where the SNR is $5\mathrm{dB}$ , the centroid algorithm gives a noisy and underestimated frequency profile. The weighted centroid, adaptive centroid, and sliding-window filter algorithms emphasize the central frequency and result in improved estimation accuracy. The estimation of phase-resolved and correlation algorithms is very noisy at this SNR level. The ANF and sliding-window algorithms are the most noise robust and provide the accurate frequency estimation profile. The CRLB analysis we discussed gives a mathematical evaluation of the ANF performance. Besides this, we can use a simple picture to explain the noise robustness of this filter: during the adaptive process, the response of the ANF to additive noise is like an additional “random walk” around the tracked frequency signal. The amplitude of the random walk is regulated by the small step-size parameter $\mu$ and will not be very large in each adaptive step. In Fig. 3 , where SNR is $25\mathrm{dB}$ , all seven algorithms show consistent increase in estimation accuracy. A notable change at this low noise level is the improved performance of phase-resolved and correlation techniques. ANF again provides the smooth and unbiased estimate of the parabolic frequency profile. To give a quantitative comparison, the estimation errors of seven algorithms with different SNR levels are plotted in Fig. 4 . The root mean square error (RMSE) between the estimated and actual parabolic frequency of the 1024-point A-line is calculated to assess the overall estimation accuracy. As we can see, the ANF achieves the best performance among all algorithms at $\mathrm{SNR}<17\mathrm{dB}$ . The sliding-window filter algorithm also shows favorable estimation accuracy at low SNR. In fact, the two algorithms have many similarities, because both of them use a filter approach in estimating Doppler frequency. The sliding-window filter uses a narrow sliding bandpass filter, while ANF uses a narrow tracking notch filter. However, the sliding-window filter needs to search and compare the signal power with a number of filtering operations in a data window to obtain each frequency estimate. The selected number of filtering operations gives a tradeoff between estimation accuracy and computation cost. On the other hand, the ANF continuously and adaptively locates the signal frequency at each estimate with only one simple filtering process. Therefore, the ANF has an overall better performance than the sliding-window filter technique. Only the phase-resolved and correlation algorithms surpass the ANF in estimation accuracy at high SNR. A real OCT signal, however, does not always possess such high SNR, especially for a deeper target. The phase-resolved and correlation algorithms give also poor estimation accuracy at low SNR, which suggests that these algorithms are more sensitive to noise.

Fig. 2

Estimated laminar flow profile (dotted line) and actual flow profile (solid line) of seven algorithms with $\mathrm{SNR}=5\mathrm{dB}$ : (a) centroid, (b) weighted centroid, (c) adaptive centroid, (d) sliding-window filter, (e) phase resolved, (f) correlation algorithm, and (g) ANF.

Fig. 3

Estimated laminar flow profile (dotted line) and actual flow profile (solid line) of seven algorithms with $\mathrm{SNR}=20\mathrm{dB}$ : (a) centroid, (b) weighted centroid, (c) adaptive centroid, (d) sliding-window filter, (e) phase resolved, (f) correlation algorithm, and (g) ANF.

Fig. 4

The estimation error (RME) of seven algorithms versus signal-to-noise ratio. The ANF algorithm has the best estimation accuracy when SNR is below $17\mathrm{dB}$ .

The signal processing is performed on a Pentium 4 processor in Matlab. The computational cost for 100 A-lines (1024-point/A-line) for each algorithm is listed in Table 1 . The ANF has the fastest processing speed among all seven algorithms. This can be attributed to two reasons. First, the ANF is a time-domain filter with only real computations. Most other algorithms listed in the table involve more intensive computations due to the Fourier transform. The sliding-window filter method does not require Fourier transform; however, it is computationally intensive with a large number of filtering calculations. Second, this IIR notch filter is a minimal-parameter filter, i.e., only a bandwidth parameter $\alpha$ and a notch frequency parameter $b$ are used in the filter structure. When using a transposed (direct form II) IIR realization, the ANF takes only five real additions, five real multiplications, plus an $\mathrm{arcos}(-b∕2)$ operation for one frequency data, while the centroid algorithm (second fastest in the table) needs $32{\log }_{2}32$ complex multiplications and additions for FFT, followed by the centroid calculations. In fact, the computation time in Table 1 is still not the actual CPU consumption on these algorithms under equal conditions. The processing of FFT-related algorithms is optimized because of the built-in FFT function in Matlab, but the ANF does not need this advantage, as it requires only simple calculations. The computational simplicity and efficiency of the ANF algorithm will be more apparent in real-time ODT systems^{10, 15, 18} where all the algorithms are implemented in “C” or a digital signal processor (DSP) instead of in Matlab.

The images of intralipid flow experiments are shown in Figs. 5a, 5b, 5c, 5d, 5e, 5f, 5g , which correspond to the centroid, weighted centroid, adaptive centroid, sliding-window filter, phase-resolved, correlation, and ANF algorithms, respectively. The time interval between adjacent A-lines at the same depth was $1∕128\mathrm{s}$ with a spatial distance of about $10\mu \mathrm{m}$ . The maximum flow speed was about $50\mathrm{mm}∕\mathrm{s}$ . Because of the phase wrapping and weak cross correlation between adjacent A-lines, the phase-resolved and correlation techniques were processed on a single A-line only, and no cross-correlation was applied. The imaging target was diluted intralipid without strong backscattering. The average SNR was approximately a few decibels, which was close to the lower side in our simulation shown in Fig. 4. The vertical color strips appearing in the image are due to the coupling of the pulsatile flow with the lateral scan. All images were obtained with no additional threshold or averaging operation. The images from centroid and sliding-window filter algorithms [Figs. 5a, 5b, 5c, 5d] are relatively smoother than the phase-resolved [Fig. 5e] and correlation algorithms [Fig. 5f], which are more subjective to noise. Many noise spikes are seen in Figs. 5a, 5b, 5c, 5d, 5e, 5f, particularly in the low SNR region outside the conduit circle. In Fig. 5g, the ANF produces the smoothest flow image with a remarkable reduction of noise. The ANF also provides fine spatial resolution and the fastest processing speed. Besides these good properties, some meteorite-like tails can be observed in the image. This is a unique feature from the ANF technique. In fact, these tails occur only in the air region outside the circular tube, where no OCT signal exists except a small carrier frequency residual. In these regions (close to white noise), the adaptivity of ANF is reduced. The line-shape structure appears along the depth direction because the tracked frequency at the last signal region is “memorized.” Several methods can be used to avoid this effect. For example, a noise threshold can remove these tails. In Fig. 6 the structural and blood flow images are obtained from the subepidermal area of a human hand. The flow area is highlighted with a color map obtained from ANF and superimposed on the grayscale structural image in Fig. 6b, where two vessel cluster spots with different flow directions can be identified. For comparison, the color map obtained from FFT-based weighted centroid is also given in Fig. 6c. Frequency and power thresholds were applied to select the flow region. This in-vivo ODT flow image demonstrates the capability of the ANF technique to estimate subsurface biological blood flows.

Fig. 5

The ODT images by centroid, weighted centroid, adaptive centroid, sliding-window filter, phase-resolved, correlation, and ANF techniques from (a) to (g), respectively. The image size is about $1.5\mathrm{mm}$ in width and $2.4\mathrm{mm}$ in depth.

Fig. 6

In-vivo B-scan image of subsurface blood flow under the skin. The blue and red spots in the right image indicate two vessel clusters with blood flow in opposite directions. From left to right: the intensity image, overlapped with ANF, and overlapped with weighted centroid.

5.

Discussion and Summary

The ANF is an IIR filter with a constrained structure: its poles are always at a given radius $\left(\alpha \right)$ from the origin. Since in most implementations $\alpha$ is not an adaptive quantity, this means that stability of the filter is not an issue. However, ANF is a tracking algorithm, and as always with tracking algorithms there is a tradeoff between the ability to track quickly and track accurately. The relevant tracking quantity is the bandwidth of the tracking filter, and in the case of the ANF, this bandwidth corresponds nicely to the bandwidth of the filter itself. That is, when $\alpha$ is close to unity, the tracking accuracy is very fine, but the ability to track quick changes is less; conversely, as $\alpha$ decreases to some more moderate value, the response improves, but there is more estimation bias. It is a matter of tuning. That said, there do exist methods from the adaptive signal processing literature that can mitigate the issue: speed is often enhanced by an update that involves the sign of the error rather than the error itself (this is called the “normalized” least mean square); and of course there are schemes that involve a time-varying $\alpha$ . We have not explored these here. In reported simulations and experiments, no significant difference in filter performance was found when the $\alpha$ was varied between 0.75 and 0.95.

The step size $\mu$ is the most important parameter; however, the selection of the best value is somewhat ad-hoc. The optimal step size is proportional to the signal intensity and related to the velocity gradients. It is also affected by SNR. In general, the magnitude of $\mu$ imposes a tradeoff between the spatial and frequency resolution. When the step size is too small, the filter may not respond fast enough to adapt to a changing frequency. When the step size is too large, noise could be introduced. By further increasing the step size beyond a certain point, the ANF may fail to converge. Although some methods such as low and high frequency limiters can be used to prevent this, it has been demonstrated that stability is generally not a problem for ANF, provided that $\mu$ is sufficiently small.^{22, 24} For the steepest descent method given in Eq. 5, the actual adaptive frequency step $\Delta b$ is determined by the interplay of both the step size $\mu$ and the gradient of the cost function. When the notch center is not tuned to the desired frequency but is sufficiently close to it, such that the signal is appreciably attenuated, the convergence of ANF is proven to be exponential.²⁴ The Doppler frequency of an ODT signal in general varies smoothly and gradually. It is unlikely that the frequency change between adjacent data points will cause the notch filter to be tuned away from the signal bandwidth. Therefore, the tracking process of ANF for ODT signals is effective. In both simulation and experimentation, $\mu$ works well in a range of values without encountering stability problems. This can be seen from Fig. 7 , where the root mean square errors with different values of $\mu$ are presented. Since a small $\mu$ reduces the process error, we can consider using smaller $\mu$ with a few iterations to improve the accuracy but also keep the adaptive power. This improvement is, however, limited by the noise level of the signal.

Fig. 7

The estimation error (RME) versus the step-size parameter $\mu$ of the ANF algorithm in simulation with no additive noise. The error increases at both sides of the curve. At the left side, the ANF cannot catch up with the changing Doppler frequency because $\mu$ is too small. At the right side, the estimation error increases suddenly because ANF starts to oscillate when $\mu$ is too large. $\mu$ works well in the middle range of the curve.

Table 1

Process time (seconds) of seven algorithms using Matlab with a Pentium 4 desktop computer for 100 A-lines with 1024 data points.

Among the algorithms we discussed, the centroid method is perhaps the most simple and effective, with relatively fast processing speed. The phase-resolved and correlation techniques are more sensitive to noise. But they are capable of achieving the highest flow sensitivity when calculation is performed on sequential A-lines. The sliding-window filter technique has been found to have good estimation accuracy and noise robustness. This technique does, however, require intensive computation. Both the sliding-window filter and the ANF technique are filter methods, but the ANF has a much more efficient structure. Regarding frequency resolution or velocity sensitivity, the ANF is quite good and it can be and has been compared to the phase locked loop. FFT-based methods are limited to the fineness of the frequency-domain sampling unless interpolation is used. A greater issue is the effect of noise: interpolated FFT methods can approach the CRLB for accuracy, but the performance of the ANF is not well understood. Nevertheless, we have found that an ANF works reliably down to $0\text{-}\mathrm{dB}$ SNR for sinusoid signals. Note that the ANF is applicable to the time-domain OCT, but may not be suitable for Fourier domain OCT.

In conclusion, we present a novel and effective frequency estimation technique suitable for ODT applications. Simulation and experimental results demonstrate several excellent properties of ANF such as high estimation accuracy, robustness to noise, and fast signal processing speed.