This paper explores a novel approach for ventricular and atrial activities estimation in electrocardiogram (ECG) signals, based on sparse source separation. Sparse decompositions of ECG over signal-adapted multi-component dictionaries can lead to natural separation of its components. In this work, dictionaries of functions adapted to ventricular and atrial activities are respectively defined. Then, the weighted orthogonal matching pursuit algorithm is used to unmix the two components of ECG signals. Despite the simplicity of the approach, results are very promising, showing the capacity of the algorithm to generate realistic estimations of atrial and ventricular activities

Divorra Escoda, O.

Granai, L.

Lemay, M.

Molinero Hernandez, J.

Vandergheynst, P.

Vesin, J.

This paper explores a novel approach for ventricular and atrial activities estimation in electrocardiogram (ECG) signals, based on sparse source separation. Sparse decompositions of ECG over signal-adapted multi-component dictionaries can lead to natural separation of its components. In this work, dictionaries of functions adapted to ventricular and atrial activities are respectively defined. Then, the weighted orthogonal matching pursuit algorithm is used to unmix the two components of ECG signals. Despite the simplicity of the approach, results are very promising, showing the capacity of the algorithm to generate realistic estimations of atrial and ventricular activities.LTS

Infoscience - École polytechnique fédérale de Lausanne

VENTRICULAR AND ATRIAL ACTIVITIES ESTIMATION THROUGHSPARSE ECG SIGNAL DECOMPOSITIONSO. Divorra Escoda, L. Granai, M. Lemay, J. Molinero Hernandez†, P. Vandergheynst, J.-M. VesinEcole Polytechnique Fe´de´rale de Lausanne, Signal Processing InstituteCH-1015 Lausanne, Switzerland, http://itswww.epfl.ch{lorenzo.granai,mathieu.lemay,pierre.vandergheynst,jean-marc.vesin}@epfl.ch,oscar.divorra@ieee.org, javiermolinerohernandez@yahoo.esABSTRACTThis paper explores a novel approach for ventricular andatrial activities estimation in electrocardiograms (ECG) sig-nals, based on sparse source separation. Sparse decompo-sitions of ECG over signal adapted multi-component dictio-naries can lead to natural separation of its components. Inthis work, dictionaries of functions adapted to ventricular andatrial activities are respectively defined. Then, the weightedorthogonal matching pursuit algorithm is used to unmix thetwo components of ECG signals. Despite the simplicity ofthe approach, results are very promising, showing the capac-ity to generate realistic estimations of atrial and ventricularactivity.1. INTRODUCTIONAtrial fibrillation (AF) is the most common type of human ar-rhythmia and it is responsible for about one third of hospital-izations for arrhythmia problems. AF is more frequent in el-derly, as its prevalence doubles with each decade of age, from0.5% at ages between 50-59 years to almost 9% at ages be-tween 80-89 years. AF is an important clinical entity becauseof the increased risk of morbidity and mortality. The mostfrequent consequences are hemodynamic function impairment(loss of atrial synchronized contraction, irregular and inade-quately rapid ventricular rate), atriogenic thromboembolicevents and tachycardia induced atrial and ventricular car-diomyopathy. AF diagnosis has been assessed for years by vi-sual inspection of the surface electrocardiogram (ECG). Onthe ECG, the AF signals are characterized by continuous,apparently disorganized, fibrillatory waves (F-waves). Dueto the much higher amplitude of the electrical ventricular ac-tivity (VA) on the surface ECG, isolation of the atrial activity(AA) component in the ECG is crucial for the study of AF.Some methods used to solve this problem are based onaverage beat subtraction. These methods are built on the as-sumption that the AA is uncoupled with the VA. An averageof the ventricular complexes (QRST complexes) is then usedto subtract VA [2]. Other methods are based on indepen-dent component analysis (ICA). These try to find indepen-∗We thank Prof. Kappenberger (Lausanne University Hospital,Lausanne Switzerland) for providing the clinical ECG database.†J. Molinero Hernandez was an exchange student fromUniversitat Polite`cnica de Catalunya (UPC), Barcelona, Spain.The work presented in this paper was performed in the frameworkof his M.S. Thesis [1].dent components in an instantaneous linear mixture [3]. Amajor difficulty in ABS approaches is the limitation imposedby the use of a small number of VA average templates for gen-eral VA approximation. In present ICA based approaches, amajor gap is that only statistical priors are considered with-out taking into account the structural nature of signals. Inorder to circumvent these problems, a possible direction to ex-plore is the use of sparse source separation approaches basedon signal adapted redundant dictionaries.During the past decade, many important advances havebeen achieved in nonlinear signal approximation methods forsparse decompositions over redundant dictionaries (e.g. [4,5, 6]). In many applications, these techniques offer betterperformances than those based on orthonormal transforms ordirect time domain processing, thanks to their good capacityfor efficient signal modeling. In this paper, we present a novelapproach for VA and AA estimation. We explore a source sep-aration approach based on sparse decomposition of ECG sig-nals on a redundant multi-component dictionary. Such an ap-proach allows for the consideration of priors on the structuralnature of the different class of signals we are willing to sepa-rate. The multi-component dictionary is composed by func-tions specially designed to match the main structural char-acteristics of VA and AA signals. We also present WeightedOrthogonal Matching Pursuits (Weigthed-OMP) [6] as a toolfor generating ECG sparse approximations for source separa-tion.2. UNDERDETERMINED SPARSE SOURCESEPARATIONLet fmixj (t) : j ∈ [0,M − 1] be a set of M signal mixturesgenerated by the weighted superposition of N source signalsfi (t) : i ∈ [0, N − 1] such that:fmixj (t) =N−1∑i=0aj,i · fi (t) + nj(t), (1)where nj(t) represents some additive noise.The problem of source separation is a classical problem inmany fields like acoustics, radio or medical signal and imageprocessing. Many applications exist where the retrieval ofthe different additive components, forming a set of complexsignals, is required. The use of sparse signal representationsin source separation problems was proposed in [7] in order toexploit some prior knowledge about the structural character-istics of each fi (t). Many signals can be sparsely representedby the superposition of a limited number of atoms from anadapted dictionary of functions (D = {gl(t) : l ∈ Ω}):fi (t) =∑l∈Λbil · gil (t) +R|Λ|fi(t), (2)where bil are the atom coefficients, Λ ⊂ Ω and RΛfi(t) is aneventual residual that depends on Λ.Combining Eq. (1) and Eq. (2), one derives that, in orderto solve the separation problem, the set of mixing coefficientsaj,i and a sparse set of expansion coefficients bil must be re-covered while nj(t) are kept as small as possible. Indeed,accurate sparse models can efficiently capture the structuralnature of signals, leading to better source separation resultsas exposed in [7].A challenging form of sparse source separation problem iswhen there are fewer mixtures than sources. An underdeter-mined source separation problem is, for example, to separatethe different components from a single ECG channel trace.This is studied in this work through the use of a novel sparsesource separation approach. In order to do this, we adapt thetwo stage separation process proposed in [7] to the particularcase of ECG components estimation: First, we a priori designan overcomplete dictionary where sources are assumed to besparsely representable. Second, the sources are unmixed byexploiting their sparse representability.3. ECG COMPONENTS SEPARATIONIn this section, we formulate the ECG activities estimationaccording to the signal models and sparse source separationstrategy described in Sec. 2. The ECG signal (fECG) is mod-eled as a noisy mixture of the two cardiac activities of ourinterest (fAA and fV A):fECG = fAA + fV A + n, (3)where n stands for the noise.The generation of good sparse models for fAA and fV A re-quires the use of basis functions fitting their particular struc-tures. As shown in the following, fAA and fV A have quitedifferent characteristics,and this is what enables us to sepa-rate them. The approach we propose is based on the decom-position of fECG on a redundant dictionary (D) composedby the union of two sub-dictionaries: DV A suited for repre-senting the ventricular activity and DAA, better adapted forrepresenting the atrial activity. In the following D,DV A andDAA stand for the synthesis matrices of D,DV A and DAArespectively, where each column represents an atom of thedictionary. Hence,fECG ' D · b = DAA · bAA +DV A · bV A. (4)Given the noisy nature of fECG and the high complex-ity of each of its components, in this work, we exclusivelyconsider sparse approximations. According to Eq. (4), b iscomposed of two parts (bAA and bV A), each one containingthe coefficients related to DAA and DV A. The approach weexplore is simple: one generates a sparse approximation of50 100 150 200 250 300−0.4−0.3−0.2−0.100.1temporal axisamplitudeOriginalApproximation0 50 100 150 200 250 300 350 400 450 50000.10.20.30.40.50.60.70.80.91temporal axisamplitudeβ=1β=1.5β=2β=2.5Fig. 1. Left: QRST VA complex and its approximation usingjust 3 atoms. Right: Effect of β on the GGF (see Eq. (5)).fECG on D and then, estimates for fAA and fV A are recon-structed by just using the components from the appropriatedictionary:fAA ' DAA · bAA and fV A ' DV A · bV A.3.1. Modeling ECG Ventricular ActivityThe sub-dictionary DV A is generated by all possible transla-tions of the Generalized Gaussian Function:gV A(t) = C1 exp(−( |t− p|α)β), (5)where C1 is a normalizing constant, α determines the scaleand β the peakiness. This waveform allows to well approxi-mate the structure of a VA complex using few atoms. Fig. 1(left) shows a QRST complex and its approximation by usingjust 3 atoms. With respect to the simpler Gaussian func-tion, β allows us to approximate with more accuracy Q andR peaks (see Fig. 1 on the right).The possible values of α and β have been chosen exper-imentally after an extensive set of tests: α ∈ {3, 4, ..., 7} ∪{49, 50, ..., 54}, the first set adapted for Q and R waves andthe second intended for T waves approximation, while β ∈{1.5, 1.6, ..., 2.2}. Together with p, this makes DV A highlycoherent, but also very flexible for VA approximation. How-ever, such dictionary is far from being optimal, and severalimprovements are still possible, mainly concerning the ap-proximation of T waves.3.2. Modeling ECG Atrial ActivityThe sub-dictionary DAA is generated by all possible transla-tions of a real Gabor function:gAA(t) = C2 exp(−(t− pα√2)2)cos(2pik(t− p)N−∆ψ),where C2 is a normalizing constant, N is the signal length,α tunes the scale, k the frequency and ∆ψ the phase. Thiswaveform is specially adapted for AA approximation. Indeed,as can be observed in Fig. 2 (left), fibrillating AA is of os-cillatory nature, which is a perfect fit for the optimal spatio-temporal frequency localization of Gabor functions (see Fig. 2on the right).The values of the Gabor function parameters have beendetermined through an extensive analysis on several ECG sig-nals (see [1]). During the design of DAA, special care in lim-iting the maximum correlation between DAA and DV A atoms0 100 200 300 400 500 600 700−0.04−0.03−0.02−0.0100.010.020.030.040.050.06temporal axisamplitude0 50 100 150 200 250−0.15−0.1−0.0500.050.10.15temporal axisamplitudeFig. 2. Left: Example of a simulated AA wave during fibril-lation. Right: Gabor atom.was required. Indeed, excessive coherence between this twosub-dictionaries turns into a complete failure of the algorithm.4. ECG SPARSE DECOMPOSITION BYWEIGHTED-OMPEq. (3) involves an underdetermined problem that has nounique solution. The search for the sparsest approximationof fECG requires the exhaustive test of all coefficient possi-bilities, i.e. a combinatorial problem. Different alternativeapproaches have been proposed in order to make computa-tionally affordable the retrieval of a solution for b (see [4] fora review of some of them). These in general do not guaran-tee the recovery of the sparsest solution. Nevertheless, recentresults show that under certain conditions on the dictionaryand the signal, even these fast sub-optimal methods find thesparsest approximation [5].4.1. Weighted Orthogonal Matching PursuitWeighted-OMP [6] is used to decompose ECG signals. It isa weighted version of Orthogonal Matching Pursuit (OMP),that iteratively builds m-term approximants by selecting ateach step the most appropriate term from D according to aselection rule. Each iteration k : k ≥ 0 can be seen as a twostep procedure:1. A selection step where an atom glk ∈ D is chosen ac-cording to:glk = arg maxgl∈D|〈rk, gl〉| · wl,where wl ∈ [0, 1] is a pre-estimated weight that reflects,according to a predefined model, the a priori likelihoodthat atom gl may belong to the good set of atoms bestapproximating the signal.2. An orthogonal projection step where an approximantfk+1 ∈ span(gip : p ∈ {0, ..., k}), and a residual rk+1 =f −fk+1 are generated (notice that r0 = f , and rk+1 ⊥fk+1 ∀k). This stage is the one in charge to update, atevery step, the set of scalar expansion coefficients.The signal representation generated by Weighted-OMP is,thus, of the form of Eq. (2).Weighted-OMP was recently proved in [6] to outperformOMP when using coherent dictionaries and reliable prior in-formation. Weighted-OMP can consider, in the decomposi-tion algorithm, a priori models about the behavior of thedictionary in use with the class of signals to decompose. Theuse of good enough a priori models can reduce the instabilityof Weighted-OMP in front of OMP when trying to recoversparse approximations/representations.4.2. Weights Generation: Relation Between ECGA Priori Knowledge and the DictionaryThanks to the structure of VA, fECG can be divided in VAperiods. In addition, each VA period can be divided in aset of intervals corresponding to the different VA waves (Q,R, S and T) and an interval without ventricular activity. VAintervals can be estimated and identified, in practice, throughthe use of QRST point estimators (e.g. see [8], used in thiswork). This prior information can thus be used to generatewl ∀l. The a priori knowledge obtained from [8] needs to berelated with D in the following way.D is divided in DAA and DV A. According to fECG na-ture, AA can be found through all the VA period. Hence,DAA atoms can not be penalized in terms of this informa-tion. This is the reason why in this study we force: wl =1 ∀l : gl ∈ DAA. To the contrary, the selection of gl ∈ DV Acan be successfully influenced by the use of the available apriori information. DV A, as seen in Sec. 5, is composed by ablock optimized for QRS waves (ventricular depolarization)and a block designed for T waves (ventricular repolarization).Depending on the VA interval, wl can be set to 1 for everygl ∈ DV A belonging to the appropriate kind for that interval.In case a gl is unsuitable for a given interval, wl can be set toa penalizing factor 0 ≤ wp < 1. Thanks to the reliability ofthe estimators used in this work, it turned out that the bestvalue for wp in our experiments is 0.5. EXPERIMENTAL RESULTS5.1. ValidationA biophysical computer model of the atria was used to obtaina realistic atrial electrical activity on the torso [9]. The AFsignals that were generated in the 12-lead ECG were addedto a clinical 4-second standard 12-lead ECG of an AF parox-ysmal patient (78 years old) in sinus rhythm in which the Pwaves (AA) were removed. The clinical ECG was selectedto represent the VA in AF as closely as possible. The ratiobetween the power of the original signal (simulated AA) andthe estimation error (estimated AA - simulated AA) was usedto evaluate the performance of our method.5.2. ResultsFirst of all, we want to underline that we validated our choiceof Weighted-OMP instead of OMP with these simulated mea-sured 4-second ECG signals. By using Weighted-OMP, weincreased the SNR in the recovery of VA (respectively, AA)by 0.81 dB (respectively, 0.65 dB). All the following resultswere obtained by approximating ECG signals with 50 atoms.Fig. 3 shows the resulting separation of VA and AA of thesimulated measured 4-second ECG signal on lead V1. Onecan see how our method succeeds in approximating each oneof the VA periods separating, at the same time, the AA withsurprising fidelity. In order to study the influence of the AAamplitude on the method, three different simulated AA sig-nals were created; 50%, 100% and 150% of the original sim-ulated AA amplitude. The ratio between the power of the1 mV1 sec(a)(b)(c)(d)(e)Fig. 3. (a) Simulated measured 4-second ECG signal on V1.(b) Original VA on V1. (c) Estimated VA on V1 (SNR : 8.69dB). (d) Simulated AA on V1. (e) Estimated AA on V1 (SNR: 6.81 dB).original activity (VA or AA) and the error on the estimatedone was evaluated on leads VR, V1 and V4 (see Table 1). Wecan observe that the quality of the AA estimation depends onthe lead and its original amplitude. The AA SNRs are muchhigher with the 150% original amplitude and the overall per-formance on lead V1 is better than those of other two leads.Of course, the SNR values are directly related to the signalamplitude and in V1, the AA amplitude is higher compared toother leads. However, we observe a decrease of VA estimationperformance in lead V1. Fig. 4 shows the resulting separa-0.5·AA+VA 1·AA+VA 1.5·AA+VAlead VRVA SNR(dB) 11.06 10.88 11.08AA SNR(dB) -6.94 -1.05 2.61lead V1VA SNR(dB) 11.13 8.69 2.41AA SNR(dB) 3.61 6.81 4.28lead V4VA SNR(dB) 12.33 11.94 11.66AA SNR(dB) -6.53 -0.8 2.4Table 1. Signal-to-noise ratio (dB) on lead VR, V1 andV4. The performance of our method is tested on 3 differentAA amplitudes (50, 100 and 150 % of the original simulatedsignal).tion of the VA and AA of the clinical 4-second ECG signalon lead V2.A part from the visually satisfying componentsseparation, these resulting signals were validated using esti-mated power spectral densities. The dominant frequency ofVA (respectively, AA) is between 1 and 2.5 Hz (respectively,between 3 and 10 Hz). The fact that there is no presence ofVA dominant frequencies in the AA estimated PSD demon-strates the quality of our clinical results. Further results canbe found in [1].6. CONCLUSIONSIn this paper, we present a new framework based on sparsesource separation that can be used for QRST cancellation.Results appear to be very promising. Complementary worksX: 1.56 HzY: 0.0100 mV2/HzX: 1.56 HzY: 0.0100 mV2/HzX: 7.55 HzY: 0.0002 mV2/Hz1 mV1 sec(a)(b)(c)(d)2.5 mV1 sec(e)(f)(g)Fig. 4. (a) Clininal 4-second ECG signal on V2 with a dom-inant frequency of 1.56 Hz (see its PSD (e)). (b) EstimatedVA on V2 with a dominant frequency of 1.56 Hz (see its PSD(f)). (c) Estimated AA on V2 with a dominant frequency of7.55 Hz (see its PSD (g)). (d) Estimated AA on V2 magnified5 times.can be spent on the improvement of the dictionary design,specially concerning the T wave modelling. Moreover, morecomplex models for a priori knowledge can be studied andthe number of atoms used for decomposing signals can beoptimized.REFERENCES[1] J. Molinero Hernandez, “Sparse decompositions for ven-tricular and atrial activity separation,” M.S. thesis, Sig-nal Processing Institute, E´cole Polytechnique Fe´de´rale deLausanne, Switzerland, August 2005.[2] M. Lemay, V. Jacquemet, A. Forclaz, J.-M. Vesin, andL. Kappenberger, “Spatiotemporal QRST cancellationmethod using separate QRS and T-waves templates,” inComputers in Cardiology 2005. September 2005, IEEE.[3] F. Castells, J. Igual, J. J. Rieta, C. Sanchez, and J. Mil-let, “Atrial fibrillation analysis based on ICA includingstatistical and temporal source information,” in ICASSP,2003, pp. V93–V96.[4] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomicdecomposition by basis pursuit,” SIAM J. Sci. Comp.,vol. 20, no. 1, pp. 33–61, 1999.[5] D. L. Donoho, M. Elad, and V. Temlyakov, “Stable recov-ery of sparse overcomplete representations in the presenceof noise,” IEEE Trans. Inform. Theory, to appear.[6] O. Divorra Escoda, L. Granai, and P. Vandergheynst, “Onthe use of a priori information for sparse signal approxi-mations,” IEEE Trans. Signal Processing, to appear.[7] M. Zibulevsky and B. A. Pearlmutter, “Blind source sep-aration by sparse decomposition in a signal dictionary,”Neural Computation, vol. 13, no. 4, pp. 863–882, 2001.[8] R. M. Rangayyan, Biomedical Signal Analysis : A Case-Study Approach, John Wiley & Sons, Inc., 2002.[9] V. Jacquemet, M. Lemay, J.-M. Vesin, A. van Oosterom,and L. Kappenberger, “A biophysical model of ECG sig-nals during atrial fibrillation used to evaluate the perfor-mance of QRST cancellation algorithms,” in Computersin Cardiology 2005. September 2005, IEEE.

Ventricular and Atrial Activity Estimation Through Sparse ECG Signal Decompositions

This paper explores a novel approach for ventricular and atrial activities estimation in electrocardiogram (ECG) signals, based on sparse source separation. Sparse decompositions of ECG over signal-adapted multi-component dictionaries can lead to natural separation of its components. In this work, dictionaries of functions adapted to ventricular and atrial activities are respectively defined. Then, the weighted orthogonal matching pursuit algorithm is used to unmix the two components of ECG signals. Despite the simplicity of the approach, results are very promising, showing the capacity of the algorithm to generate realistic estimations of atrial and ventricular activities.LTS2IT

VENTRICULAR AND ATRIAL ACTIVITY ESTIMATION THROUGHSPARSE ECG SIGNAL DECOMPOSITIONSO. Divorra Escoda, L. Granai, M. Lemay, J. Molinero Hernandez†, P. Vandergheynst, J.-M. VesinEcole Polytechnique Fe´de´rale de Lausanne, Signal Processing InstituteCH-1015 Lausanne, Switzerland, http://itswww.epfl.ch{lorenzo.granai,mathieu.lemay,pierre.vandergheynst,jean-marc.vesin}@epfl.ch,oscar.divorra@ieee.org, javiermolinerohernandez@yahoo.esABSTRACTThis paper explores a novel approach for ventricular and atrial activities estimation in electrocardiogram (ECG)signals, based on sparse source separation. Sparse decompositions of ECG over signal-adapted multi-componentdictionaries can lead to natural separation of its components. In this work, dictionaries of functions adapted toventricular and atrial activities are respectively defined. Then, the weighted orthogonal matching pursuit algorithmis used to unmix the two components of ECG signals. Despite the simplicity of the approach, results are verypromising, showing the capacity of the algorithm to generate realistic estimations of atrial and ventricular activities.1. INTRODUCTIONAtrial fibrillation (AF) is the most common type of human arrhythmia and it is responsible for about one third ofhospitalizations for arrhythmia problems. AF is more frequent in elderly, as its prevalence doubles with each decadeof age, from 0.5% at ages between 50-59 years to almost 9% at ages between 80-89 years. AF is an important clinicalentity because of the increased risk of morbidity and mortality. The most frequent consequences are hemodynamicfunction impairment (loss of atrial synchronized contraction, irregular and inadequately rapid ventricular rate),atriogenic thromboembolic events and tachycardia induced atrial and ventricular cardiomyopathy. AF diagnosishas been assessed for years by visual inspection of the surface electrocardiogram. On the ECG, the AF signals arecharacterized by continuous, apparently disorganized, fibrillatory waves (F-waves). Due to the much higher amplitudeof the electrical ventricular activity (VA) on the surface ECG, isolation of the atrial activity (AA) component in theECG is crucial for the study of AF.Some methods used to solve this problem are based on average beat subtraction (ABS). These methods are built onthe assumption that the AA is uncoupled with the VA. An average of the ventricular complexes (QRST complexes) isthen used to subtract VA [2]. Other approaches are blind source separation methods based on independent componentanalysis (ICA). They try to find independent components in an instantaneous linear mixture [3]. A major difficultyin ABS approaches is the limitation imposed by the use of a small number of VA average templates for general VAapproximation. In present ICA based approaches, a major gap is that only statistical priors are considered withouttaking into account the structural nature of signals. In order to circumvent these problems, a possible direction toexplore is the use of sparse source separation approaches based on signal adapted redundant dictionaries.During the past decade, important advances have been achieved in nonlinear signal approximation methods forsparse decompositions over redundant dictionaries (e.g. [4, 5, 6]). In many applications, these techniques offer betterperformances than those based on orthonormal transforms or direct time domain processing, thanks to their goodcapacity for efficient signal modeling. In this paper, we present a novel approach for VA and AA estimation. Weexplore a source separation method based on sparse decomposition of ECG signals on a redundant multi-componentdictionary. Such an approach allows for the consideration of priors on the structural nature of the different class ofsignals we are willing to separate. The multi-component dictionary is composed by functions specially designed to∗We thank Prof. Kappenberger (Lausanne University Hospital, Lausanne, Switzerland) for providing the clinical ECG database.†J. Molinero Hernandez was an exchange student from Universitat Polite`cnica de Catalunya (UPC), Barcelona, Spain. The workpresented in this paper was performed in the framework of his M.S. Thesis [1].match the main structural characteristics of VA and AA signals. We also present Weighted Orthogonal MatchingPursuit (Weigthed-OMP) [6] as a tool for generating ECG sparse approximations for source separation.2. UNDERDETERMINED SPARSE SOURCE SEPARATIONLet fmixj (t) : j ∈ [0,M − 1] be a set of M signal mixtures generated by the weighted superposition of N sourcesignals fi (t) : i ∈ [0, N − 1] such that:fmixj (t) =N−1∑i=0aj,i · fi (t) + nj(t), (1)where nj(t) represents some additive noise.Source separation is a classical problem in many fields like acoustics, radio or medical signal and image processing.Many applications exist where the retrieval of the different additive components, forming a set of complex signals, isrequired. The use of sparse signal representations in source separation problems was proposed in [7] in order to exploitsome prior knowledge about the structural characteristics of each fi (t). Many signals can be sparsely representedby the superposition of a limited number of atoms from an adapted dictionary of functions (D = {gl(t) : l ∈ Ω}):fi (t) =∑l∈Λbil · gil (t) +R|Λ|fi(t), (2)where bil are the atom coefficients, Λ ⊂ Ω and RΛfi(t) is an eventual residual that depends on Λ.Combining Eq. (1) and Eq. (2), one derives that, in order to solve the separation problem, the set of mixingcoefficients aj,i and a sparse set of expansion coefficients bil must be recovered while nj(t) are kept as small aspossible. Indeed, accurate sparse models can efficiently capture the structural nature of signals, leading to bettersource separation results as exposed in [7].A challenging form of sparse source separation problem is when there are fewer mixtures than sources. An exampleis the separation of the different components from a single ECG channel trace. This is studied in this work throughthe use of a novel sparse source separation approach. In order to do this, we adapt the two stage separation processproposed in [7] to the particular case of ECG components estimation: First, we a priori design an overcompletedictionary where sources are assumed to be sparsely representable. Second, the sources are unmixed by exploitingtheir sparse representability.3. ECG COMPONENTS SEPARATIONIn this section, we formulate the estimation of ECG activities according to the signal models and sparse sourceseparation strategy described in Sec. 2. The ECG signal (fECG) is modeled as a noisy mixture of the two cardiacactivities of interest (fAA and fV A):fECG = fAA + fV A + n, (3)where n stands for the noise.The generation of good sparse models for fAA and fV A requires the use of basis functions fitting their particularstructures. As shown in the following, fAA and fV A have quite different characteristics, and this is what enables usto separate them. The approach we propose is based on the decomposition of fECG on a redundant dictionary (D)composed by the union of two sub-dictionaries: DV A suited for representing the ventricular activity and DAA, betteradapted for representing the atrial activity. In the following D,DV A and DAA stand for the synthesis matrices ofD,DV A and DAA respectively, where each column represents an atom of the dictionary. Hence,fECG ' D · b = DAA · bAA +DV A · bV A. (4)Given the noisy nature of fECG and the high complexity of each of its components, we exclusively consider sparseapproximations in this work. According to Eq. (4), b is composed of two parts (bAA and bV A), each one containingthe coefficients related to DAA and DV A. The approach we explore is simple: one generates a sparse approximation offECG on D and then, estimates for fAA and fV A are reconstructed by just using the components from the appropriatedictionary:fAA ' DAA · bAA and fV A ' DV A · bV A.50 100 150 200 250 300−0.4−0.3−0.2−0.100.1temporal axisamplitudeOriginalApproximation0 50 100 150 200 250 300 350 400 450 50000.10.20.30.40.50.60.70.80.91temporal axisamplitudeβ=1β=1.5β=2β=2.5Fig. 1. Left: QRST VA complex and its approximation using just 3 atoms. Right: Effect of β on the GGF (seeEq. (5)).3.1. Modeling ECG Ventricular ActivityThe sub-dictionary DV A is generated by all possible translations of the Generalized Gaussian Function (GGF):gV A(t) = C1 exp(−( |t− p|α)β), (5)where C1 is a normalizing constant, α determines the scale and β the peakiness. This waveform allows to well approx-imate the structure of a VA complex using few atoms. Fig. 1 (left) shows a QRST complex and its approximation byusing just 3 atoms. With respect to the Gaussian function, fine tuning β gives us higher flexibility to approximateQ, R and S peaks (see Fig. 1 on the right).The possible values for α and β have been chosen experimentally after an extensive set of tests: α ∈ {3, 4, ..., 7}∪{49, 50, ..., 54}, the first set adapted for Q, R and S waves and the second intended for T wave approximation, whileβ ∈ {1.5, 1.6, ..., 2.2}. Together with p, this makes DV A highly coherent, but also very flexible for VA approximation.However, such dictionary is far from being optimal, and several improvements are still possible, mainly concerningthe approximation of T waves.3.2. Modeling ECG Atrial ActivityThe sub-dictionary DAA is generated by all possible translations of a real Gabor function:gAA(t) = C2 exp(−(t− pα√2)2)cos(2pik(t− p)N−∆ψ),where C2 is a normalizing constant, N is the signal length, α tunes the scale, k the frequency and ∆ψ the phase.This waveform is specially adapted for AA approximation. Indeed, as can be observed in Fig. 2 (left), fibrillatingAA is of oscillatory nature, which is a perfect fit for the optimal spatio-temporal frequency localization of Gaborfunctions (see Fig. 2 on the right).0 100 200 300 400 500 600 700−0.04−0.03−0.02−0.0100.010.020.030.040.050.06temporal axisamplitude0 50 100 150 200 250−0.15−0.1−0.0500.050.10.15temporal axisamplitudeFig. 2. Left: Example of a simulated AA wave during fibrillation. Right: Gabor atom.The values of the Gabor function parameters have been determined through an extensive analysis on several ECGsignals (see [1]). During the design of DAA, special care in limiting the maximum correlation between DAA and DV Aatoms was required. Indeed, an excessive correlation between this two sub-dictionaries translates into a completefailure of the algorithm.4. ECG SPARSE DECOMPOSITION BY WEIGHTED-OMPEq. (4) involves an underdetermined problem that has no unique solution. The search for the sparsest approximationof fECG requires the exhaustive testing of all coefficient possibilities, i.e. it is a combinatorial problem. Variousalternative approaches have been proposed in order to make the retrieval of a solution for b computationally affordable(see [4] for a review of some of them). These, in general, do not guarantee the recovery of the sparsest solution.Nevertheless, recent results show that under certain conditions on the dictionary and the signal, these sub-optimalmethods find the sparsest solution [5, 6]. In our case, Greedy Algorithms appear to be the most appropriate familyof decomposition algorithms to handle the very large dictionary we use.4.1. Weighted Orthogonal Matching PursuitWeighted-OMP [6] is the greedy algorithm used to decompose ECG signals. It is a weighted version of OrthogonalMatching Pursuit (OMP), that iteratively builds m-term approximants by selecting at each step the most appropriateterm from D according to a selection rule. Each iteration k : k ≥ 0 can be seen as a two-step procedure:1. A selection step where an atom glk ∈ D is chosen according to:glk = arg maxgl∈D|〈rk, gl〉| · wl,where wl ∈ [0, 1] is a pre-estimated weight that reflects, according to a predefined model, the a priori likelihoodthat atom gl may be a correct component of the signal.2. An orthogonal projection step where an approximant fk+1 ∈ span(gip : p ∈ {0, ..., k}), and a residual rk+1 =f − fk+1 are generated (notice that r0 = f , and rk+1 ⊥ fk+1 ∀k). This stage updates, at every step, the set ofscalar expansion coefficients.The signal representation generated by Weighted-OMP is, thus, of the form of Eq. (2).Weighted-OMP was recently proved in [6] to outperform OMP when using coherent dictionaries and reliable priorinformation. Weighted-OMP can consider, in the decomposition algorithm, a priori models about the behavior ofthe dictionary in use with the class of signals to decompose. The use of good enough a priori models can reduce theinstability of Weighted-OMP with respect to OMP when trying to recover sparse approximations/representations.4.2. Weight Generation: Relation Between ECGA Priori Knowledge and the DictionaryThanks to the structure of VA, fECG can be divided in VA periods. In addition, each VA period can be divided in aset of intervals corresponding to the different VA waves (Q, R, S and T) and an interval without ventricular activity.VA intervals can be estimated and identified in practice through the use of QRST point estimators (e.g. see [8], usedin this work). This prior information can thus be used to generate wl ∀l. The a priori knowledge obtained from [8]needs to be related with D in the following way.D is divided in DAA and DV A. Due to dynamics of AF, AA can be found through all the VA period. Hence, DAAatoms cannot be penalized. This is the reason why in this study we force: wl = 1 ∀l : gl ∈ DAA. To the contrary,the selection of gl ∈ DV A can be successfully influenced by the use of the available a priori information. DV A, asseen in Sec. 3, is composed of a block optimized for QRS waves (ventricular depolarization) and a block designed forT waves (ventricular repolarization). Depending on the VA interval, wl can be set to 1 for every gl ∈ DV A belongingto the appropriate kind for that interval. In case a gl is unsuitable for a given interval, wl can be set to a penalizingfactor 0 ≤ τ < 1. Thanks to the reliability of the estimators used in this work, it turned out that the best value forτ in our experiments is 0.5. EXPERIMENTAL RESULTS5.1. ValidationA biophysical computer model of the atria was used to obtain a realistic atrial electrical activity on the torso [9].The AF signals that were generated in the 12-lead ECG were added to a clinical 4-second standard 12-lead ECG ofan AF paroxysmal patient (78 years old) in sinus rhythm in which the P waves were removed. The clinical ECGwas selected to represent the VA in AF as closely as possible. The ratio between the power of the original signal(simulated AA) and the estimation error (estimated AA - simulated AA) was used to evaluate the performance ofour method.5.2. ResultsFirst of all, we want to underline that we validated our choice of Weighted-OMP instead of OMP with these simulatedmeasured 4-second ECG signals. By using Weighted-OMP, we increased the SNR in the recovery of VA (respectively,AA) by 0.81 dB (respectively, 0.65 dB). All the following results were obtained by approximating ECG signals with50 atoms. Fig. 3 shows the resulting separation of VA and AA for the simulated measured 4-second ECG signal onlead V1. One can see how our method succeeds in approximating each one of the VA episodes and in separating, atthe same time, the AA with surprising accuracy.1 mV1 sec(a)(b)(c)(d)(e)Fig. 3. (a) Simulated measured 4-second ECG signal on V1. (b) Original VA on V1. (c) Estimated VA on V1 (SNR: 8.69 dB). (d) Simulated AA on V1. (e) Estimated AA on V1 (SNR : 6.81 dB).In order to study the influence of the AA amplitude on the method, three different simulated AA signals werecreated; 50%, 100% and 150% of the original simulated AA amplitude. The ratio between the power of the originalactivity (VA or AA) and the error on the estimated one was evaluated on leads VR, V1 and V4 (see Table 1). Wecan observe that the quality of AA estimation depends on the lead and its original amplitude. The AA SNRs aremuch higher with the 150% original amplitude and the overall performance on lead V1 is better than those on theother two leads. Of course, the SNR values are directly related to the signal amplitude. In V1, the AA amplitude ishigher compared to other leads. However, we observe a decrease of VA estimation performance in lead V1.Fig. 4 shows the resulting separation of the VA and AA of the clinical 4-second ECG signal on lead V2. Apartfrom the visually satisfying component separation, the resulting signals were validated using estimated power spectraldensities (PSD). The dominant frequency of VA (respectively, AA) is between 1 and 2.5 Hz (respectively, between 3and 10 Hz). The fact that there is no presence of VA dominant frequencies in the AA estimated PSD demonstratesthe quality of our clinical results. Further results can be found in [1].0.5·AA+VA 1·AA+VA 1.5·AA+VAlead VRVA SNR(dB) 11.06 10.88 11.08AA SNR(dB) -6.94 -1.05 2.61lead V1VA SNR(dB) 11.13 8.69 2.41AA SNR(dB) 3.61 6.81 4.28lead V4VA SNR(dB) 12.33 11.94 11.66AA SNR(dB) -6.53 -0.8 2.4Table 1. Signal-to-noise ratio (dB) on lead VR, V1 and V4. Our method performance is tested on 3 different AAamplitudes (50, 100 and 150 % of the original simulated signal).X: 1.56 HzY: 0.0100 mV2/HzX: 1.56 HzY: 0.0100 mV2/HzX: 7.55 HzY: 0.0002 mV2/Hz1 mV1 sec(a)(b)(c)(d)2.5 mV1 sec(e)(f)(g)Fig. 4. (a) Clininal 4-second ECG signal on V2 with a dominant frequency of 1.56 Hz (see its PSD (e)). (b)Estimated VA on V2 with a dominant frequency of 1.56 Hz (see its PSD (f)). (c) Estimated AA on V2 with adominant frequency of 7.55 Hz (see its PSD (g)). (d) Estimated AA on V2 magnified 5 times.6. CONCLUSIONSIn this paper, we present a new framework based on sparse source separation that can be used for QRST cancellation.Results appear to be very promising. Additional works are planned for the improvement of the dictionary design,specially concerning T wave modeling. Moreover, more efficient use of a priori knowledge can be studied and thenumber of atoms used for decomposing the signals can be optimized.REFERENCES[1] J. Molinero Hernandez, “Sparse decompositions for ventricular and atrial activity separation,” M.S. thesis, SignalProcessing Institute, Ecole Polytechnique Fe´de´rale de Lausanne, Switzerland, August 2005.[2] M. Lemay, V. Jacquemet, A. Forclaz, J.-M. Vesin, and L. Kappenberger, “Spatiotemporal QRST cancellationmethod using separate QRS and T-waves templates,” in Computers in Cardiology 2005, September 2005.[3] F. Castells, J. Igual, J. J. Rieta, C. Sanchez, and J. Millet, “Atrial fibrillation analysis based on ICA includingstatistical and temporal source information,” in ICASSP, 2003, pp. V93–V96.[4] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition by basis pursuit,” SIAM J. Sci. Comp.,vol. 20, no. 1, pp. 33–61, 1999.[5] D. L. Donoho, M. Elad, and V. Temlyakov, “Stable recovery of sparse overcomplete representations in the presenceof noise,” To appear in IEEE Trans. Inform. Theory.[6] O. Divorra Escoda, L. Granai, and P. Vandergheynst, “On the use of a priori information for sparse signalapproximations,” To appear in IEEE Trans. Signal Processing.[7] M. Zibulevsky and B. A. Pearlmutter, “Blind source separation by sparse decomposition in a signal dictionary,”Neural Comp., vol. 13, no. 4, pp. 863–882, 2001.[8] R. M. Rangayyan, Biomedical Signal Analysis : A Case-Study Approach, John Wiley & Sons, Inc., 2002.[9] V. Jacquemet, M. Lemay, J.-M. Vesin, A. van Oosterom, and L. Kappenberger, “A biophysical model of ECGsignals during atrial fibrillation used to evaluate the performance of QRST cancellation algorithms,” in Computersin Cardiology 2005, September 2005.

http://infoscience.epfl.ch/record/87185/files/Divorra_Escoda2005_1401.pdf

Ventricular and Atrial Activity Estimation Through Sparse ECG Signal Decompositions

Abstract

Similar works

Full text

Available Versions

Infoscience - École polytechnique fédérale de Lausanne

Infoscience - École polytechnique fédérale de Lausanne