An EMD-Based Algorithm for Emboli Detection in Echo Doppler Audio Signals

Divers’ health state after underwater activity can be assessed after the immersion using precordial echo Doppler examination. An audio analysis of the acquired signals is performed by specialist doctors to detect circulating gas bubbles in the vascular system and to evaluate the decompression sickness risk. Since on-site medical assistance cannot always be guaranteed, we propose a system for automatic emboli detection using a custom portable device connected to the echo Doppler instrument. The empirical mode decomposition method is used to develop a real-time algorithm able to automatically detect embolic events and, consequently, assess the decompression sickness risk according to the Spencer’s scale. The proposed algorithm has been tested according to an experimental protocol approved by the Divers Alert Network. It involved 30 volunteer divers and produced 37 echo Doppler files useful for the algorithm’s performances evaluation. The results obtained by the proposed emboli detection algorithm (83% sensitivity and 76% specificity) make the system particularly suitable for real-time evaluation of the decompression sickness risk level. Furthermore, the system could also be used in continuous monitoring of hospitalized patients with embolic risks such as post surgery ones.


Introduction
Diving related illness has become a public health concern [1] as diving is worldwide growing in popularity [2]. Amateur diving has changed decompression sickness (DCS) from being a specialized problem to becoming a health risk for many thousands of people. DCS is a well-known risk for divers equipped with self-contained underwater breathing apparatus (compressed air or artificial gas mixture), although its reported incidence in the diving community is generally low (about 1%) [3]. DCS is thought to result when rapid pressure reduction (e.g., during ascent from a dive, ascent to altitude or in aerospace-related events) causes gas previously dissolved in blood or tissues to form bubbles in blood vessels [4]. Gas bubbles may originate from inert gas supersaturation or the traumatic injection of gas into the arterial circulation following pulmonary barotrauma [5]. The main objective is to provide a system capable of acquiring, processing and transmitting remotely the signal of interest. Thus, the trained staff in the emergency central can continuously monitor a multiplicity of people and can organize the intervention in the most appropriate way [6].
The present paper focuses on decompression associated with the sudden decrease in pressure during underwater ascent, usually occurring during free or assisted dives. During a dive, the body tissues absorb nitrogen from the breathing gas in proportion to the surrounding pressure. As long as the diver remains at pressure, the gas presents no problem. If the pressure is reduced too quickly, nitrogen comes out of solution and forms bubbles in the tissues and bloodstream. When high levels of bubbles occur, complex reactions can take place in the body, usually in the spinal cord or brain. If large numbers of bubbles enter the venous bloodstream, congestive symptoms in the lung and circulatory shock can then be occurred. When this happens, the quicker the treatment begins, the better the chance for a full recovery. The early management of DCS is the recompression [7]. Since no specific tests for DCS exist, acute DCS is a purely clinical diagnosis that requires a fair amount of clinical suspicion to avoid missing cases [8]. Many preventive studies have been conducted [9] but further efforts focused on early detection of DCS are necessary.
In several previous studies on the correspondence between venous gas emboli (VGE) and the risk of developing DCS for various modes of decompression [10], it has been shown that a diver can have a large quantity of bubbles without any symptoms of DCS. As a consequence, the presence of gas bubbles alone is of no diagnostic value in individual cases [11]. On the other hand, most of the studies also show that absence of detectable bubbles is a good predictor of decompression safety. Thus, if an association between gas bubbles and DCS risk could be established with some degree of accuracy, bubble detection could be used as a tool for validation of the safety of decompression procedures.
Gas bubbles in liquids are strong reflectors of sound so ultrasound-based methods are well suited for detection of circulating vascular gas bubbles. Doppler [12] and in particular ultrasound Doppler [13,14] is employed both in medical and in industrial applications for estimating fluid velocity. Doppler systems are most commonly used [15] and ultrasonic detection of circulating VGE after diving is considered a useful index for safe decompression [16]. Doppler ultrasound (DU) monitoring provides an audio signal that, through aural analysis, allows a trained observer to detect the number of embolic event in divers. Only moving bubbles, i.e., intra-vascular bubbles, may be detected by this method. Doppler shift signals also arise from other moving objects in the 3D volume of blood sensed by the ultrasound probe, including blood cells, heart valves and muscle contraction.
Aural detection of decompression-induced bubbles is usually performed in the precordial region (the preferred location being the pulmonary artery) or in peripheral veins such as the subclavian or femoral [16].
Precordial DU is the most sensitive non-invasive monitoring [17]. Theoretically, precordial monitoring allows for an estimate of bubble production in the entire venous system, as all the returning blood from the body must flow through the pulmonary artery. However, in the precordial region background noise is ever present, created by the heart wall producing heartbeats, valves and greater blood flow. On the other side, peripheral monitoring only provides information regarding the bubble production in one part of the body showing a significant decrease in background noise [18]. Despite the reliability of the human ear, aural scoring is known to be observer-dependent [19]. In order to train and facilitate the tasks of an observer, an automatic tool for the objective analysis of echo Doppler signals is necessary to avoid interpretation errors [20].
Doppler signals from blood vessels contain both the echo related to bubble events and background noise. Noise in echo Doppler signal can be eliminated through dedicated hardware solutions such as the proposal in [14]. The background signal may be classified as a non-stationary but almost periodic process, whereas the bubble events are transient non-stationary random processes [21][22][23][24]. The Finite impulse response (FIR) technique with windowed Fourier series method (also called short-time Fourier transform, STFT) has been proposed for subclavian Doppler ultrasound VGE signals [25]. One of the pitfalls of STFT is that both time and frequency resolutions become fixed for all frequencies and times respectively [26]. Multiresolution analysis (MRA) methods, such as Wavelet decomposition, overcome most of the limitations of the Fourier transform. Wavelet transform is suitable for non-stationary analysis as it allows us to separate the frequency content of a signal without losing the information in the time domain [27,28]. Fast discrete wavelet transform (DWT) has been used to analyze precordial [29] and transcranial [22,30] DU audio signals. From experimental results [31], it was found that almost all gas bubble signals present in transcranial Doppler audio can be detected by the wavelet coefficients in the fourth scale. One drawback of using DWT for gas bubble detection is the reduced frequency resolution at lower scales, in which gas bubble signal is found mostly.
Another method for analyzing non-stationary, non-linear signals is the empirical mode decomposition (EMD) [32,33]. EMD is an adaptive data-driven technique allowing to localize any event in the time as well as in the frequency domain. Like the Fourier or wavelet transforms, EMD reduces a time signal into a set of basis signals; unlike STFT or DWT, however, basis functions are derived from the signal itself. Each basis function of the EMD, known as an intrinsic mode function (IMF), captures the energy associated with a particular time scale. Real world signals such as DU signals have often multiple causes and each of these causes may happen at specific time intervals. This type of information is evident in a set of IMFs, but quite hidden in STFT or DWT analysis. EMD has been applied for the off-line detection of VGE in the precordial region [15]. Such an approach relies on the regular occurrence of artifacts in the background signal compared to the transient nature of the embolic signals. EMD is more suitable for detecting single bubbles of big entity (i.e., greater than 25 microns) than other algorithms, but small or showers of bubbles may be undetectable. Moreover, although excellent sensitivity has been achieved with this method when run on precordial Doppler audio recordings, specificity (false positives) was reported as a problem [15]. To overcome such limitations, a new EMD-based approach for analyzing echo Doppler audio recorded in precordial region is presented. The proposed method employs an adaptive threshold aiming at providing quantitative information about individual bubble occurrence in DU data, unlike the semi-quantitative information available from aural grading.

Materials and Methods
Since a unique clinical acquisition protocol for the divers in post immersion does not exist, the Divers Alert Network (DAN) decided to set up one [34]. The DAN has started an extensive acquisition campaign with the purpose to obtain the diver's profiling for security reasons, the identification of individual risk, more detailed decompression tables and in order to monitor health state of the diver after the dive. According to the DAN, a DU monitoring has widespread application as it is non-invasive, relatively inexpensive and does not use ionizing radiation. In this study we start from a dataset of DU recordings provided by the DAN. The dataset consists of 37 audio files acquired by 30 divers, 60% male and 40% female from 25 to 65 years old. DU audio recordings were obtained from an extensive diving research project where the safety of the decompression profiles was assessed through the ultrasound registration. Such an acquisition campaign provided DU recordings from Maldive and Madagascar areas. Sensitive data of the dives have been acquired for each diver, both professional and amateur. All individuals who volunteered for the campaign gave their informed consent for inclusion before they participated in the study. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by DAN medical team.
Divers were monitored after surfacing using a 2 MHz Doppler probe (FD1, Huntleigh Ltd., Cardiff, UK). Uncompressed audio recordings were stored in linear pulse-code modulation (LPCM) format using a digital recorder (Tascam DP-004, TEAC America Inc., Santa Fe Springs, California, USA) not affected by audio saturation during the acquisition. Such devices were chosen for their relatively low price-performance ratio. This device is well suited for real-time monitoring of bubble events thanks to long life battery (about 500 min for both, Huntleigh FD1 and Tascam DP-004) and small size (140 mm × 70 mm × 27 mm-Huntleigh FD1, 155 mm × 33.5 mm × 107 mm-Tascam DP-004) and weight (295 g-Huntleigh FD1, 360 g-Tascam DP-004). The device is easy to use as it does not require any preliminary setup or configuration. Its technical specifications (sensitivity 92 dB/mW and frequency response in the 10-22,000 Hz range) guarantee optimal performances for use with both musical instruments and audio playback. The acquisition protocol consisted of two 45 s long acquisitions separated by a 10 s pause. Most of the bubbles were expected to be detected during the first acquisition, performed when the diver surfaces. Since showers of bubbles may remain trapped in tissues and thus cannot be detected by the echography, divers were asked to perform a series of three compressions bending the legs. Indeed, physical effort allows the bubbles to break away from the tissues and to flow into the blood stream. The second 45 s Doppler acquisition is performed in order to detect such bubbles. Each audio recording provided by the DAN [34] consists of the three acquisition phases described.

Bubbles Detection Algorithm
In this study, a signal processing algorithm for automatic emboli detection from echo Doppler audio recorded in precordial region is presented. The proposed algorithm is based on the EMD technique which is able to separate the contribution of bubble movement from other sources of sound in DU signal [15]. The EMD technique is based on the decomposition of a signal into a series of non-stationary elemental components called IMFs. IMFs are all of the same length and form a complete and nearly orthogonal basis for the original signal [32]. The decomposition takes place in the time domain then each IMF allows for frequency information to be preserved. Decomposition results in a set of ordered IMF components. Each successive IMF contains lower frequency oscillations than the previous one. Given a signal x(t), the effective algorithm of EMD can be summarized as [32,35]: 1. identify all extrema of x(t); 2. interpolate between minima (and maxima) ending up with envelope e min (t) (and e max (t)); 3. compute the mean r(t) between e min (t) and e max (t); 4. extract the detail signal c(t) = x(t) − r(t); 5. iterate on the residual signal r(t).
Interpolation in step 2 was performed using cubic splines [32]. The detail signal c(t) is referred to as an IMF if it satisfies the following requirements: • the number of IMF extrema (the sum of the maxima and minima) and the number of zero-crossings must either be equal or differ at most by one; • at any point of an IMF, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima shall be zero.
The above procedure has to be refined by a sifting process [32] which consists of iterating steps 1 to 4 S times, until c(t) can be considered as zero-mean according to a stopping criterion. Once this is achieved, the corresponding residual is computed and step 5 applies. Huang empirically established that an optimal sifting was obtained when 3 ≤ S ≤ 5 [36] and Chappell used S = 3 [15]. The signal x(t) can then be represented in the form: where n is the number of IMFs extracted from x(t). Theoretically, the decomposition process ends when the residual signal becomes a monotonic function. The EMD of one of DU signals used for this work is reported in Figure 1.
Although the EMD method was not designed as a sparse decomposition method, it can be used for this purpose [37]. In a sparse representation of a signal, only a few coefficients are significantly different from zero [38]. In order to separate these nonzero coefficients from the others, the proposed method employs an adaptive threshold chosen to minimize Stein's unbiased risk estimate (SURE) [39]. The term "adaptive" refers to the fact that SURE threshold is data-driven. A block diagram representing the algorithm for automatic detection of gas bubbles is shown in Figure 2.

Pre-Processing of Input Signal
Pre-processing has been performed on the original audio recordings in order to automatically filter out the features that can negatively affect the detection performances. Useless parts in Doppler audio files are the central pause and first and last seconds of acquisitions. Initial and final seconds of audio are related to probe positioning so they are characterized by sudden peaks of noise. The two acquisitions have been isolated from the central pause by identifying the noise peaks due to probe placement. For this purpose the signal was windowed using non-overlapping rectangular pulse functions of duration 0.1 s and selecting the maximum value within each window. A zoom on original Doppler audio signal is shown in Figure 3 together with the windowing function.  In order to distinguish the initial transient from the first acquisition, the differences between consecutive windows in normalized audio signal were calculated. The first useful window in the first part of DU signal was identified as the one for which the difference exceeds a threshold of 0.3. Such a threshold has been empirically selected by analyzing the initial transient of all 37 audio files provided by the DAN. The proposed algorithm did not take into account the last two seconds of audio files since no shower of bubbles is expected to be detected after 15 s following physical effort [34]. The frame with maximum amplitude within 50 s from the beginning of the audio signal has been identified as the end of the first acquisition. Samples contained in such a frame and in the previous two frames were considered as useless signal as they were related to probe motion. Analogously, the beginning of useful signal in second acquisition was identified as the first sample in the third frame following the one with maximum amplitude. Finally, the two segments of audio signal were sequentially merged into a single synthetic signal. Figure 4 shows an example of such a signal to be processed by the detection algorithm.

EMD Implementation
In the present study we chose parameters S and n in an empirical way. We performed EMD on a randomly chosen set of 25 DU signals with n = 3 using S from 1 to 20 for each IMF. Aural detection was also performed by 8 expert observers on such signals. By comparing the results of the tests, we found that large values for S (S ≥ 10) cause a reduction both in the number of bubbles detected and in the related Spencer level [40,41]. Therefore, with the aim of finding a unique value for S, applicable to the entire data set, various experimental tests have been carried out with S between 1 and 9. From the analysis of the different outputs it emerged that the optimal value for S is equal to 5. So we choose to use S = 5 as number of iteration that represents the termination rule of the algorithm. Moreover, we chose to calculate EMD with n = 1 as c 1 (t) contains all bubbles which are clearly distinguishable from heartbeats by aural analysis. We also noted the almost total absence of background noise between adjacent peaks in c 1 (t).

Threshold Calculation and IMF Windowing
In order to distinguish air bubbles from heartbeats and background noise, the computation of the embolic signal intensity to average background ratio (EBR) parameter is needed. We evaluated EBR as the peak to threshold ratio (P2TR) as suggested by Aydin [22]: The threshold A th was set by borrowing ideas from Donoho and Johnstone [42] who studied the problem of nonlinear estimation of signals under a sparse representation. In particular, we used a data-driven threshold starting from input signal x(t) [43]: where N is the length of x(t) and σ B is the noise parameter estimated by means of the standardized median absolute deviation (MAD) of x(t) [44] and defined as in (4): Assuming that the minimum duration of a bubble event is 10 ms [25], the signal c 1 (t) has been windowed with non-overlapping 10 ms long rectangular pulses. The signal A peak was then found by considering the maximum value of windowed c 1 (t) in each frame.

Bubble Detection Refining and Determination of Embolism Risk
For each frame in A peak , a bubble event was detected each time P2TR exceeded the threshold A th and then stored in a two-dimensional array structure together with the corresponding time instant of occurrence. Since bubble duration can be longer than a single frame length, in this way embolic events occurred in adjacent frames can be easily grouped as a single bubble event. Shower identification and determination of individual gas embolism risk have then been performed on DU signals. A shower may be defined as a continuous flow of bubbles occurring as a result of some form of muscular contraction. In this work, each pair of embolic events occurred in distinct frames with maximum spacing in time of 20 ms were grouped to form a shower of bubbles. The last part of the proposed detection algorithm is about the determination of Spencer level and the corresponding individual embolism risk [18,40,41]. It is possible as the proposed algorithm has been designed to store the instant of occurrence as well as the duration of all bubble events and showers.

Results
The proposed signal processing algorithm has been evaluated using the dataset provided by DAN. In particular, the performance has been verified by comparing the provided output with the notes handwritten by experts listening the same echo Doppler files. File's annotations report are provided by independent blind teams. DAN medical experts stated that the estimate provided by each blind team can differ of ±0.5 degrees around the presumed real value. Therefore, ±0.5 represents the real error margin in manual annotations made by each blind team. It represents a compromise that does not invalidate the reliability of the measure and still provides a reliable indication on the decompression sickness risk. According to a joint analysis with DAN medical experts, we assumed a report was acceptable if the Spencer level estimated by automatic algorithm did not differ more than ±0.5 degrees on the Spencer level from the assessment made by the blind teams. The comparison between the Spencer level manually performed by DAN and the automatically performed by the algorithm is shown in Table 1. The manual reports, provided by DAN, contain, in addition to the Spencer level, the instant of occurrence of each bubble in each files, through which we calculated the percentage of specificity and sensitivity of the algorithm based on the Table 2. The Sensitivity S e and the Specificity S p were obtained using the Equations (5) and (6). The average S e obtained for all tested file is 83% and the average S p is 76%.

Discussion
The results show that the proposed algorithm performance are suitable for the automatic estimation of decompression sickness risk. In one case it misclassifies a DU signal with a Spencer level of zero. Analyzing this particular case emerged that the overestimation of the Spencer level is due to the noises introduced during the re-positioning of the echo Doppler probe after the compressions bending the legs in the second part of the test. These noises are identified as bubbles by the algorithm leading to an incorrect estimate of decompression sickness risk. In other analyzed cases where there is no complete correspondence with the assessments provided by the three blind teams (matching with threshold) the algorithm overestimates the Spencer level due to the noises and the manual positioning of the probe.
In order to validate the obtained results from the EMD technique, the power spectral density (PSD) of DU signals was calculated. The main idea of using PSD is to detect a change in frequency content for signals with and without gas bubbles. In order to calculate a reference PSD (PSD re f ) we need to choose a frame with no gas bubbles within DU signals. Starting from the whole signal x(t) we applied a 10 ms-long (441 samples) Hamming window to x(t) with 256 overlapping samples and computed FFT on such a signal. The presence of gas bubbles was evaluated using signals x(t), c 1 (t) and c 3 (t), where c 3 (t) represents the IMF obtained from EMD with n = 3 and S = 5. c 3 (t) was chosen since no bubbles are expected to be found in this signal while heartbeats are still audible. The three signals under test were windowed using non-overlapping rectangular pulse functions of 10 ms and the PSD was evaluated for each frame (PSD f rame ). Each PSD f rame was compared with PSD re f . When all components of PSD f rame exceeded PSD re f an embolic event was detected and the instant of occurrence was stored. Figures 5-7 show a comparison between PSD f rame and PSD re f for the three signals x(t), c 1 (t) and c 3 (t), respectively.
From the conducted tests, we conclude that c 1 (t) can be considered as the best base signal for bubble detection in DU recordings since it represents an optimal compromise in terms of number of detected bubbles and noise.
The increase in PSD f rame compared to PSD re f is mainly due to the background noise rather than the presence of bubbles. On the other hand, c 3 (t) is not to be considered for detection purposes since it provided a number of bubbles even greater than that obtained using signal c 1 (t).
Further tests have been performed using different widths for the windowing function of the base signal. In particular, we chose 20 ms, 50 ms, 0.1 s, 0.5 s and 1 s. We obtained that a value of 20 ms provides the best agreement between results of PSD-based algorithm and aural grading scores.   Comparison between PSD f rame and PSD re f for signal c 3 (t).

Conclusions
The automatic detection of bubbles is a particularly interesting problem that needs to be analyzed because its solution would open new ways for prevention and treatment of decompression sickness symptoms that affects divers and not only. Using the proposed EMD algorithm we properly identified the degree of embolism risk for over 83% of the set of available audio files. In order to have a reliability level of the EMD algorithm, we compared the results with those of a second algorithm based on PSD. This has led to an accurate match between results obtained from EMD and PSD algorithms and DAN annotations. The EMD algorithm is based on an adaptive threshold calculated using SURE method. Further analysis is needed to investigate the use of adaptive thresholds based on heuristic variance techniques or on the minimum and maximum mean square value. We found that PSD technique is able to identify signal regions in which there is more energy due to the presence of bubbles. This is influenced by the presence of heartbeats, and then it would be appropriate to test a series of algorithms that allow an accurate peak's detection of pre/post physical effort and exploit such data as an additional information to take in account. The algorithm is currently used and tested by the DAN Europe during its acquisition campaigns. Further acquisition campaigns have been planned by DAN to increase the dataset of files available with the aim of improving the performance of the developed algorithms. An embedded implementation of the algorithms described was also developed and provided to DAN for testing the whole system composed by acquisition unit, elaboration unit and transmission unit. An embedded version of the algorithm could be used in a telemedicine system to assist divers remotely. The algorithm results (Extended Spencer Scale, Spencer Scale or the number of bubbles and showers) can be send to an hyperbaric medical that can re-examine the echo Doppler file and confirm the diver therapy before field intervention. Funding: This research received no external funding.

Acknowledgments:
The authors are grateful to the DAN Europe Foundation for the opportunity to carry out this research work. They would like also to thank WiSense s.r.l. for its insight and expertise that greatly assisted this study.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: