A Time-Frequency Approach for Cerebral Embolic Load Monitoring

Objective: To enable reliable cerebral embolic load monitoring from high-intensity transient signals (HITS) recorded with single-channel transcranial Doppler (TCD) ultrasound. Methods: We propose a HITS detection and characterization method using a weighted-frequency Fourier linear combiner that estimates baseline Doppler signal power. An adaptive threshold is determined by examining the Doppler signal power variance about the baseline estimate, and HITS are extracted if their Doppler power exceeds this threshold. As signatures from multiple emboli may be superimposed, we analyze the detected HITS in the time-frequency (TF) domain to segment the signals into individual emboli. A logistic regression classification approach is employed to classify HITS into emboli or artifacts. Data were collected using a commercial TCD device with emboli-detection capabilities from 12 children undergoing mechanical circulatory support or cardiac catheterization. A subset of 696 HITS were reviewed, annotated, and split into training and testing sets for developing and evaluating the HITS classification algorithm. Results: The classifier yielded 98% and 96% sensitivity for 100% specificity on training and testing data, respectively. The TF approach decomposed 38% of candidate embolic signals into two or more embolic events that ultimately account for 69% of the overall embolic counts. Our processing pipeline resulted in highly accurate emboli identification and produced emboli counts that were lower (by a median of 64%) compared to the commercial ultrasound system's estimates. Significance: Using only single-channel, single-frequency Doppler ultrasound, the proposed method enables sensitive detection and segmentation of embolic signatures. Our approach paves the way toward accurate real-time cerebral emboli monitoring.


I. INTRODUCTION
A CUTE neurological complications remain an important clinical problem in patients undergoing extracorporeal membrane oxygenation (ECMO) [1]- [3] and ventricular assist device (VAD) support [4]. One cause of acute brain injury in these populations is cerebral embolism, which may be detected clinically in real-time by transcranial Doppler (TCD) ultrasonography as high intensity transient signals (HITS) within the Doppler spectrum [5]- [7]. HITS, representing cerebral emboli, may be composed of air, thrombi, atheromatous plaque, lipid, or platelet aggregates. Cerebral emboli can occlude the cerebral vasculature, potentially causing transient ischemic attacks, stroke, or other acute neurologic injury. A clear understanding of the prevalence and clinical significance of HITS in patients on mechanical circulatory support (ECMO, VAD) or undergoing cardiac catheterization, and at high risk of cerebral embolic events is lacking. In a previous study in children with congenital heart disease undergoing cardiac catheterization, we found the process of visual review and manual annotation of HITS and their classification into emboli and artifacts to be prohibitively time consuming and essentially impossible when HITS occurred in clusters (often designated as curtains or showers) [8]. We also found that commercial TCD emboli-detection software generated excessive false positive events.
Typical ultrasound-based emboli detection methods use baseband (Doppler) ultrasound signals from one or two depths and one or two simultaneous insonation frequencies [9]- [11]. The signals may first be prefiltered, for example using wavelet transforms [12]- [15], to help differentiate embolic signals from artifacts and background blood signatures. HITS may then be detected using the embolus-to-blood ratio (EBR), defined as the ratio of backscattered power from an embolic source, normalized by the power calculated over data segments not containing any emboli. Embolic sources tend to have a high EBR because of their size and acoustic impedance mismatch relative to surrounding blood [16], [17].
Robust emboli detection using EBR is a challenging task, however. A baseline Doppler power level of the normal (nonembolic) blood flow must first be established. This baseline This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ power estimate should vary with the cardiac cycle since the backscattered Doppler power due to pulsatile blood flow is modulated by an order of magnitude between systole and diastole. A dynamic detection threshold must also be determined, so that EBR excursions above that threshold can be flagged as candidate emboli. A subsequent artifact rejection stage is required since tissue or ultrasound probe motion can generate large excursions in EBR that should not be counted as emboli. Finally, multiple emboli may flow through the ultrasound sample volume simultaneously, and the corresponding Doppler signal should be decomposed into individual embolic signatures for accurate counting. Lipperts et al. [18], for instance, reported that existing commercial TCD systems do not accurately estimate the number of cerebral emboli in such situations. To our knowledge, the problem of automatically separating signatures from multiple simultaneous emboli using single-depth, single-frequency TCD systems has not been addressed in the literature.
In this paper, we describe a signal processing pipeline that enables real-time HITS detection and classification into likely emboli and artifact using single-channel, single-frequency Doppler data. We model Doppler baseline power as a Fourier series, and propose a weighted-frequency Fourier linear combiner (WFLC) [19] to adaptively estimate the Fourier coefficients in real-time. Variance of the Doppler power about this baseline leads to an adaptive HITS detection threshold. Disabling WFLC adaptation during HITS allows us to retain estimates of the signal background during prolonged periods of HITS showers or artifact. We then propose an algorithmic separation of detected HITS into signatures from individual emboli by time-frequency (TF) analysis. Finally, logistic regression classification is used to reject artifacts. The method was evaluated on data from twelve pediatric patients undergoing ECMO, VAD support, or cardiac catheterization.
We first outline the data collection and annotation steps in Section II. We then describe our emboli detection and TF-based separation approach in Section III. The artifact rejection classifier is described in Section IV, and we present and discuss the results of applying our processing pipeline in Sections V and VI, respectively.

II. DATA COLLECTION AND ANNOTATION
The study was approved by the Boston Children's Hospital Institutional Review Board. Written informed consent was obtained for all subjects from the legally authorized representative, and patient assent was obtained whenever possible. Children on mechanical circulatory support (MCS), i.e. ECMO or VAD, or undergoing cardiac catherization were eligible for study inclusion. Subjects who lacked an acoustic window to permit TCD ultrasound examination of the middle cerebral artery (MCA) were excluded after enrollment. Subjects underwent emboli monitoring of the right or left MCA with a dual frequency (2 + 2.5 MHz), range-gated, pulsed-wave TCD system (DWL Doppler-BoxX, Compumedics Germany GmbH, Singen, Germany). The ultrasound probe was handheld, or secured in place with a soft elastic headband, over the right or left temporal window. Emboli monitoring began once an adequate Doppler signal was obtained from the M1 segment of the MCA at the level of the bifurcation of the MCA and anterior cerebral artery. Data were collected from eight patients on MCS (3F, 5M, ages: 3 weeks to 14 years), and four patients undergoing cardiac catheterization (1F, 3M, ages: 4 months to 14 years). Recording durations ranged from 9 to 118 minutes for a total 625 minutes (10.5 hours) of data. Further clinical details of our patient cohort are provided in the appendix.
The comparatively large volume of ultrasound data collected precluded exhaustive manual HITS annotation and classification into embolic and artifact events. We therefore first extracted candidate HITS using an automated approach reported previously [20], and two expert annotators (KLL, BDK) were presented with candidate HITS so identified from a subset of seven MCS patients. Each annotator independently assessed each candidate HITS using previously published criteria for emboli detection [21] and indicated whether each identified segment was judged to be an embolic event, an artifact, or the annotator was unsure which of the two categories to assign. Only HITS marked by both annotators as either emboli or artifacts were used for training and testing. A 60% cohort of annotated data segments was randomly selected and used for training of the artifact-rejection classifiers (Section IV). The remaining annotated data from the MCS patients were used for testing classifier performance. To determine the robustness of our emboli classification approach, we retained the data from the four cardiac catheterization patients as an independent hold-out validation cohort that was neither used for classifier training nor testing.

A. Data Preprocessing
The DWL Doppler-BoxX exports Doppler data in binary format along with timestamps of the emboli detected by the device's proprietary software. The device exports the inphase, r i , and quadrature, r q , demodulated signals for the selected target depth from one insonation frequency (2 MHz) [22]. From the exported signals we form the complex signal r n = r i n + jr q n , where n is a discrete sampling index with samples recorded at the pulse repetition frequency, PRF. Since the DWL system generates separate binary files whenever the acquisition parameters are modified during a recording session, we concatenated the Doppler streams from each file by rescaling the signals to a common signal gain and by using MATLAB's resample function to resample all segments to the highest PRF used during the recording session. In accordance with prior work [9], [10], we computed the signal power, P, in non-overlapping data windows of 2 ms duration. For the m th non-overlapping window of length N p , the power was computed as Since the Doppler signal power can vary by an order of magnitude during the cardiac cycle, we base our HITS detection and segmentation approaches on the log-transformed power P m = 10 log 10 (P m ).

B. HITS Detection
For HITS detection, we propose an adaptive baseline power estimation approach that uses a modified WFLC [19]. The WFLC was originally developed for canceling physiological tremor in robotic surgery applications; it models a quasi-periodic signal as a Fourier series, estimating the Fourier series coefficients and the harmonic frequency in real-time and in an adaptive manner. The original WFLC method was designed to continually update its parameters. In our approach, we update the parameters only during baseline flow conditions and forgo updating when a candidate HITS is detected. For each 2 ms data window, we compute the difference, e m , between the logtransformed power estimate, P m , from a predicted background power, P m , for that window. A HITS is detected if e m > γ m , where γ m is an adaptive threshold. The WFLC parameters and γ m are retained (i.e. not updated) if a HITS is detected and updated otherwise. The resulting algorithm architecture is illustrated in Fig. 1.
More concretely, given initial estimates for the filter weights, w 1 and w 0,1 , and the fundamental frequency, ζ 0,1 , a prediction at a later sample is computed as To update the parameters from one window to the next, we define the prediction error where the latter condition occurs during a HITS. Setting the associated error term to zero prevents parameter adaptation to the embolic or artifact signal properties. The WFLC parameters are then updated by performing a gradient-descent step in which where μ, μ 0 , and μ ζ are preset adaptation parameters [19]. To initialize the computation, we provide estimates of the fundamental frequency (heart rate), and the corresponding Fourier series coefficients by computing the discrete Fourier transform of the first 10 seconds of P m and by analyzing the dominant frequencies, amplitudes, and phases.
To determine the detection threshold, we examine the standard deviation of prediction errors. (A similar strategy was previously proposed in [23].) In our method, the detection threshold is set to where α is a tunable parameter and γ m is not allowed to exceed 15 dB for system stability. We empirically set α = 3 to strike a balance between high probability of detecting emboli and acceptable probability of false alarm. SD(e m ) is a recursively low-pass filtered version of the standard deviation of the prediction errors, SD(e m ), in HITS-free segments (7) where α lp was set to 0.9. The WFLC parameters used in our analysis are summarized in Table I, and the resulting baseline and threshold estimates for a representative data segment are shown in Fig. 2.
The WFLC parameters are reinitialized if a HITS segment longer than T r einit = 10 s is detected. This is to prevent changes in signal quality or probe position from being falsely detected as embolic signatures. The detected candidate HITS segments can be quite long in duration, which can lead to significant computation load in the subsequent TF analysis. To reduce this load, we split HITS into sub-segments of at most T max = 0.25 s.

C. HITS Separation
HITS detected by the WFLC-based method can often appear as consisting of multiple individual embolic signatures that are temporally merged. We therefore further examined the fine structure of the detected HITS and separated HITS into constituent embolic events via TF analysis.
To perform the TF analysis, we used a discrete-time approximation of the continuous wavelet transform. Specifically, a detected HITS is passed through a filter bank that enables realtime computation. Each filter is a Gaussian kernel modulated to a center frequency, f v . Each center frequency corresponds to a Doppler velocity, v, according to the Doppler equation f v = 2f 0 v/c, where we assume c = 1540 m/s as the speed of sound, and f 0 = 2 MHz is the transmitting frequency. We denote the resulting TF decomposition at time sample n and Doppler velocity v as R n,v where the convolution is performed in time and h n,v is a bandpass filter of the form The filter has a temporal spread governed by σ v , and h 0 is a scaling constant. In our approach, we select the center veloci- m/s, and 200 center velocities are used. We set σ v = SV/βv, where β = 10 is a scaling constant, and SV is the value of the sample volume selected during data acquisition. (The term sample volume is a misnomer since it represents the axial length of the insonated region and not a volume; we retain its use since the term is widely accepted.) Filters for higher center velocities therefore have narrower temporal spread, allowing finer temporal localization of embolic signals. A given TF image may then be inverted back to the time domain as We first segment HITS in the TF domain before conducting a linkage step to merge signatures that may correspond to the same embolus. The resulting merged signatures are then inverted back to the time domain as illustrated in Fig. 3. The segmentation and merging steps proceed as follows: 1) TF Segmentation: For each TF domain image, we generate a threshold and segment the absolute value of the TF image of the selected HITS. The threshold is generated by applying Otsu's method [24] on log-compressed absolute values of the TF representation (MATLAB's graythresh function), and taking the anti-log of the resulting threshold. Log-compression is used since the TF pixel values can vary by several orders of magnitude. Applying the thresholding method on the raw TF images may therefore lead to unsuitably high thresholds. Regions of the absolute TF representation that are higher than the threshold are segmented into patches. First, a rescaled TF image, RS n,v is generated according to where R min and R max are the minimum and maximum absolute values in the TF representation, respectively. Rescaling allows the application of the H-minima transform [25] to RS n,v in order to suppress local minima; we used an empirically determined suppression threshold of 0.001. The Watershed image segmentation algorithm [26] is then used on the resulting image to extract patches that are above the detection threshold. For each patch, we compute the location of the highest intensity, (n max , v max ), and the normalized traveled distance, ND. The latter is computed by first determining the instantaneous velocity IV n for each sample, n, and subsequently integrating the velocity. The absolute of the resulting displacement is normalized by the sample volume, SV. The instantaneous velocity, IV n , is estimated by computing a weighted average of the TF image for each n, such that IV n = v |vR n,v |/ v |R n,v |. The metrics n max , v max , and ND are used subsequently to merge patches that potentially correspond to the same embolus.
2) TF Merging: The segmentation process may result in separate patches that belong to the same embolic signal. To avoid such spurious fragmentation and overcounting of embolic events, a merging step is necessary. We designed a set of rules to determine if such merging is necessary. Patches are merged if they are close in speed and time, have not individually traversed a sizable fraction of the sample volume, and do not lead to large traveled distances when combined together. Specifically, two patches i and j are merged on the basis of 1) the time between their intensity maxima (|n max,i − n max,j |/PRF < T min ), 2) the absolute difference between their velocity maxima, the normalized displacement of the union of the patches, (ND < ND max ). Here, T min , ND min , ND max , and Δv are predefined thresholds set to 6 ms, 0.85, 1.25, and 0.5 m/s, respectively, and all conditions must be met for a merger. In our approach, we consider all possible pairs of HITS until we merge a pair that fits these criteria. The process is then repeated for the new set of patches until no further matches can be made. The algorithm reverts to the segments fed originally to the TF-based segmentation stage if the merging process does not converge within a maximum number of passes, set to 100. Finally, the merged segments are converted to the time domain. Artifacts in the remaining segments are then removed using a feature-based classifier described below.

IV. ARTIFACT REJECTION
Since embolic signals are generally longer than 8 to 10 ms in duration [8], [21], we first rejected any detected HITS from further analysis if their duration was less than 6 ms (or three 2 ms data windows). To classify the remaining HITS into likely emboli or artifact, we applied a feature-based logistic regression classifier. We computed and evaluated six candidate HITS features and selected a subset of three features for our final artifact rejection classifier. The final classifier is applied twice, once after the initial WFLC-based HITS detection step and then after the final TF emboli separation step.

1) Unidirectionality:
Emboli are known to move in the direction of blood flow [21] leading to a single-sided Doppler frequency spectrum. A quantitative measure of such unidirectional flow is the ratio P ≥0 /P < 0 [9], where P ≥0 and P < 0 are the power of the HITS in the positive and negative frequency bands, respectively, and blood flow is assumed to be in the positive direction. It is possible, however, to simultaneously insonate two vessels with opposite flow directions, and thus a dominant blood flow direction cannot be assumed a priori. Also, the ratio can assume arbitrarily large values. Thus, we define the (nondimensional) unidirectionality, U , as where u = max P ≥0 P < 0 , where we set u max = 1000 and computed P ≥0 and P < 0 by summing the squared magnitude values of the Blackman-windowed discrete Fourier transform, F (ω), computed over the duration of each candidate HITS (Fig. 4).

2) Spectral Concentration:
We expected emboli to travel at a finite range of velocities, leading to frequency spectra concentrated around a center frequency. We computed a measure of such spectral concentration as max ω |F (ω)| ω |F (ω)| with values close to unity indicating a high degree of spectral concentration and values close to zero indicating a broad frequency spectrum.
3) Speed: In contrast to emboli, artifacts commonly have bidirectional frequency spectra [9], [21]. Thus, we expected artifacts to have average Doppler speeds close to zero, and computed the instantaneous signal frequency, IF n , by numerically differentiating the unwrapped instantaneous phase [27], arg{r[n]}, of the Doppler signal According to the Doppler equation [28], the instantaneous velocity, IV n , is then where we assume that the insonation direction is parallel to the flow direction. We then define the HITS speed for the i th HITS as where the time index n spans the duration of the detected HITS.
(Here, the definition of instantaneous velocity introduced earlier in Section III-C was not used in order to bypass the need to convert HITS into their TF representations.) 4) Normalized Distance: Motivated by the work of Smith et al. [29], we note that emboli tend to traverse a significant fraction of the target SV. Thus, we integrate IV n over time, and normalize the absolute value of the resulting displacement by SV. In our implementation we use trapezoidal integration to compute the HITS displacement before normalizing the absolute value of the result by SV.

5) Temporal Skewness:
We observed that artifacts tend to have a significantly skewed temporal envelope (Fig. 4). We therefore defined skewness as the time from the start to the peak of the envelope divided by the total HITS duration. A value of 0.5 indicates no skew; artifacts tend to have a temporal skewness value that is small compared to this reference.

6) Measured/Expected Duration:
In our visual review of sample data, we found artifacts to have a short duration compared to embolic signatures that are expected to have durations corresponding to their speed and SV [29]. Thus the ratio of measured to expected duration may provide a means of separating artifacts from embolic events. We computed the expected duration as d i = SV/s i .

B. Classifier Design
We employed logistic regression in our artifact rejection classifier. Emboli were assigned the value of 1 and artifacts the value −1. Classifiers were trained on emboli and artifacts; HITS classified as unsure were excluded from our analysis. For the i th HITS, the classification function is of the form where y i is the algorithm-assigned label, g i = [1, g i1 , ..., g iJ ] is the vector of J features augmented by a bias term, h = [h 0 , h 1 , ..., h J ] is the vector of model parameters, and η is the classification threshold. Features were first converted to Zscores by subtracting the respective feature mean, and dividing by the feature standard deviation in the training data. The model parameters were obtained by minimizing the l 2 -regularized logistic loss function where y i is the training label assigned to the i th HITS, I is the number of training samples, and the regularization parameter λ was empirically set to 1.

C. Classifier Evaluation
We analyzed feature statistics and trained logistic regression classifiers on the individual features to assess their artifactrejection performance. We then selected the three top performing features in this univariate analysis for inclusion in the final classifier. The final three-feature classifier was applied after WFLC-based HITS detection and again after the final TF-based emboli separation step.
We evaluated classifier performance by computing classification sensitivity and specificity. By varying the classification threshold we obtained the full receiver operating characteristic (ROC) curve for each individual classifier, and for the final three-feature classifier. To select the detection threshold for each classifier, we computed the distance from each point on the ROC to the (0,1) point on the ROC plot, thereby giving equal weight to both sensitivity and specificity. We selected the threshold value corresponding to the point on the ROC that minimized that distance.
A randomly selected 60% subset of the agreed-emboli and agreed-artifacts HITS annotations from the seven MCS patients was used for classifier training and threshold selection, and the remaining 40% were used for testing classifier performance. To ensure robustness of our approach, we applied the classification rule to annotated HITS from an independent hold-out validation data set, consisting of 500 emboli and 133 artifact annotations from the four patients in our study cohort undergoing cardiac catheterization.

A. Data Annotation and Inter-Rater Variability
Each annotator reviewed and scored 696 detected HITS events from seven MCS patients. Per-patient annotation counts ranged from 50 to 200. Notably, all annotated emboli events came from just two patients. The annotation results are summarized as a confusion matrix in Table II, and Cohen's kappa metric [30] for inter-rater agreement was 72%. The annotation accuracy, or fraction of annotations on the main diagonal of the confusion matrix, was 83%. We trained and tested our classifiers on the 482 annotated HITS events of agreed-embolic and agreed-artifact events.

B. Artifact Rejection
Box-plots of the feature values obtained on the training data for each of the six candidate features are shown in Fig. 5. Good separation of the median feature values for emboli and artifacts were achieved for unidirectionality, duration ratio, and normalized distance. The single-feature classification performance is summarized in Table III. Based on their classification performance, we selected unidirectionality, duration ratio, and normalized distance for inclusion in the three-feature Fig. 6. Scatter plot and 2D projections for three features of emboli (red) and artifacts (blue) for the validation data. Optimal artifact rejection decision boundary for the classifier is shown in magenta and was based on the optimal decision thresholds derived from the training data.
The trained classifier assigned weights of 2.5, 1.3, and 1.2 to Z-score-normalized unidirectionality, duration ratio, and normalized distance, indicating that the classification is driven strongly by the unidirectionality feature. Emboli and artifacts from the validation cohort are shown in the scatter plot of Fig. 6 along with the classifier boundary derived from the training data. The projections in Fig. 6 demonstrate that duration ratio and normalized displacement show strong collinearity for emboli and may only add incrementally over the classification performance of the unidirectionality feature. When we performed a formal sequential feature selection approach using MATLAB's sequentialfs function, only the unidirectionality and normalized distance features were selected. The two-feature sensitivity and specificity were 99.0% (95% CI 97.1-100.0) and 99.5% (95% CI 98.5-100.0) for the training set, 100.0% (95% CI 100.0-100.0) and 98.4% (95% CI 96.1-100.0) for the testing set, and 96.4% (95% CI 94.8-98.0) and 99.2% (95% CI 97.8-100.0) for the validation set. The slight loss in sensitivities on the validation data set is to be expected given that the algorithm was not trained on any of the validation data and the fact that the validation data were captured from a different clinical intervention from the ones represented in the training data.

C. Patient Embolic Loads
We applied the final emboli detection pipeline-consisting of WFLC-based adaptive HITS detection, TF emboli separation, and artifact rejection-to the entirety of all twelve patient recordings. Of the WLFC-derived embolic HITS, 38% were further segmented into two or more embolic events by the TF-based HITS separation approach. This decomposition accounted for 69% of the final emboli count (Fig. 7), emphasizing the need to incorporate an emboli segmentation step into EBR-based HITS detection approaches.
Representative cumulative embolic counts for three recordings are shown in Fig. 8, and the embolic loads for all records are summarized in Table IV. The table also lists embolic counts as derived by applying the three-feature artifact-rejection classifier to the WFLC-derived HITS (i.e. without applying the  TF HITS separation). Relative to the manufacturer-provided counts, we observed a median percentage reduction of 64% in emboli counts.

VI. DISCUSSION
Accurate, real-time emboli monitoring remains an open problem in the pediatric population. In adults, real-time emboli monitoring during carotid endarterectomy can alert the surgeon to incorporate cerebral protection measures [31]. In cardiopulmonary bypass, it led to the change from bubble to membrane oxygenators and the introduction of arterial line filters [32]. Approximately 10% of neonates and infants have seizures (clinical or subclinical) following heart surgery [33]. As seizures have been associated with adverse neurodevelopmental outcome [34], correlating the burden of emboli with post-operative seizures may lead to new strategies for their prevention.
Several limitations of existing Doppler-based embolus detection methods have been reported in the literature. These include requiring computations that operate over large signal blocks, thereby limiting real-time operation [11], generation of excessive false positive events [8], and an inability to distinguish multiple emboli that flow through the insonation region simultaneously [18]. We have developed a novel single-depth, singleinsonation-frequency embolus detection method that attempts to address these problems.
We introduced a WFLC framework to generate baseline power estimates of received Doppler data. Segments whose power exceeds an adaptively estimated threshold were selected as candidate emboli. We integrated a time-frequency segmentation step into our algorithm that attempts to separate signatures from emboli that flow into the ultrasound beam concurrently. When compared to the embolus detection performance of a commercially available two-depth, dual-frequency device, our method led to a median reduction in embolic counts by 64% in a pediatric patient cohort.
Computation requirements: Our system does not utilize information from future signal values, thereby allowing it to function in real-time, albeit with latency inherent in the internal computations. Preliminary HITS detection is performed with minimal delay since signal power computation introduces only a 2 ms latency (as 2 ms nonoverlapping data windows are used), and because the WFLC algorithm does not introduce additional delay-it was designed for zero-phase cancellation of periodic disturbances [19]. An artifact-rejection classifier is applied to minimize subsequent computation burden. The classification procedure itself uses three easy-to-compute features in a simple logistic regression model. We designed our finite impulse response filter bank such that each filter has the same group delay [27]. Thus, these filters may be implemented as a set of parallel, causal delay lines to generate time-frequency representations of candidate HITS with a latency equal to the group delay. In our subsequent TF analysis, we employ commonly used image processing tools, optimized implementations of which are readily available for target deployment platforms.
Artifact rejection performance: Our classification performance is predicated on HITS training labels provided by our expert annotators, who achieved an inter-rater reliability (Cohen's kappa) of 72%. Our reference annotations may thus be interpreted as reliable [35]. Similar kappa values (72%, and 90%) have previously been reported in the literature for embolus annotations by human experts [36], [37].
Our artifact rejection scheme uses a logistic regression classifier that allows interpretation of the factors driving high classification sensitivity and specificity. Specifically, upon examining the classifier weights, it may be seen that there is greater emphasis on unidirectionality. The attained high classification sensitivity and specificity are on par with those reported in prior literature [9], [13], [38], [39]. For instance, Darbellay et al. [13] reported embolus classification sensitivity of 95% and associated specificity of 97% on a testing data set comprising 600 emboli and 530 artifacts. Using seven classification features, Sombune et al. [39] recently reported an average classification sensitivity of 91.5%, average specificity of 90.0%, and average accuracy of 90.5%, outperforming the work by Karahoca and Tunga [38]. Brucher and Russel [9] previously proposed using four features in a decision tree: difference in Doppler shift due to dual-frequency insonation (2 and 2.5 MHz), a measure of expected signal duration, emboli presence in a second depth, and unidirectionality. They reported that 99% of all artifacts and emboli were classified correctly by their system in a data set comprising 554 emboli and 800 artifacts. In our approach, we found that while the HITS speed (or equivalently, the signal frequency) is different between artifacts and emboli (Fig. 5), the attained classification accuracy in our data set for this feature was not strong. Also, our classification approach only uses information from one depth and one insonation frequency. During our initial experiments, we found that several emboli may flow simultaneously, making it difficult to reliably match their signatures across different depths. The traveled distance feature has previously been shown to discriminate between gaseous and solid emboli [29]. We found this metric to be useful in separating artifacts from emboli as well.
Embolus separation using TF analysis: Multiple emboli can often be generated simultaneously, such as during catheter manipulation during cardiac catheterization or aortic cross-clamp release in cardiac surgery [8], [18]. Single-channel Doppler devices have been reported to be incapable of reliably detecting emboli in such circumstances [18]. Instead, methods have been proposed that use information from multiple depths (M-mode imaging) [40] or raw radio-frequency (RF) data [18]. Lipperts et al. [18], for instance, proposed an image processing approach using successively received RF ultrasound signals to improve the estimation of the number of emboli encountered in embolic showers during cardiac surgical procedures. They claim that existing TCD systems do not accurately estimate the number of cerebral emboli during such showers. Using RF data is akin to processing information from a range of depths, and allows the authors to separate signals from multiple emboli more easily, albeit at the expense of processing requirements.
To our knowledge, the problem of automatically separating signatures from multiple simultaneous emboli using single-depth, single-frequency TCD systems has not been addressed in the literature. Colantonio and Salvetti [41] extracted HITS patches from Doppler TF images using a line-tracking procedure, but did not explicitly attempt to separate close HITS. Moreover, in their approach, the authors use a segmentation threshold determined via a pre-trained neural network. Likewise, in [11], the authors extract a region of interest in HITS spectrograms by examining asymmetric (unidirectional) flow regions, without attempting to separate individual HITS. In our study, 38% of HITS detected by the WFLC stage were subsequently split into two or more emboli by the TF processing stage. Emboli split in this fashion accounted for 69% of the total embolic load, suggesting the potential need to incorporate such emboli segmentation into emboli detection systems.
We believe that extracting individual emboli signatures is important, not just to establish accurate emboli statistics, but also for subsequent characterization of embolic signal properties (for example, their material composition). At present, we employ a set of simple heuristic rules that determine how TF patches are merged by analyzing the net traversed distance of the patches and the difference in velocities of those patches. In doing so, we implicitly assume that the underlying emboli do not have a wide size range (by constraining normalized displacements between ND min and ND max ). It has been reported that in adults, particulate emboli with diameters below 100 μm are unlikely to be detected via Doppler ultrasound owing to the diameter of the MCA [13], [16]. Likewise, particulate emboli with sizes above 240 μm are thought to cause stroke [13]. None of the patients in our cohort suffered from a clinically apparent stroke, and hence it is plausible that the particulate emboli in our data were within a narrow size range. Future work, however, can focus on assessing particle size to further guide the TF patch merging process.
Comparison with the DWL: The DWL device exports its detected HITS with a timestep granularity of 10 ms, thereby preventing a segment-by-segment comparison between embolic counts. We found, however, that our embolic counts exhibit greater sensitivity during embolic showers, as exhibited by larger steps in the cumulative counts in Fig. 8. At the same time, we found that in several recordings, the device's cumulative counts exhibit linear trends, suggesting a constant background embolic rate. Our method does not show such linear trends, and this difference could be due to both different detection sensitivities and embolic classification steps. On the whole, our method reduced the embolic counts by a median 64%, potentially suggesting that the device may be generating excessive false positive events.
Contributions: We have proposed a single-depth, singlefrequency Doppler based approach to detect, classify, and separate closely opposed emboli. The initial detection is performed via a WFLC-based method. This is attractive because it enables modeling the pulsatile nature of blood flow and also computes an adaptive detection threshold in real-time using simple computations for high detection sensitivity in both systolic and diastolic segments. We integrated our simple and interpretable logistic regression based artifact-rejection scheme into a TF processing approach in order to separate HITS into individual embolic events that may overlap in both time and frequency (velocity) using a single Doppler channel. The proposed approach, when applied to data from pediatric patients ranging in age from 3 weeks to 14 years, reduced the median embolic counts by more  than a factor of two, thereby warranting further exploration of accuracies of commercial devices. Future work should also focus on examining differences between embolic signatures of gaseous and particulate emboli. Likewise, integrating the ability to size emboli will enable better separation of HITS that occur simultaneously.
Limitations of current work: Currently, our method's performance has been evaluated on a small data set in which ambiguous HITS were excluded and ground truth information about the type, number, and size of emboli was missing. Further work is needed to test the classifier on more heterogeneous test sets, potentially in flow phantoms where embolic composition and size can be controlled and analyzed (or more reliably determined). Likewise, our TF method will need to be tested in a variety of flow environments on a range of embolic sizes and compared against ground truth data in order to further assess its detection ability.

VII. CONCLUSION
Patients with a variety of clinical conditions are susceptible to embolic events and stroke. Single-channel Doppler devices are commonly used to detect emboli, but current commercial TCD systems seem to overestimate embolic load. Our proposed embolus detection approach advances single-channel Doppler emboli monitoring by: 1) introducing a novel emboli-detection algorithm, coupled with artifact rejection stages that use simpleto-compute features; and 2) by separating embolic signatures through time-frequency processing. Our method paves the way for more reliable embolic load assessment so that appropriate clinical trials can be designed that may lead to improved patient care and neurologic outcomes.

ACKNOWLEDGMENT
The authors would like to thank Rachel Bernier, MPH, and Erin Halpin, RN for supporting the study administration and data collection at Boston Children's Hospital and would also like to thank the children who assented and the parents and legal guardians who consented for their children to participate in this study. APPENDIX Patient information is summarized in Table V.