Convolutive Blind Source Separation on Surface EMG Signals for Respiratory Diagnostics and Medical Ventilation Control

TL;DR

This paper introduces a convolutive blind source separation method to effectively distinguish respiratory and cardiac muscle activities in surface EMG signals, enhancing respiratory diagnostics and ventilation control.

Contribution

It presents a novel broadband BSS algorithm capable of separating muscle activities in sEMG signals without restrictive assumptions, improving clinical analysis.

Findings

Successful separation of inspiratory, expiratory, and cardiac muscle signals

Broadband BSS algorithm achieves clear component discrimination

No restrictive assumptions on the demixing system

Abstract

The electromyogram (EMG) is an important tool for assessing the activity of a muscle and thus also a valuable measure for the diagnosis and control of respiratory support. In this article we propose convolutive blind source separation (BSS) as an effective tool to pre-process surface electromyogram (sEMG) data of the human respiratory muscles. Specifically, the problem of discriminating between inspiratory, expiratory and cardiac muscle activity is addressed, which currently poses a major obstacle for the clinical use of sEMG for adaptive ventilation control. It is shown that using the investigated broadband algorithm, a clear separation of these components can be achieved. The algorithm is based on a generic framework for BSS that utilizes multiple statistical signal characteristics. Apart from a four-channel FIR structure, there are no further restrictive assumptions on the demixing system.

Equations

$$\mathbf{x}(n)=\mathbf{H}\mathbf{s}(n),$$

$$\mathbf{x}(n)=[x_1(n),\dots ,x_P(n){]}^T\in {\mathbb{R}}^{P\times 1},$$

$$\mathbf{H}=\left(\begin{array}{ccc}{\mathbf{h}}_{1,1}^T& \cdots & {\mathbf{h}}_{1,Q}^T\\ ⋮& \ddots & ⋮\\ {\mathbf{h}}_{P,1}^T& \cdots & {\mathbf{h}}_{P,Q}^T\end{array}\right)\in {\mathbb{R}}^{P\times MQ},$$

$$\mathbf{s}(n)=[{\mathbf{s}}_1^T(n),{\mathbf{s}}_2^T(n),\dots ,{\mathbf{s}}_Q^T(n){]}^T\in {\mathbb{R}}^{MQ\times 1},\text{where}$$

$${\mathbf{s}}_q(n)=[s_q(n),s_q(n-1),\dots ,s_q(n-M+1){]}^T\in {\mathbb{R}}^{M\times 1}.$$

$$x_p(n)=\sum _{q=1}^Q{\mathbf{h}}_{p,q}^T{\mathbf{s}}_q(n)\approx \sum _{q=1}^{\overline{Q}}\overline{\mathbf{h}}{}_{p,q}^T\overline{\mathbf{s}}{}_q(n)\forall p,n.$$

$$y_q(n)=\sum _{p=1}^P{\mathbf{w}}_{q,p}^T{\mathbf{x}}_p(n),$$

$$\mathbf{ED}{\mathbf{E}}^{\mathrm{T}}$$

$${\left[x_1(n),\dots ,x_P(n)\right]}^{\mathrm{T}}\leftarrow \mathbf{E}{\mathbf{D}}^{-1/2}{\mathbf{E}}^{\mathrm{T}}\mathbf{x}(n)$$

$$X_p^{(\nu )}(m)\Leftarrow x_p(n),p=1,\dots ,P$$

$$X_p^{(\nu )}(m)\leftarrow X_p^{(\nu )}(m)-\frac{1}{N}\sum _{m^{\prime }=0}^{N-1}{X}_p^{(\nu )}(m^{\prime })$$

$${\mathbf{Y}}^{(\nu )}(m)$$

$$b_p(m)=\sqrt{\frac{1}{M}\sum _{\nu =0}^{M-1}|Y_p^{(\nu )}(m){|}^2}$$

$${\mathbf{\Phi }}^{(\nu )}(m)={\left[Y_1^{(\nu )}(m)/b_1(m),\dots ,Y_P^{(\nu )}(m)/b_P(m)\right]}^{\mathrm{T}}$$

$${\mathbf{W}}^{(\nu )\ell }$$

$${\mathbf{W}}^{(\nu )\ell }\leftarrow \mathrm{diag}\left\{{\left({\mathbf{W}}^{(\nu )\ell }\right)}^{-1}\right\}{\mathbf{W}}^{(\nu )\ell }$$

$$\mathbf{W}(n)\Leftarrow {\mathbf{W}}^{(\nu ){\ell }_{\max }},\nu =0,\dots ,M-1$$

Full text

1 Introduction

Blind Source Separation (BSS) plays an important role in processing surface electromyographic (sEMG) measurements [1] and can be employed to detect and distinguish crosstalk between adjacent muscles [1]. Furthermore, several research groups recently have reported the successful removal of electrocardiographic (ECG) interference from EMG signals by means of BSS techniques [2]. Both features are crucial in respiratory applications since they allow to discriminate between the activity of inspiratory and expiratory muscles and to remove the strong interference resulting from the electrical activity of the heart.

The mixing model assumed in a particular BSS algorithm is a major determinant of separation success. For the separation of EMG signals, mostly instantaneous mixing models have been employed in the past, despite the fact that these models are only applicable in very specific settings [3]. Convolutive mixing models, on the other hand, appear to provide a much more realistic representation of EMG signal propagation paths [1] but have only rarely been used in applications [4, 5], and to the author’s knowledge never been applied to the analysis of respiratory EMG. In contrast to the muscles usually considered for sEMG signal processing, respiratory muscles such as the diaphragm and the abdominal muscles are rather large and flat. This observation distinguishes the proposed respiratory application from various other EMG-based signal processing settings and entails important consequences for the choice of the mixing model, which will be the topic of section 2 of this article.

Following the discussion of the different mixing models for BSS of respiratory EMG signals in section 2, we describe a broadband algorithm for convolutive BSS that is based on a generic framework for the simultaneous exploitation of several statistical signal characteristics (TRINICON, TRIple-N Independent component analysis for CONvolutive mixtures) [6, 7] in section 3. Results of the application of this algorithm to a clinical data set are presented in section 4, before section 5 concludes this article.

2 Mixing Models for Source Separation of Surface EMG Signals

The contraction of a muscle fiber is triggered by an initial depolarisation of the muscle fiber membrane at the neuromuscular junction (NMJ), which then propagates along the fiber in the form of an intracellular action potential (IAP) at a speed of $3$–$5$\mathrm{m}\mathrm{/}\mathrm{s}$$. [8]. During the IAP’s propagation along the fiber, the fiber represents a distributed current source and sink that evokes a change in the skin surface potential. The electric field induced by this current source can be considered quasi-static [9], i.e., changes are propagated nearly instantaneously from the fiber to the skin. Figure 1 shows a simulation of the surface potentials resulting from a single activation of two muscle fibers at different depths, as measured by surface electrodes at different locations along the fibers.

Notably, the increased damping induced by the detection of a deeper fiber affects the end-of-fiber artifacts much less than the propagating signal components. This is due to the spatial low-pass behavior of the volume conductor and the end-of-fiber artifacts’ low spatial frequency [8].111The additional end-of-fiber signal at around $18\text{\,}\mathrm{m}\mathrm{s}$ is due to the two fiber halves being of different lengths - at that point in time, the wave front travelling along the opposite direction from the NMJ has reached the other fiber end. This discrepancy between the damping of the different signal components with increasing distance can not be represented by means of an instantaneous mixing model. A convolutive model, on the other hand, can precisely represent this behavior when we consider the input signal $s(t)$ to be an impulse train, i.e., a digital signal indicating whether the motor neuron fires an impulse at instant $t$. This mathematical formulation allows to represent the distributed current source by means of a point source, while all information on the spatial behavior is contained in the coefficients of the mixing filter.

In the most general setting, a $P$-channel EMG signal can be modelled as the sum of the contributions of all Motor Units222A MU denotes the ensemble of a motor neuron and all muscle fibers innervated by that neuron. Note that these fibers always fire synchronously. (MUs) belonging to the muscles of interest. In terms of the Motor Unit Action Potentials (MUAPs), i.e., the changes in surface potential that a single contraction of a particular MU evokes at a particular electrode, this model can be formulated in discrete time as

[TABLE]

with the measurement vector

[TABLE]

the mixing matrix

[TABLE]

where $\mathbf{h}_{p,q}=[h_{p,q,0}\cdots h_{p,q,(M-1)}]^{T}\in\mathbb{R}^{M\times 1}$ contains the time course of the length-$M$ MUAP evoked at sensor $p$ by motor unit $q$, and with the source vector

[TABLE]

[TABLE]

In all practical applications, $P\ll Q$ and hence the problem is strongly underdetermined. One therefore usually wishes to approximate the system by a much smaller number of sources $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu\ll Q$. This is equivalent to assuming

[TABLE]

Now, on the one hand, due to the large spatial extension of the respiratory muscles, different MUs belonging to the same muscle can produce MUAPs of strongly varying wave shape, depending on the relative position of the MU and the recording electrode. On the other hand, the muscles are relatively well-separated, in contrast to, e.g., forearm muscles. Hence, the MUAPs of MUs belonging to one respiratory muscle may vary in shape, but they will all be rather different from those of the other respiratory muscles. This argument provides an informal justification of approximation (5) when considering the $\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{h}\mkern-1.5mu}\mkern 1.5mu_{p,q}$ representative of all MUs belonging to a particular muscle, and the $\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{s}\mkern-1.5mu}\mkern 1.5mu_{q}(n)$ a superposition of the source signals of all MUs represented by $\mkern 1.5mu\overline{\mkern-1.5mu\mathbf{h}\mkern-1.5mu}\mkern 1.5mu_{p,q}$.

For source separation, i.e., for eliminating the crosstalk between the channels, a suitable MIMO demixing system needs to be designed, yielding the output signals

[TABLE]

where $\mathbf{w}_{p,q}\in\mathbb{R}^{L\times 1}$ are the coefficients of the demixing filters. The $P\cdot\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu\cdot L$ filter coefficients $w_{p,q,\kappa},\kappa=0,\dots,L-1$ are to be blindly estimated by the BSS algorithm. It can be shown that for an ideal separation the number $P$ of sensors must be at least equal to the assumed number $\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$ of sources [7]. Hence, in this study we choose $P=\mkern 1.5mu\overline{\mkern-1.5muQ\mkern-1.5mu}\mkern 1.5mu$.

3 TRINICON-based Blind Source Separation for Convolutive Mixtures

The most common method for BSS is the so-called Independent Component Analysis (ICA), for which several different algorithms have been developed in recent years. However, the original ICA model is based on instantaneous mixtures (i.e., $L=1$) [11]. To cover convolutive mixtures, i.e., to blindly estimate all filter coefficients $w_{p,q,\kappa}$, our examinations are based on a more general approach. The employed concept called TRINICON (TRIple-N ICA for CONvolutive mixtures) is a framework for broadband adaptive MIMO signal processing algorithms [6, 7] and allows the simultaneous exploitation of all fundamental statistical signal properties, namely nongaussianity, nonwhiteness and nonstationarity. For convolutive BSS these signal properties are used to minimize the mutual information known from information theory, i.e., the crosstalk between all output signals of the demixing system. It can be shown that in contrast to the instantaneous case $L=1$, the convolutive case $L>1$ requires the simultaneous exploitation of at least two of the three signal properties. In this study, we apply a specially adapted and efficient realization of the concept using fast Fourier transforms (FFTs).

In addition to the modeling by convolutive mixtures, another challenge in the present biomedical application is the significant level-difference (typically three orders of magnitude) between the ECG component and the relevant EMG signals. We therefore apply a two-step approach which consists first of a so-called sphering, followed by the actual source separation as described above. The sphering stage, consisting of an eigenvalue decomposition of the spatial correlation matrix and a power normalization, turns out to be a crucial step in handling this level imbalance so that it leads to a decisive numerical improvement of the subsequent separation task.

Both the sphering stage as well as the separation stage can be derived rigorously from the TRINICON framework. In the present study, they are both carried out in an offline fashion, i.e., calculations are performed with the complete dataset already available. The efficient frequency-domain realization shown here is based on the use of a multivariate probability density function so that cross-correlations between all frequency components are taken into account and the commonly known permutation ambiguity in traditional narrowband frequency-domain BSS algorithms is avoided. It follows directly from the generic broadband frequency-domain coefficient update equation (10.72) in [6] by approximating the constraint matrices $\mathbf{L}$ and $\mathbf{G}_{\cdot}^{\cdot}$. The coupling of the frequency components is ensured by using a spherical multivariate Laplacian density. (See also [7], p. 130, second column in fig. 4.10, for a classification of this realization.)

Algorithm 1 Convolutive Blind Source Separation based on Independent Vector Analysis

1:

2:Spatial prewhitening

3::

4:Principal component analysis

[TABLE]

5:Sphering:

[TABLE]

6:

7:Main algorithm

8::

9:Short-time Fourier transform by windowing:

10:

[TABLE]

11:Centering:

12:

[TABLE]

13:for $\ell=1,\ldots,\ell_{\mathrm{max}}$ do:

14:    Circular convolution to obtain (preliminary) outputs:

[TABLE]

15:    Computation of broadband normalization factors:

[TABLE]

16:    Multivariate score function:

[TABLE]

17:    Coefficient update:

[TABLE]

18:    Minimum distortion principle:

[TABLE]

19:Transformation of MIMO filter coefficients to time domain:

[TABLE]

20:($n=$ time index, $m=$ (temporal) block index, $\nu=$ frequency index, $p=$ channel index, $N=$ number of blocks, $M=2L=$ number of frequency bands, $P=$ number of channels, $\mu=$ step size, $\ell=$ iteration index)

4 Experimental Results

The BSS algorithm described in the previous section has been applied to actual sEMG measurement data.333Data collected by Philippe Jolliet and his team during a clinical study entitled “Comparative Effects on Diaphragmatic Electrical Activity and Respiratory Pattern of Various Levels of Assistance”, approved by the local ethics committee at Lausanne University Hospital (identifier: NCT01248845). Figure 2 shows the four EMG channels that have been recorded - note the extremely dominant ECG component in the upper three channels.

Figure 3 shows the output signals generated by the BSS algorithm.

As desired, the ECG component (channel 1) seems to be reasonably well separated from the other signal sources. Moreover, output channel 2 seems to capture expiratory muscle activity - supposedly from the abdominal muscles - while channels 3 and 4 show inspiratory muscle activity. One might furthermore speculate that channels 3 and 4 most likely represent intercostal and diaphragmatic activity. Note that there is no reference signal that could be used to validate the separated source signals (except for pneumatic signals, as shown), since the real level of activity of the different muscles is unknown.

The necessity of convolutive mixtures could also be confirmed experimentally on this data set. For the results shown here, we empirically determined a necessary filter length of $L=64$.444The filter length was varied in powers of two in order to obtain an efficient BSS implementation using FFTs. In particular, it turned out that in the simple case of an instantaneous mixing model (i.e., $L=1$) it was not possible to obtain any significant separation. Note that the filter length $L=64$ appears reasonable from a physiological point of view: Each demixing filter effectively represents a propagation path of total length $\Delta s=vLT_{\mathrm{a}}\approx$25\text{,}\mathrm{c}\mathrm{m} (assuming $v=$4\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}, see above). This is consistent with the dimensions of the muscle-sensor setup.

5 Conclusion

In this article we have presented evidence that the use of convolutive BSS may significantly advance the preprocessing of surface EMG measurements of respiratory muscles. Despite the limited validation data set, it can be concluded that the proposed algorithm has the potential of separating activity between inspiratory and expiratory muscles and removing ECG artifacts, both significantly better than previously proposed algorithms for instantaneous BSS. While in the present study an offline algorithm was used to compute the coefficients of the demixing filter, an adaptive variant of the proposed algorithm is immanent: instead of only computing the coefficients of the demixing filter initially before performing the actual demixing, an online update of the filter coefficients using block or incremental updates can be performed, allowing for real-time applications. The performance of the proposed algorithm will be evaluated on simulated data as well as on a wider set of clinical data in a subsequent study.

Acknowledgments

We thank Philippe Jolliet and his team at the Lausanne University Hospital for providing the clinical data used for the evaluation of the proposed algorithm.