Blind Analysis of Atrial Fibrillation Electrograms: a Sparsity-aware Formulation

The problem of blind sparse analysis of electrogram (EGM) signals under atrial fibrillation (AF) conditions is considered in this paper. A mathematical model for the observed signals that takes into account the multiple foci typically appearing inside the heart during AF is firstly introduced. Then, a reconstruction model based on a fixed dictionary is developed and several alternatives for choosing the dictionary are discussed. In order to obtain a sparse solution, which takes into account the biological restrictions of the problem at the same time, the paper proposes using a Least Absolute Shrinkage and Selection Operator (LASSO) regularization followed by a post-processing stage that removes low amplitude coefficients violating the refractory period characteristic of cardiac cells. Finally, spectral analysis is performed on the clean activation sequence obtained from the sparse learning stage in order to estimate the number of latent foci and their frequencies. Simulations on synthetic signals and applications on real data are provided to validate the proposed approach.


Introduction
The clinical term atrial fibrillation (AF) refers to a family of common heart disorders characterized by fast and uncoordinated activations in the atrium.The mechanisms causing the initiation and maintenance of AF include a set of heterogenous interactions at different levels (cells, tissues and the whole heart) chang-ing along time and resulting into different states of AF [7,23].Several theories about the physiological causes underlying AF initiation and maintenance have been formulated over the last 50 years [21,22], and numerous efforts have been made in the implantation of algorithms for the automatic detection of AF states [29].One of the most prominent hypothesis considers multiple uncoordinated activation foci placed at different locations inside the atrium.These fast and asynchronous activations cause a disordered global electrical activity that contributes to AF maintenance [13].In contrast, during normal heart operation conditions (sinus rhythm) we observe a single activation focus, placed at the sinus node, acting as a pacemaker for the whole heart and leading to a regular global electrical activity.
In order to understand the pathophysiology of AF, dominant frequency analysis (DFA) has been traditionally used to analyze the data collected from electrocardiograms (ECGs) or electrograms (EGMs).DFA is useful for identifying the areas corresponding to the highest activation frequencies that may be the drivers maintaining AF, and therefore the targets of ablation therapy for AF termination [28].However, DFA provides very limited information about the signal's structure, since it is based on the implicit assumption that the underlying signal consists of a single quasi-periodic component plus an irregular component [6].Hence, the only spectral parameter required is the dominant frequency (DF), which characterizes the periodicity of the signal, but is very sensitive to distortions and can provide misleading information [24,25].
More recently, organization analysis techniques have been introduced, and additional parameters, such as the regularity index (RI) and the organization index (OI), have been used to describe the signals [6,8].Many other linear and non-linear measures have been proposed for the characterization of AF [16,18,26]: the cross-correlation index, the non-linear association measure, the fractionation index, etc.However, all of them are based on the same implicit assumption that the observed signals are modeled by a single regular component with an additional distortion and noise.
This paper presents a novel formulation which takes into account the presence of multiple foci and the sparse nature of EGM signals.This new formulation includes two major theoretical contributions.Firstly, it takes into account the multiple activation foci to perform spectral analysis, detecting the number of foci and their frequencies.Secondly, based on the sparse nature of the recorded signals, a sparsity-aware learning technique based on the Least Absolute Shrinkage and Selection Operator (LASSO) is applied, followed by an additional stage that gets rid of spikes that violate the biological restrictions, to obtain an activation sequence on which the spectral analysis is performed.
In this paper, bipolar intracardiac electrograms (EGMs), obtained placing a set of electrodes in direct contact with the heart muscle during heart surgery [24,28], are used for the analysis.The resulting algorithm applied to the signals consists of four steps: 1. Pre-processing to eliminate potential artifacts, especially outside of the desired frequency range.2. Inferring the spike trains associated to the activation times using a LASSO [33], followed by a post-processing step to ensure that biological restrictions are met [20].3. Performing a sparse spectral analysis of that activation sequence, using an iterative deflation approach to detect the number of foci and their frequencies.4. Applying a late post-processing step in order to eliminate harmonics and subharmonics.
A preliminary version of this paper has appeared in [20].When compared to [20], this contribution includes several theoretical and practical modifications: the pre-processing stage has been enhanced, the choice of the reconstruction dictionary is briefly discussed and a novel wavelet-based dictionary developed, a method for estimating the noise variance has been incorporated and the spectral analysis stage has also been improved.However, the most important extensions in the paper are in the application side: more realistic simulations on synthetic data have been performed using a typical activation shape for EGMs [8], and applications on real data have been performed for several patients under sinus rhythm and AF conditions.
The paper is structured as follows.Section 2 describes the prevalent approach for the analysis of EGMs: dominant frequency analysis.Section 3 details the problem formulation used throughout the paper, defines the novel mathematical model (based on a set of unobserved latent signals) proposed for describing the recorded EGMs, the sparsity-aware formulation introduced for solving it, and the dictionary considered for modeling the unknown latent signals.Section 4 describes the approach proposed for inferring the sparse activations.Then, Section 5 shows how the sparse spike train inferred using this approach can be used to perform sparse spectral analysis (SSA), thus inferring the number of latent foci as well as their activation frequencies.Results are presented in Section 6, both for synthetic data and for real EGMs, obtained both under sinus rhythm and atrial fibrillation (AF) conditions.Finally, Section 7 provides conclusions and perspectives.

Background: Dominant Frequency Analysis
Dominant frequency analysis (DFA) is the prevalent approach for the analysis of EGMs.DFA assumes implicitly that the observed signals are composed of a single regular component (i.e., a quasi-periodic signal) plus an irregular component that contains the remaining noise and distortion.Hence, from a mathematical Algorithm 1 Dominant frequency analysis (DFA) for the q-th signal [8]. 1. Band-Pass filtering from 30 Hz to 400 Hz. 2. Rectification of the resulting signal, recovering near direct current (DC) spectral components.3. Low-Pass filtering with a 15 Hz cut-off frequency.4. Computation of the spectrum using a localized Fast Fourier Transform (FFT) with a Hanning window of Λ = 4 s duration, resulting in a resolution f Λ = 1/Λ = 0.25 Hz in the frequency domain.5. Search for the peak with the maximum amplitude in the frequency domain.The frequency associated to this peak is the dominant frequency (DF) of the q-th EGM, fq .point of view, the q-th EGM, 1 ≤ q ≤ Q with Q denoting the number of outputs, can be modeled as [8] y where φ q (t) idenotes the average shape of the regular component of the signal, T q its period, τq the delay for k = 0, and w q (t) is used to represent the irregular components.The goal of DFA is characterizing that quasi-periodic signal through its average period, T q , or equivalently its average frequency, fq = 1/ Tq , which is the so called dominant frequency (DF).Occasionally other parameters, such as the organization or the regularity indexes, are also obtained to determine whether the estimated DF is reliable or not [6,8].
The DF is usually obtained separately for each channel using standard spectral analysis techniques.A typical signal processing approach includes the five steps for each EGM shown in Algorithm 1 [8].Several segments can be averaged in order to improve the estimation of the dominant frequency.However, the ability of the DF to reflect the average atrial activation rate depends on the accuracy of Eq. (1) in representing the true observed signal.Unfortunately, several characteristics of atrial activation, such as the complexity of the electrogram morphology, can alter the power spectrum.In these cases, fq is more related to the complexity of the signal than to the atrial activation rate, thus providing misleading information [24,25].

Problem Formulation
In this section a novel formulation is proposed as an alternative to the standard DFA formulation of Eq. ( 1) First of all, a more realistic mathematical model based on the assumption that the observed signals are the result of several unobserved latent functions (the unknown activation foci that we want to estimate) propagating through the heart is introduced.Then, since the real shapes of these latent signals are not precisely known, a sparsity-aware formulation, based on an overcomplete dictionary, is introduced.

Signal Model
Assume that the recorded EGMs are composed of the sum of several periodic or quasi-periodic signals plus distortion and noise.Each of these observed periodic signals stem from a set of sparse activation foci (spike trains) that propagate through the atrium and reach the sensors.Hence, these unobserved activations play the role of latent signals, providing a way of describing the correlation between the outputs.The primary goal of this paper is detecting the number of activation foci, as well as their frequencies.
From a mathematical point of view, we consider a model with Q correlated outputs, y q (t), obtained from a set of bipolar electrodes.Furthermore, we assume that these observations are generated by R unobserved activation foci (latent signals) propagating inside the atrium, plus noise and interference.Hence, the output of the q-th channel (1 ≤ q ≤ Q) is modeled as where p r (t) (1 ≤ r ≤ R) denotes the r-th focus, w q (t) models all the elements in the q-th output that cannot be explained by the model (i.e., noise, interferences and distortion), h r,q (t) is the impulse response of the channel between the r-th focus and the q-th output EGM and * denotes the standard linear convolution operator.Note that h r,q (t) includes the response of the sensor and can be slowly time-varying.However, since the sparse learning and the subsequent spectral analysis are performed using short time windows, the channel can be considered time-invariant in practice.Fig. 1 shows the underlying signal model corresponding to Eq. ( 2).Since we are not interested in recovering the precise shape of the activations, but only in their num- ber and frequencies, we model them as periodic spike trains (this is not a limitation however, since the shape of the activations can be included inside the channel's impulse response, h r,q (t)): with δ(t) denoting Dirac's delta, T r = 1/f r the average period of the r-th spike train (with f r denoting its associated average frequency) and τ r its shift w.r.t. the origin (0 ≤ τ r < T r ).Substituting Eq. (3) into Eq.( 2), the q-th output becomes The discrete-time version of this model, obtained assuming a uniform sampling frequency, f s = 1/T s , is where h r,q [n, k] = h r,q (nT s −kT r −τ r ) is the discretetime equivalent channel and w q [n] = w q (nT s ) are the noise plus distortion and interference samples at the sampling instants.In the sequel, this discrete-time model is used for inferring the global spike train (i.e., the spike train resulting from the sum of the R foci) and estimating R and f r = 1/T r for r = 1, . . ., R.

Reconstruction Model
Let us denote the (N + 1) × 1 vector with the samples from the q-th EGM as y q = [y q [0], . . ., y q [N ]] .Now, let us define the N × 1 vector containing the discrete-time differentiation of the q-th output, Since the precise shape of the activations is irrelevant in this case, and the number of latent foci is still unknown, z q [n] is approximated by a mixture of shifted smooth generic curves, i.e., where ε q [n] is additive white Gaussian noise (AWGN) with zero-mean and unit variance (i.e., ε q [n] ∼ N (0, 1)), σ q denotes the actual noise variance, assumed to be known or estimated from the data, and β m,q [n] are the coefficients of the q-th output associated to the m-th sampled activation shape, φ m [n] = φ m (nT s ) for 1 ≤ m ≤ M and 1 ≤ q ≤ Q. Eq. ( 5) is the discrete-time version of the continuous-time model Indeed, the reconstruction models assumed by the sparsity-aware formulation, given by Eqs. ( 5) and ( 6), are very similar to the assumed underlying signal models, given by Eqs. ( 2) and ( 4), although there are two important differences: 1. Focusing on the discrete-time models, Eq. ( 5) describes the first-order time-difference of the sampled EGM signals, whereas Eq. ( 4) models the EGM signals themselves.Regarding the equivalent continuous-time models, this corresponds to working with the first derivative of the signals instead of the signals themselves, which is a common approach to remove their baseline.
2. Since the number of activations and their shapes (i.e., the impulse responses associated to the Q channels) are unknown, in the reconstruction model a set of M ≥ R activations, constructed using generic smooth curves φ m (t), is used.The same activation shapes are used for all the channels, and the relevant activations for each case are selected automatically through a sparse learning process based on LASSO with a post-processing stage.This allows us to effectively remove the subindex q from the original activations, h r,q (t), moving it to the set of coefficients, β m,q (t).

Ideal Dictionary
The optimum reconstruction model would be composed of a set of Q dictionaries tailored to the characteristics of each of the Q outputs.More specifically, the dictionary for the q-th output would be given by the following N × RN matrix, where Φ n,q is the N × R matrix composed of the R impulse responses between all the latent signals and the q-th observation at the n-th time instant, i.e., Using this dictionary, Eq. ( 5) can be expressed as where β q is an RN × 1 column vector composed of N subvectors of size R: with ] .This dictionary leads to the sparsest possible solution, containing at most P = N R r=1 T s /T r non-zero elements for each channel (P Q overall) that match the amplitudes of the activations, A r [k].In fact, as some activations will be masked due to the refractory period characteristic of cardiac cells, the number of non-zero coefficients will actually be lesser than P .

Alternative Dictionaries
Unfortunately, the ideal dictionary requires either knowledge of the h r,q (t) or a reliable estimation, something which may be unattainable for the application considered.However, it provides us with a criterion for comparing different dictionaries: the best dictionary will be the one that attains the sparsest representation, while obtaining a good reconstruction error.
In practice, the activation shapes are constructed using a family of smooth functions centered around the origin, φ m (t) with 1 ≤ m ≤ M and M ≥ R, such that their support is approximately −T m ≤ t ≤ T m (i.e., φ m (t) ≈ 0 for |t| > T m ).Furthermore, this family is selected in such a way that the time dispersion of φ m (t) increases with m (i.e., T 1 < T 2 < . . .< T M ).The discrete-time sequences used to build the dictionary are then obtained through uniform sampling of the truncated and time-shifted continuous-time functions, i.e., where T m /T s being the integer part of the ratio T m /T s .Hence, all the discrete-time activation elements suffer a delay of N M samples (i.e., N M T s seconds) that must be taken into account when interpreting the results obtained in the sequel. 1hus, the sparse model in Eq. ( 5) can be expressed more compactly in matrix form by defining a set of Concatenating all these matrices, we obtain the following global N × M N dictionary: Since Φ is a non-square matrix with more columns than rows, this dictionary is said to be overcomplete.
Hence, many elements in Φ k (1 ≤ k ≤ N ), and thus also in Φ, will actually be zero.Therefore, Eq. ( 5) can be expressed by means of Eq. ( 13) as where ] is an N × 1 column vector with the noise samples associated to each sample of z q [n], and β q is an M N × 1 column vector composed of N subvectors of size M with a structure identical to the one shown in Eq. ( 10), but In [20], the activation shapes were obtained sampling truncated and time-shifted Gaussian functions,2 As a better alternative, in this paper we propose using the Mexican hat wavelet, also known as Ricker wavelet, which is the negative normalized second derivative of a Gaussian function.Wavelets are frequently used in bioinformatics, biometrics and image applications [30,32] due to their ability to represent data compactly [11,19,27] and to estimate precisely details of ECG signals [18].Furthermore, thresholding methods can be easily derived to enhance the sparsity of the data representation [3] and characterize biomedical data [2,9].In particular, the mexican hat wavelet is often used in electrocardiographic applications due to its similarity to activations observed in real data [1,5,10,14].Fig. 2 displays the wavelet-based overcomplete dictionary used in the simulations, constructed using M = 10 and σ m = 0.4m for m = 1, . . ., M .

Proposed Sparse Solution: post-processed LASSO
A procedure for obtaining a sparse solution β q for Eq. ( 14), which fulfills the biological constraints is developed in this section.The proposed method, which is a variation of the algorithm introduced in [34], follows a two step approach: an initial sparse solution obtained by applying a LASSO regularization is followed by a greedy procedure for selecting only the largest coefficients that respect the biological constraints.
In order to obtain a sparse regressor, from which the information on the arrival times can be retrieved, β q is firstly estimated by means of LASSO [33].Namely, βL1 q is given by βL1 q = arg min ]] , β q 1 denotes the L 1 norm of β q and λ q is the regularization parameter, that indicates the trade-off between sparsity and estimation precision: the higher the value of λ q the more emphasis will be placed on obtaining a sparse solution, although at the expense of an increased quadratic error in the approximation.
However, in order to obtain an even sparser representation that takes into account the physiological restrictions imposed on the signals, an additional step after the computation of βL1 q is introduced.The samples associated to the arrival times of the spikes are estimated recursively as follows: nk,q = arg max n=1,...,N where I(L(x)) is an indicator function, i.e., a function that takes a value equal to one if the logical condition L(x) is fulfilled and zero otherwise, with the logical condition in this case being and η q and N min user-defined thresholds.The first one, η q , is used to discard the βL1 q [n] with a small L 1 norm, which contribute to improve the signal reconstruction but provide no information on spike localization.A value η q = 3σ q is used, as it has been found out empirically to provide good results.The second one, N min , accounts for the fact that consecutive pulses cannot overlap.Thus, in practice N min is chosen in such a way that N min /f s ≈ T ref , with T ref being the minimum possible value of the refractory period: T ref = 100 ms for sinus rhythm and T ref = 50 ms for AF.
The final procedure used in practice to implement Eq. ( 18) for the q-th channel is an iterative greedy approach that follows the steps shown in Algorithm 2. Following this procedure, a set of P arrival times are obtained and used to construct an activation sequence (also called spike train) composed of Kronecker deltas at the locations of the activations,3 This clean sequence is used in the sequel to perform the spectral analysis, as it is free from the effect of Algorithm 2 Iterative greedy approach for selection of the final spikes for the q-th signal.Initialization: set k = 1 and β max = ∞.while β max > η q do Select the index corresponding to the largest coefficient: nk,q = arg max n=1,...,N Set the unknown channels, h r,q (t), the particular dictionary used, given by φ m (t), and the noise, perturbations and interferences, represented by w q [n].
As an example of the behavior of the proposed spike detection procedure, Figs. 3, 4 and 5 show representative EGM signals and the corresponding detected activation sequences for three cases: a synthetic signal generated following the procedure described in Section 6; a real signal recorded under sinus rhythm condition; and a real AF signal.In all cases the detected spikes correspond to true activations, although some activations seem to be missed for the real AF signal.

Deflation Approach for Peak Estimation
The spectral analysis is based on applying an iterative deflation approach to the discrete Fourier transform (DFT) of π q [n], extracting peaks with decreasing amplitudes up to a user defined threshold.Hence, since the spectral analysis is applied to the inferred sparse activation sequence, this approach is called sparse spectral analysis (SSA) [20].The number of peaks extracted (after the post-processing described in Section  5.2) is an estimate of the number of existing foci, and their locations an estimate of their frequencies. 4he first step in the SSA algorithm consists in splitting π q [n] into J windows containing N s samples with or without overlap, as done in DFA (see Section 2).Then, after bandpass filtering to focus on the desired frequency band, Algorithm 3 is applied to the DFT of each windowed segment, Π j q (f ) = F π j q [n]w[n] for 1 ≤ j ≤ J, with w[n] denoting an appropriate windowing function (the Hamming window in this case).Algorithm 3 follows a deflation approach, searching iteratively for the highest peak of |Π j q (f )| within the Algorithm 3 Iterative Spectral Analysis for the q-th signal.
for j = 1, . . ., J do 3. Filter the signal: end while end for frequency range that is physiologically interpretable (0.5 ≤ f r ≤ 2 Hz for sinus rythm and 2 ≤ f r ≤ 10 Hz for AF) and adding it to the set of potential activation frequencies, f j q .After each iteration a secondorder IIR digital notch filter is applied to the signal centered around the detected frequency with bandwith B 3 dB = 2f Λ = 0.5 Hz to remove the detected peak before searching for a new one.The algorithm stops when the highest peak detected is below a threshold, Γ j q = γ q max |Π j q (f )|, with γ q user defined.Figs. 6, 7 and 8 show the spectrum obtained iteratively for the single segments in Figs.3-5.All these figures display the spectrum obtained after each iteration of Algorithm 3 (i.e., i = 1, . . ., R).

Post-Processing: Discarding Harmonics
The post-processing stage takes the set of potential activation frequencies detected inside each window, f j q , and determines whether they belong to different activation foci or not applying the following steps: 1. Elimination of repeated frequencies.Two frequencies, f 1 and f 2 , correspond to the same fo- If this happens, the one associated to the smallest peak is deleted.2. Analysis of 2/3 frequency relationships.Due to the frequency range used in the analysis, given a single frequency, f 0 , in practice at most two harmonics can be found: f 1 = 2f 0 and f 2 = 3f 0 .Thus, if the first and second harmonics of a given frequency, f 0 , have been detected, their relationship will be f 1 = 2 3 f 2 .This relationship is checked here, keeping only the frequency associated to the highest amplitude in the spectrum.3. Discovery of harmonics and subharmonics.If two detected frequencies have a harmonic or subharmonic relationship, only the one detected first in the spectral analysis is kept, deleting the other.4. Discovery of cross-modulation frequencies.Analyze whether each new element in fq is a crossmodulation product of two previously detected frequencies, i.e. whether f 3 = mf 1 + nf 2 for any two integers m and n.In this case, f 3 will be deleted.Given the frequency range of interest, and assuming f 2 > f 1 , only f 2 ± f 1 has to be checked (i.e., m = ±1 and n = 1).
With this analysis, the number of activation foci present in the EGMs, Rj q , as well as their frequencies, f j q , can be estimated.Continuing with the examples shown in Figs.6-8, the post-processing produces the following results: in Fig. 6, R = 2 with f1 = 8.18 Hz and f2 = 3.20 Hz for AF synthetic data (note that the second estimated frequency corresponds to first true one, whereas the first estimated frequency is the first harmonic of the second true one, which has a higher amplitude than the true frequency in this case); in Fig. 7, R = 1 with f1 = 1.67 Hz for sinusal real data (expected, since there should be a single focus for sinus rhythm with a frequency in the range 0.5-2 Hz); and in Fig. 8, R = 4 with f = [2.30, 7.25, 8.68, 3.28] Hz for AF real data (note that 5 spurious frequencies have been discarded).

Results
This section shows the results obtained after applying the spike extraction and spectral estimation algorithms described in the previous sections.We first provide a summary of the overall algorithm, including the important pre-processing stage, and describe the experiments on synthetic data, which are used to validate the proposed approach and expose the limitations of spectral estimation.Results on real data obtained both under sinus rhythm and AF conditions are then presented.

Overall Algorithm
Algorithm 4 summarizes all the steps for the analysis of both the synthetic and the real data.Firstly, a preprocessing stage is performed.Each EGM, obtained through uniform sampling with f s = 977 Hz, is decimated by a factor L = 2 (AF signals) or L = 4 (sinus rhythm) after filtering using a 12 tap low-pass digital filter with cut-off frequency ω c = π/L rad. to avoid aliasing.The decimation factors used can be seen as a compromise between performance (increasing L leads to important activations not being detected) and computational cost.This is followed by a discrete-time differentiation for baseline removal, as detailed in Section 4. After the pre-processing stage, Algorithm 2 is applied on a channel-by-channel basis in order to obtain a set of Q clean activation sequences.The sparsity parameter, λ q , is set to an appropriate value obtained after extensive simulations with synthetic data, 5  whereas σ q is assumed to be known for the synthetic data and estimated from the EGMs (see Section 6.3) 5 The results obtained are not too sensitive to the value of λq as long as it remains within a proper range, i.e., as long as it is not too small (thus resulting in a non-sparse least squares solution) or too large (leading to small activations not being detected).See Table 2.

Spike extraction:
(a) Set λ q .(b) Set σ q or estimate it from the data as shown in Section 6.3.(c) For each segment of each channel: obtain βL1 q using LASSO (as described in Section 4) and apply Algorithm 2 to discard invalid activations.
3. Discard invalid channels as shown in Section 6.3.

Spectral analysis of the clean activation sequence
given by Eq. ( 22): (a) Divide the signal into segments composed of N samples with or without overlap.(b) Apply Algorithm 3 to each segment of each valid channel.(c) Perform the post-processing described in Section 5.2 to discard spurious frequencies.
for the real data.Invalid channels (i.e., channels with too few or too many activations) are discarded following the procedure described in Section 6.3.Finally, the spectral analysis summarized in Algorithm 3 (followed by the post-processing described in Section 5.2) is applied to each segment of each valid channel's activation sequence.As a result, we obtain Rj q and f j q for q = 1, . . ., Q and j = 1, . . ., J, where Q ≤ Q is the number of valid channels.

Synthetic Data
The synthetic data are generated following the approach shown in Algorithm 5. First of all, an activation sequence is constructed for each of the R latent foci (R ∈ {1, 2, 3, 4}), whose frequencies, f r for 1 ≤ r ≤ R are selected uniformly within the physiologically meaningful range: f r ∈ [0.5, 2] Hz for sinus rhythm (i.e., R = 1) and f r ∈ [2, 10] Hz for AF (i.e., R > 1).Their phases τ r are also selected uniformly within [0, 1/f r ].Furthermore, for R > 1 the newly generated frequencies (f 2 , . . ., f r ) cannot be equal to previous ones or have any harmonic or cross-product relationship with them, since this is unlikely to occur Algorithm 5 Generation of synthetic EGMs.
if sinus rhythm then Set R = 1, f min = 0.5 Hz, Discard activations violating the refractory period, obtaining taken from [8]. in real situations.A global activation sequence is obtained by concatenating all the activations from the different foci, and those activations placed within a time distance less than T ref (T ref = 100 ms for R = 1 and T ref = 50 ms for R > 1) w.r.t. a prior activation are eliminated, since they would be masked by the refractory period characteristic of cardiac cells in real situations.Finally, the synthetic EGMs are obtained convolving the pruned global activation sequence with the typical activation shape given in [8] with a length equal to 32 samples (i.e., time duration ≈ 32.75 ms for f s = 977 Hz).An example of an AF synthetic EGM with R = 2, and the activations detected in that case, can be seen in Figure 3.
Algorithm 4 is then applied to the synthetic data.In a first stage, the spike extraction procedure described in Section 4 is evaluated, focusing on the detection probability.Table 1 shows the detection probability, P d , as a function of the signal to noise ratio (SNR) for R ∈ {1, 2, 3, 4} using the optimum values of the sparsity parameter: λ = 5 • 10 −5 for sinus rhythm (i.e., R = 1) and λ = 5•10 −4 for AF (i.e., R > 1).The good behaviour of the proposed algorithm can be clearly appreciated, with detection probabilities above 90 % for SNR ≥ 30 dB when R = 1 and for SNR ≥ 25 dB when R > 1. Table 2 displays P d as a function of λ for R ∈ {1, 2, 3, 4} when SNR = 30 dB, showing that the results are not very sensitive to λ as long as it lies within a proper range.In a second stage, the performance of the iterative deflation approach for spectral analysis is evaluated.Here, the goal is detecting the correct number of foci and the results are presented in terms of the false alarm probability, P f a = Pr{ R > R}, and the non-detection probability, P nd = Pr{ R < R}.Fig. 9 displays the results for sinusal synthetic data (i.e., R = 1) as a function of the SNR and the relative threshold γ, showing that P f a is largely independent of the SNR (as long as SNR ≥ 25 dB) and can be reduced to levels below 0.1 by selecting γ ≥ 0.6.Note that, since at least one peak (the dominant frequency) is always detected, P nd = 0 when R = 1.Regarding AF synthetic data (i.e., R > 1), Fig. 10 displays ROC-type curves (i.e., 1 − P nd vs. P f a ) for R ∈ {2, 3, 4} and SNR = 30 dB. 6 The good performance in all cases can be clearly appreciated.Note that these results are homogeneous with the simpler simulations performed in [20], where artificial spike trains (i.e., without an activation shape or masking due to the refractory period) and a fixed loss probability for the spikes, P loss , were used.The probability of a spike being masked due to the refractory period of cardiac cells is P loss ≈ 0.28 for R = 2, P loss ≈ 0.43 for R = 3 and P loss ≈ 0.51 for R = 4.

Real Data
In this section the proposed approach is tested on real data from 5 different patients for both basal rhythms: sinus rhythm (induced by a pacing device) and AF (paroxysmal type).The study is composed of 24 signal sets, each one containing 9-10 channels recorded during ablation surgery inside the Right Atrium, using a lasso type sensor.Table 3 summarizes all the relevant information about the signals in the sudy.In all cases the measurements have been taken in several points of the Right Atrium close to the Pulmonary Veins, where the reentry circuits are more common and ablation therapy usually takes place.
Algorithm 4 is applied to all the channels in the signal set on a channel-by-channel basis.The parameters used (obtained from the extensive study performed on the synthetic data) for sinus rhythm are λ q = • 10 −5 , T ref = 100 ms and γ q = 0.6; whereas for AF they are λ q = 5 • 10 −4 , T ref = 50 ms and γ q = 0.6. 7The noise variance is estimated from the data following a simple procedure: construct a histogram of the EGMs amplitudes and estimate σ q as the standard deviation of the equivalent zero-mean Gaussian that fits the lowamplitude levels of the histogram.
Furthermore, thresholds are established on the minimum and maximum number of activations inside a segment that can be expected from a physiological point of view.If those thresholds are exceeded the signal is discarded as invalid, since it can provide misleading results.This may happen when a sensor is not in direct contact with the heart or due to distortions produced by the movement of the measurement (lasso) catheter, as it contacts the ablation catheter (a very common situation during ablation surgery).Fig. 11 shows an example of this type of distortions both for sinus and AF rhythms.
Examples of the signals and the detected spikes can be seen in Figs. 4 and 5 for sinus rhythm and AF respectively, and the corresponding single segment spectra are shown in Figs.7 and 8.In the sequel, two examples of the spectrograms obtained for sinus rhythm and AF are shown in Figs. 12 and 13 respectively.On the one hand, Fig. 12 shows the regular nature of sinus rhythm, with a single focus and a stable frequency.On the other hand, Fig. 13 shows the fragmented nature of AF EGMs, with a variable number of foci per seg- Finally, a summary of all the information extracted from the real data is provided in Table 4.This table shows, for each signal of each patient, the number of useful channels, Q, the number of segments available, J, the average number of foci detected and the standard deviation, R ± σ R , and the vector with the average detected frequencies, f , for f1 , . . ., f R with R denoting the result of rounding R to the closest integer number.First of all, from Table 4 it can be seen that many channels are labelled as not valid, since channels where a single segment is wrong are discarded in order to avoid erroneous conclusions.Concentrating on the valid patients, the algorithm provides the expected results for sinus rhythm: it usually detects a single focus (with two foci detected only occasionally in some segments of some channels) with a fre-    the results are harder to interpret, since there are many fluctuations in R and f across different channels and even segments, as shown in Fig. 13.However, looking at the global picture shown by Table 4 provides some interesting insights for the three patients with valid data.For instance, patient 1 shows a large num-ber of foci, R = 3, with high activation frequencies, fr ∈ [6.662, 8.005] Hz.In contrast, for patient 2 we have R = 1 with f1 = 5.9 Hz on average, a much lower activation frequency.Finally, patient 4 has R = 2 with fr ∈ [5.766, 6.983], which corresponds to an intermediate situation.

Conclusions and Perspectives
A completely novel signal processing methodology for the blind analysis of atrial fibrillation (AF) electrograms (EGMs) has been presented in this paper.The proposed approach is based on a new and more realistic mathematical model for EGM signals that takes into account the multiple unobserved latent foci (i.e., quasiperiodic spike trains) typically present during AF.In order to infer the number of foci and their frequencies, a sparse reconstruction model (using a waveletbased overcomplete dictionary) has been introduced.The reconstruction problem is then solved using an initial LASSO regularization followed by a second postprocessing stage which enforces the biological restric-tions imposed by the refractory period of cardiac cells.The algorithm has been validated through a set of elaborate experiments with synthetic data, as well as a preliminary study using 24 real data sets, both under sinus rhythm and AF conditions.Future work will focus on improving the spectral analysis (taking into account lost activations due to the masking effect caused by the refractory period and introducing eigen-value-based methods) and extending the proposed framework to the multi-channel case using some type of Group LASSO formulation.

Fig. 7 .
Fig. 7. Example of the SSA for a single segment of real data under sinus rhythm conditions with γ = 0.4: R = 1 and f1 = 1.67 Hz.

Algorithm 4
Blind analysis of EGMs: overview of the full sparsity-aware estimation algorithm.1. Pre-processing stage: (a) Low-pass filtering and decimation by a factor L = 2 (AF signals) or L = 4 (sinus rhythm).(b) Discrete-time differentiation.

Fig. 12 .
Fig. 12. Example of a spectrogram for sinusal real data with detected frequency in black (signal 5, channel 4).

Fig. 13 .
Fig. 13.Example of a spectrogram for AF real data with detected frequencies in black (signal 14, channel 1).

Table 3
Information for Real Data signals with Sinus and AF Rhythm.RIPV = Right Inferior Pulmonary Vein.LIPV = Left Inferior Pulmonary Vein.RSPV = Right Superior Pulmonary Vein.LSPV = Left Superior Pulmonary Vein.LPVA = Left Pulmonary Veins Atrium (the place where both Superior and Inferior Left Pulmonary Veins come together, where ablation isolation is generaly performed).

Table 4
Results obtained after processing all the real data.