Blind Source Separation in Persistent Atrial Fibrillation Electrocardiograms Using Block-Term Tensor Decomposition With Löwner Constraints

The estimation of the atrial activity (AA) signal from electrocardiogram (ECG) recordings is an important step in the noninvasive analysis of atrial fibrillation (AF), the most common sustained arrhythmia encountered in clinical practice. This problem admits a blind source separation (BSS) formulation that has been recently posed as a tensor factorization, using the Hankel-based block term decomposition (BTD), which is particularly well-suited to the estimation of exponential models like AA during AF. However, persistent forms of AF are characterized by short R-R intervals and very disorganized (or weak) AA, making it difficult to model AA directly and perform its successful extraction through Hankel-BTD. To overcome this drawback, the present work proposes a tensor approach to estimate QRS complexes and subtract them from the ECG, resulting in a signal that, ideally, only contains the AA component. Such an approach tackles the problem of blind separation of rational functions, which models QRS complexes explicitly. The data tensor admitting a BTD is built from Löwner matrices generated from each lead of the observed ECG. To this end, this paper formulates a variant of the recently proposed constrained alternating group lasso (CAGL) algorithm that imposes Löwner structure on the decomposition blocks. This is done by performing an orthogonal projection, which we explicitly derive, at each iteration of CAGL. Results from experiments with synthetic data show the consistency of the proposed Löwner-constrained AGL (LCAGL) in extracting the desired sources. Experimental results obtained on a population of 20 patients suffering from persistent AF show that the proposed variant outperforms other tensor-based methods in terms of atrial signal estimation quality from ECG records as short as a single heartbeat.


I. INTRODUCTION
A RRHYTHMIAS are heart diseases characterized by an irregular rate and/or rhythm of the heartbeats. Being the most frequent sustained arrhythmia diagnosed in clinical practice, decreasing life quality and increasing healthcare costs, atrial fibrillation (AF) is a supraventricular tachyarrhythmia distinguished by an uncoordinated and irregular atrial depolarization [1]. Persistent AF represents a particularly complex form of this arrhythmia, where extensive atrial remodelling has taken place due to sustained AF, significantly affecting atrial activity (AA) and causing AF perpetuation. This cardiac rhythm disturbance is considered an economical burden, as a patient suffering from AF costs, anually, approximately $8 700 more to healthcare providers than a patient without AF in the USA. It is also associated with a higher morbidity, representing a major health concern, since about half million hospitalizations with AF as the primary diagnosis are registered every year, also in USA [1]. If the number of new diagnosis stays stable, it is estimated that AF will increase from 2.3 to at least 12.1 million patients in the USA by 2050 [2], and from 8.8 to about 17.9 million in the European Union by 2060 [3], becoming then a new epidemic. Worldwide, the number of AF patients was estimated to be around 33.5 million in 2010 [4]. Considering these worrisome statistics and the fact that the electrophysiological mechanisms of AF are not totally understood, it is not surprising that this challenging cardiac condition has been a topic of intense research in the past few years, and is expected to be investigated even more deeply.
During AF episodes, the atrial depolarization, represented by the P wave during normal sinus rhythm, is replaced by fibrillatory waves, or f waves, which are present continuously throughout the electrocardiogram (ECG) recording. The f waves are masked by the QRS-T complex that corresponds to the ventricular depolarization and repolarization, i.e., the ventricular activity (VA), in each heartbeat. Fig. 1 illustrates an example of a persistent AF ECG recording with short R-R intervals and fine AF, i.e., when the AA has a very low amplitude. These two characteristics are quite common in persistent AF ECGs.
In order to provide a better understanding of the complex mechanisms behind AF, a detailed analysis of the AA signal is necessary. A cost-effective noninvasive method to perform this task is to extract the AA from the standard 12-lead ECG. During AF, the AA and the VA signals are typically assumed uncoupled, which makes the extraction of AA from the ECG admit a blind source separation (BSS) formulation [5]. The block-term decomposition (BTD) built from Hankel matrices, proposed as a technique to solve BSS problems in [6], was used to noninvasively extract the AA signal from AF ECG recordings, and shown to outperform the matrix-based techniques in this particular application [8]- [11]. Also, tensor-based techniques do not require orthogonality or mutual independence between the sources, enjoying essential uniqueness under mild constraints, and can perform BSS from rather short signal records, on a beat-to-beat basis.
The BTD tensor model is flexible, as it can accommodate sources with different characteristics, i.e., with diferent number of poles (or degrees, in the case of rational functions). Two commonly used methods for computing this decomposition are the nonlinear least squares (NLS) method implemented in the Tensorlab MATLAB Toolbox [12] and the alternating least squares method with enhanced line search (ALS-ELS) [13]. However, their performance depends considerably on their initialization and model structure, i.e., the rank of the matrix factors and the number of blocks of the tensor.
To overcome these limitations, an algorithm to compute the BTD called alternating group lasso (AGL) and a variant called constrained AGL (CAGL) that deals with linear constraints over block matrices were recently proposed [14]. The principle of AGL is to compute BTD as the minimization of a least squares (LS) fit term plus a regularization term in order to penalize the rank of the blocks. AGL is able to find an approximate BTD with weak assumptions on its model parameters (number of block terms and multilinear ranks), since only an upper bound is imposed on these parameters. In CAGL, approximate projections are performed at each iteration to ensure that the block matrices have low rank and belong to a specified subspace, a problem known as structured low-rank approximation (SLRA) [15]. Indeed, the estimation of low-rank block matrices is a crucial condition to extract the desired information contained in the (unique) model parameters. CAGL was also assessed as an AA extraction tool in [14], showing its ability to extract more information of the AF ECG signal than other existing techniques, as this method could extract two possible atrial source signals, probably having their origin in different parts of the atria.
The Hankel-BTD has proven to be a useful and powerful AA extraction tool in AF analysis, but a satisfactory performance is only achieved for segments with long R-R intervals and with well defined AA, visible in most of the segment. However, as stated earlier in this section, recordings with short R-R intervals (< 0.75 s) and fine AF (amplitude of the f waves lower than 0.1 mV) are quite more common during persistent AF episodes. When assessed in these challenging cases, the Hankel-BTD as well as the matrix-based techniques do not provide satisfactory results [16]. Automatically selecting the AA signal among the estimated sources after performing BSS is also an unsolved issue, since no systematic method is reported in the literature for this purpose.
Aiming to avoid such limitations in these common and challenging persistent AF scenarios, the present work proposes the BTD built from Löwner matrices as a solution for BSS of rational functions [17] to model the VA and separate it from the AA. This strategy suits the characteristics of VA in ECG recordings, since the QRS complex can be well approximated by rational functions [18], [20], and when mapped onto Löwner matrices, the degree of the rational function matches the rank of the Löwner matrix [17]. Modeling VA instead of AA is a reasonable strategy in these difficult cases (weak AA and/or short R-R intervals), common in persistent AF episodes, as the AA signal becomes very difficult to model [16]. The VA estimated by the Löwner-BTD is then subtracted from the ECG, resulting in a signal that contains mainly AA. Another advantage of the proposed method is that no technique for atrial source selection is needed, since the resulting signal is already the desired signal.
The main contributions of the paper can be summarized as follows: r We put forward a new tensor-based method for noninvasive AA estimation in the challenging scenarios of low amplitude f-waves or short R-R intervals, typical of persistent AF episodes.
r The new technique is based on a BTD variant based on Löwner matrices modeling the VA interference present in the recording. To compute the BTD with Löwner constraints, we develop a new method that we denote Löwner-constrained AGL (LCAGL), assuring the Löwner structure of its matrix factors.
r To better gauge its capabilities in the context of this biomedical application, the proposed method is compared with state-of-the-art tensor-based algorithms such as the Hankel-constrained AGL and the NLS-based Gauss-Newton algorithm to compute the Löwner-BTD.
r Experiments with realistic synthetic AF signals are performed to validate the proposed Löwner-based BTD algorithm in simulated scenarios where the ground truth is available.
r The clinical relevance of the proposed method is evaluated on a database composed of 20 different patients suffering from persistent AF. This work is a significant extension of the preliminary results reported in the conference paper [16]. Indeed, while the NLS method is used to compute the Löwner-BTD and tested in a reduced database in [16], we provide herein an improved method for Löwner-BTD computation and statistically more significant results by doubling the size of the AF patient database. This paper formulates a variant of CAGL that deals with Löwner constraints, imposing a low-rank Löwner structure. We show how this constraint can be enforced by an alternating projection approach, for which we derive an explicit expression of the projection onto the subspace of Löwner matrices. The resulting LCAGL is compared to the Löwner-BTD computed with the NLS method (used in the experiments reported in [16]), as well as the recently proposed CAGL dealing with Hankel constraints. Experimental results with synthetic and real AF ECG segments show the superiority of the proposed technique in challenging cases of this complex arrhythmia. Successful AA estimation is achieved from observation records as short as a single heartbeat, making the new method potentially very attractive in real-time AF analysis.
The rest of this paper is organized as follows. Section II presents the notation and some basic definitions of tensor algebra. Section III recalls BTD as a tensor approach to solve BSS problems and how to build a tensor from an ECG data matrix. Section IV describes AGL for computing BTD, while Section V formulates its variant with Löwner constraints and derives the linear projection onto the Löwner subspace. A numerical evaluation of the proposed approach with synthetic signals is shown in Section VI. Section VII presents the database setup and results of experiments with real AF data. Section VIII discusses the advantages and limitations of the present work and, finally, Section IX concludes the paper and provides some prospects of future research.

II. NOTATION AND ALGEBRA PREREQUISITES
To ease the comprehension of the paper, this section summarizes the notation used in the present work and introduces some definitions that are not common in research fields outside tensor algebra.
Scalars, vectors, matrices and tensors are represented by lower-case (a, b,...), boldface lower-case (a, b, ...), boldface capital (A, B,...) and calligraphic (A, B,...) letters, respectively. The matrix transpose operator is represented by (·) T , ||·|| F is the Frobenius norm and • is the outer product. Diag(a) builds a diagonal matrix by placing the elements of a along the diagonal and I N represents an identity matrix of size N × N . A third-order tensor A ∈ C I 1 ×I 2 ×I 3 , with scalar entries a i 1 ,i 2 ,i 3 , has its frontal slices represented by A ..i 3 ∈ C I 1 ×I 2 . A matrix A ∈ C I 1 ×I 2 , with scalar entries a i 1 ,i 2 , has its i th 1 row and i th 2 column represented by a i 1 . and a .i 2 , respectively. Symbol ||·|| denotes the l 2 -norm and ||·|| 2,1 denotes the mixed l 2,1 -norm, defined for an arbitrary matrix A with I 2 columns as: Given a vector y of length N , one can build a Hankel matrix with i = 1, . . ., I, j = 1, . . ., J.

III. TENSOR-BASED BSS OF ECG DATA
This section discusses the tensor-based BSS formulation and its connection with AF ECGs. First, Section III-A describes the matrix-based BSS formulation with its constraints to guarantee uniqueness, motivating the formulation of a tensor-based BSS problem with uniqueness under mild conditions. The Hankel-BTD for modeling the AA is recalled in Section III-B, whereas its Löwner-constrained version to model the VA is introduced in Section III-C.

A. BSS Formulation
An ECG recording composed of K leads and N time samples can be modeled as a matrix factorization: where M ∈ R K×R is the mixing matrix, modeling the propagation of the cardiac electrical sources from the heart to the body surface, S ∈ R R×N is the source matrix that contains the atrial and ventricular sources and R is the number of sources. Since the goal is to estimate M and S from matrix Y only, one can see that source estimation in ECG recordings is a BSS problem. As previously stated, matrix-based methods require some constraints in order to guarantee the uniqueness of the factorization in (4). For example, principal component analysis (PCA) [19] requires the columns of the mixing matrix M to be orthogonal, whereas independent component analysis (ICA) [7] methods require statistical independence between the sources and non-Gaussianity. Although mathematically coherent, such constraints may lack physiological grounds, hindering the interpretation of the results. By contrast, tensor-based techniques do not require such constraints and can ensure uniqueness under relaxed conditions.
As it will be detailed later, the BTD with constrained matrix factors can well suit the characteristics of ECG signals, which makes this tensor factorization technique a powerful tool for AA extraction from AF ECGs. The BTD of an arbitrary third-order tensor T ∈ R I×J×K is written as: where R is the number of blocks, c r ∈ R K is nonzero and E r ∈ R I×J is a structured matrix that has rank L r , admitting the factorization E r = A r B T r , where A r ∈ R I×L r and B r ∈ R J×L r have rank L r . We may then rewrite (5) as: One can see that the BTD is a decomposition of T in multilinear rank-(L r , L r ,1) terms [6], represented by a sum of the outer product of its matrix and vector factors, as illustrated in Fig. 2. Note that the number of sources R equals the number of blocks in the BTD model (5)- (6). Several results can be invoked to guarantee the essential uniqueness of this tensor decomposition. For example, in [6, Theorem 2.2], it is shown that the BTD is essentially unique if the following conditions are satisfied: 1) The matrix factors contain proportional columns. For the processed ECG segments, we typically have that Also, matrix C corresponds to the mixing matrix M, whose columns represent the contribution of each source to the ECG leads, which are unlikely to be proportional. Essentially unique means that the structure of the decomposition does not change if the r th 1 and r th 2 terms are permuted, as long as they have equal ranks, and if E r is rescaled, as long as c r is counterscaled. Milder uniqueness conditions can be found in [6].
If the matrices E r are linearly structured in the sense that they belong to a given subspace U S , then the frontal slices of T , given by T ..k , also belong to U S . For the BSS application of this decomposition, two commonly used structures are the Hankel structure, that focuses on separating complex exponential functions, and the Löwner structure, focusing on the separation of rational functions. In the following, we will describe the principles behind these choices.

B. AA Modeling via Exponential Functions
The Hankel-BTD is especially adapted to source signals composed of a sum of complex exponentials: where z f ∈ C is a possibly complex pole, λ f the associated scale factor (that can also be complex-valued), 1 ≤ f ≤ F , and t denotes continuous time. This model is indeed a good representation of the quasi-periodic behaviour of the AA signal in AF episodes [26].
It is shown in [6] that a BTD built from Hankel matrices solves the BSS problem (4) when the underlying sources follow the exponential model (7). This results relies on the fact that a Hankel matrix associated with an F -pole exponential model has a rank bounded by F , which follows from the Vandermonde decomposition of the Hankel matrix. To exploit this idea to perform BSS, a tensor is first obtained from the data matrix Y by mapping each of its K rows into a Hankel matrix H (k) ∈ R I×J , 1 ≤ k ≤ K, according to (2). Next, the tensor is built by stacking each Hankel matrix along the third dimension, as frontal slices, of a third-order tensor Y H ∈ R I×J×K , that is: S is a Hankel matrix built from the r th source of S. Because the described mapping is linear and each row of Y is a linear combination of the sources with weights determined by one row of M, it follows that the third-order tensor Y H admits a BTD model and can be written as: The Hankel-BTD suits the characteristics of AA in AF episodes since atrial signals can be approximated by a sum of complex exponentials like (7) and mapped into Hankel matrices with rank equal to the number of poles [8]. However, in persistent AF, ECG recordings present short R-R intervals, significantly masking the AA in the recording. Also, f-waves during AF are characterized for being disorganized with very low amplitude, sometimes lower than the noise. Such issues make it difficult to fit the AA signals to this exponential model, compromising the AA extraction via Hankel-BTD.

C. VA Modeling via Rational Functions
To circumvent the aforementioned difficulty in the persistent AF scenarios considered in this work, we turn to another alternative: VA modeling. Indeed, the QRS morphology is not significantly changed during AF, as compared to normal sinus rhythm (NSR), and hence it is reasonable to focus on VA estimation. When subtracting the VA estimate from the ECG, the resulting signal will ideally contain the AA signal.
Previous works notice that QRS complexes can be approximated by rational functions of low degree [17], [18], [20]: where a(t) is a polynomial of degree A, F is the number of complex poles, denoted p f , D f is their multiplicity, and c f, Interestingly, a second kind of BTD suits the characteristics of signal model (10). Instead of Hankel matrices, this alternative BTD is based on Löwner matrices built as in (3). The key result underying this decomposition is that when a signal following model (10) is mapped into a Löwner matrix, the degree of the rational function matches the rank of the matrix [17]. Hence, the idea is to map each row of the data matrix Y into a Löwner matrix L (k) ∈ R I×J as in (6), 1 ≤ k ≤ K. Next, the tensor is built by stacking each Löwner matrix along the third dimension, as frontal slices, of a third-order tensor Y L ∈ R I×J×K . Similarly as with Hankel matrices, one can show that: where L (r) S ∈ R I×J is a Löwner matrix built from the r th row of S. One can see that the procedure for constructing Y ..k in (11) is a linear mapping and, for each r, the outer product between matrix L (r) S and the r th column of M, i.e., m .r , is performed to build a third-order tensor containing the contribution of the r th source to the ECG tensor. Putting together the contribution of all sources, the third-order tensor Y L admits a BTD model and can be written as: In short, the Löwner-BTD tensor model (12) is directly related to the rational function model (10), and is therefore well-suited to VA estimation in AF ECGs. To do so, one first needs to estimate the BTD model factors from third-order tensor Y L , an important step that is addressed next.

IV. ALTERNATING GROUP LASSO ALGORITHM FOR BTD
In this section, the classical approach for the computation of the BTD model factors and its limitations are recalled. Then, we review the AGL algorithm, which was conceived to overcome such limitations.
In general, an approximate BTD is computed by minimizing the Euclidean distance between the observed data tensor Y ∈ C I×J×K and a model of fixed structure with respect to the model components: In the special case of interest shown in eqn. (12), L r must belong to the subspace of Löwner matrices with dimensions ( N 2 × N 2 ), denoted S L . The mode-3 slices Y ..k of the observed tensor are Löwner by construction. However, a solution (Â,B,Ĉ) of (13) may not satisfyÂ rB T r ∈ S L , due to noise and modeling imperfections. Note that even if the sum R r=1 (Â rB T r )ĉ k,r is Löwner, there is no guarantee that so are matricesÂ rB T r , r = 1, 2, . . ., R. Also, algorithms based on (13) are strongly dependent of the initialization of their matrix factors and do not estimate the model parameters, i.e., the number of blocks and their ranks.
To overcome such limitations, instead of using a fixed BTD structure as (13), an algorithm called AGL and its constrained version for Hankel matrices called CAGL, are introduced in [14]. This method includes penalization terms promoting low-rank blocks and controlling the number of blocks as: where f (·, ·, ·) is the same as in (13), γ > 0 is a regularization parameter and g is a regularization function of the form: Due to the geometric properties of the mixed 2,1 -norm, solutions where A, B and C have null columns (for sufficiently high γ values) will be induced, allowing one to select the relevant low-rank blocks. This method is called group lasso and is a generalization of the the lasso estimator principle [22].
Since problem (14) is nonconvex (and nonsmooth), but convex by blocks, a block coordinate descent (BCD) approach is employed in [14]. BCD consists in partitioning the set of optimization variables and sequentially solving convex subproblems in each subset of variables, fixing the others.
Analogously, strictly convex subproblems can be derived for B and C. AGL solves these subproblems alternatively, through the updates ofÂ (ν) ,B (ν) andĈ (ν) , in this order, at each iteration ν. When all subsets are updated, one iteration of the algorithm is completed. Iterations are repeated until convergence.

V. LÖWNER-CONSTRAINED ALTERNATING GROUP LASSO
In order to ensure the Löwner structure of the matrix factors at the end of iterations, this section explains how Löwner constraints can be imposed in CAGL, and then details the linear projection that can be implemented in SLRA algorithms.

A. Imposing Löwner Constraints in L r
We propose a formulation to deal with Löwner constraints in CAGL, called here Löwner-constrained AGL (LCAGL). We highlight that LCAGL follows the same principle of CAGL, which is to ensure a specified structure over the block matrices. However, instead of using a projection on the Hankel subspace as described in [14], LCAGL uses an optimal projection on the Löwner subspace that will be explicitly given later.
Löwner constraints can be enforced by solving a SLRA problem for each block matrix at each iteration of LCAGL. For this purpose, Cadzow's Algorithm (CA) [21] is used at the end of each iteration. CA performs alternating projections, where the one that leads to the set of matrices with rank up to L r is performed by the truncated singular value decomposition (SVD), whereas the one that leads to the Löwner subspace is formulated below. It is important to highlight that, in practice, the rank used in the SVD truncation step is computed by counting the current number of simultaneously nonzero columns of A r and B r . The algorithm is summarized in Table I.

B. Orthogonal Projection Onto the Löwner Subspace
In this work, we consider for simplicity signals with an even number of samples, i.e., I = J = N/2, and the interleaved partitioning method. In this way, we can rewrite the matrix L (k) ∈ R N 2 × N 2 of (3) as: Considering that the signals are regularly sampled with a sampling period Δ = t n − t n−1 , for n = 2, . . ., N, (17) becomes: where each element is given by: Note that this formulation can be easily adapted to signals that are irregularly sampled. For this other scenario, the denominator (19) is replaced by Δ ij = t 2i−1 − t 2j . Since in practice the signals are regularly sampled, the formulation of (19) will be kept in the sequel.
For a given matrix E (k) ∈ R I×J , its projection onto the Löwner subspace can be obtained by minimizing the LS cost function given by: where a i =ŷ k,t 2i−1 and b j =ŷ k,t 2j , for i, j = 1, . . ., N/2, witĥ y k. being the estimate of y k. . Taking the derivative of F with respect to an arbitrary a p we have: Then, In this way, the solution for a p , for p = 1, . . ., N/2, can be obtained as: Similarly, for b p we have: The necessary conditions for optimality can be expressed as: for p = 1, . . ., N/2, where: and Hence, a system of linear equations can be solved in order to estimateâ = [â 1 , . . .,â N/2 ] andb = [b 1 , . . .,b N/2 ]: where matrices U and V, and vectors α, β, z and w are constructed from u p,j , v i,p , α p , β p , z p and w p in (29), (30), (31), (32) and (33), respectively. Note that this linear system has an infinite number of solutions because the block matrix of (30) has rank N − 1, as a consequence of the vertical shift invariance of the Löwner structure. Indeed, the Löwner construction is invariant to a vertical shift of the signal because the Löwner matrix built from signal y k. is equal to the one built from y k. = y k. + c, for any constant c ∈ R. To fix this ambiguity, without loss of generality, we can assume that the first sample is null, i.e., a 1 = 0, so that the linear system has now N − 1 linearly independent equations in N − 1 variables, thus providing a unique solution. A matrix with Löwner structure can then be reconstructed from estimated vectorsâ andb.

VI. NUMERICAL EVALUATION ON SYNTHETIC DATA
This section presents some numerical results on synthetic data, which are convenient for assessing the performance of LCAGL since in these scenarios the ground-truth is known. This assessment is carried out using two different synthetic signals that simulate persistent AF recordings. Both signals contain 12 leads, generated by a random mixing matrix and 3 sources: AA, VA, and noise.
To simulate the VA signals, two synthetic QRS complexes modeled by rational functions were generated according to the model proposed in [25]: where Θ is the coefficient parameter, responsible for the symmetric/asymmetric behavior of the model, R(·) is the real part of its complex argument, and r m a [·] is a basic normalized rational function given by: whereā is the reciprocal of the pole a and m is a multiplicity parameter. Two synthetic QRS complex models with parameters a = 0.8, m = 2, and Θ = {π/2, 0} are generated in this manner.
To simulate the AA signal during AF, the model proposed in [26] that mimics the f waves is used. This model with P harmonics is given by: with modulated amplitude and phase respectively given by: and where x is the sawtooth amplitude, Δx is the modulation peak amplitude, f x is the amplitude modulation frequency, f s is the sampling frequency, f 0 is the frequency value upon which θ[n] depends linearly, Δf is the maximum frequency deviation and f f is the modulation frequency. Two synthetic f-wave models were generated with the parameters presented in Table II.
To simulate the interference typically present in an ECG, additive white Gaussian noise (AWGN) with zero-mean and variance σ 2 , respectively, is introduced. Finally, the mixing matrix is also generated according to a Gaussian distribution, with scaling factors chosen to obtain an average power ratio between the components consistent with clinical ECGs.
One generated ECG signal mimics the challenging scenario where the AA is weak, while the other synthetic signal mimics the scenario where the R-R interval is short. Both synthetic AF ECG signals are illustrated in Fig. 3.
The normalized mean square error (NMSE) between the estimated and original VA source signals is computed to evaluate estimation quality. As LCAGL is relatively robust to initialization, no Monte Carlo runs were performed and the value of γ between 8 × 10 −4 and 10 −2 that provided the best solution was chosen. Fig. 4 shows the VA estimates (blue dashed line) and the respectives ground-truths (gray solid line) for the two synthetic signals along with the computed NMSE described previously and the rank of the block. The low NMSE shows that a successful source extraction is achieved while the Löwner structure of the block is guaranteed. Also, one can see that the rank of the block, rank B , increases with the number of QRS complexes in the segment, as expected, since the signal has more transient components, needing more poles to be well modeled.

VII. EXPERIMENTAL RESULTS IN REAL AF ECGS
In order to validate the proposed tensor approach, experimental results using real AF ECG recordings, comparing LCAGL with the Löwner-BTD computed using the NLS method, and with CAGL imposing Hankel constraints, are reported and discussed in this section. First, the indices used to measure AA quality and the observed database are summarized. Then, the experimental results from the two challenging scenarios of this arrhythmia, which were previously described in Section I, are shown and discussed.

A. Signal Quality Measurement in AF Episodes
Measuring estimation quality (or the AA content) of real signals is a difficult task. Since in practice there is no ground truth for comparison, one needs to take advantage of some features which are typically present in AA during AF episodes and use them as surrogate performance criteria. For example, in the frequency domain, the AA during AF has a peak between 3 and 9 Hz. The position of this peak is called dominant frequency (DF). In this section, two parameters used to measure AA estimation quality are presented. The first one, used to measure AA extraction quality, is the spectral concentration (SC), defined as the relative amount of energy around the DF. The SC is computed as in [23]: where f p is the value of the DF, F s is the sampling frequency, f i is the discrete frequency and P S is the power spectrum of the source signal computed using Welch's method, as in [23].
The second parameter, that also provides a quality measurement of AA extraction, is the power contribution to lead V1, denoted P (r), given by [9]: in mV 2 , where m V 1,r is the contribution of the r th source to lead V1 and s r. is the r th source in time domain. The P (r) of an AA source is expected to be relatively strong (> 10 −4 mV 2 ), since lead V1 is the one that typically best reflects AA in AF ECGs, as this lead strongly correlates with the AA from the right atrium and moderately correlates with that from the left atrium [24].

B. AF Database and Experimental Setup
Real ECG data from 20 patients suffering from persistent AF are used in this study. All the recordings belong to a database provided by the Cardiology Department of Princess Grace Hospital Center, Monaco. The recordings are acquired at a 977 Hz sampling rate and are preprocessed by a zero-phase forward-backward type-II Chebyshev bandpass filter with cutoff frequencies of 0.5 and 40 Hz, in order to suppress high-frequency noise and baseline wandering.
Out of the 20 different segments of ECG recordings, 10 segments have a disorganized and/or weak AA, and the other 10 have very short R-R intervals. Because in a standard 12-lead ECG only 8 leads contain linearly independent signals, the other ones being combinations thereof, we process only these independent measurements (I, II, V1-V6) in order to reduce computing cost while still exploiting all available spatial diversity. The two types of segments used in these experiments are illustrated on lead V1 in Fig. 5 (solid lines): short R-R intervals (bottom line) and disorganized/weak AA (top line).
The 10 segments with disorganized and/or weak AA, from patients P1 to P10 are composed of one heartbeat, i.e., the QRS complex followed by the T wave and the visible f waves, and have between 0.8 and 1.7 seconds. The other 10 segments with short R-R intervals from patients P11 to P20 have 1.5 seconds of duration with at least 2 QRS complexes. All the 20 segments are downsampled by a factor of 10, since this significantly reduces the computational complexity without any noticeable loss of quality, as the signals' energy is concentrated in lower frequencies and no spectral aliasing is generated [14].
The NLS method used for comparison in the experiments is available in Tensorlab MATLAB Toolbox and is set up as in [16].

C. Segments With Disorganized and/or Weak AA
In Fig. 6, we can see the observed recording as well as the estimated VA and AA of the segment of Patient P2 processed by LCAGL. The segments are shown now on lead V1 for a better clarity of the estimated AA, in the time domain on the left and in the frequency domain on the right, with a DF equal to 6.67 Hz, typical of AA. It can be seen that the SC provided by LCAGL is 58.73% and P (r) value equal to 2.02 × 10 −4 mV 2 , while the ones provided using the Löwner-BTD computed by the NLS method, choosing the best match (R, L r ), are 5.36% and 1.59 × 10 −4 mV 2 , respectively, for the same segment. It should be noted that the Hankel-BTD computed by CAGL was not able to extract the AA from this recording, nor from any other used recording characterized by a disorganized and/or weak AA. Indeed, Hankel-BTD attempts to model the AA present in the ECG, which is a difficult task under these conditions. Due to the weak AA amplitude, its estimate is basically a straight line as shown in Fig. 5 (Top).
Expanding the previous result to the population of 10 patients with segments presenting this particular characteristic, Fig. 7   (Top) shows the SC values of the AA source estimated by the Löwner-BTD computed by the NLS method and LCAGL, for each patient. Regarding this AA quality index, LCAGL outperforms the NLS method in all but patient P4, where the SC of NLS is 1% higher than the one provided by LCAGL. In addition, the NLS method was not able to successfully extract the AA from patient P7, and CAGL could not successfully extract the AA from any of the 10 observed patients, as stated before. Also, the average SC over the 10 patients is 45.48% for LCAGL and 27.12% for the NLS method. Fig. 8 shows how the P (r) values vary through the population of patients for both Löwner-based methods. No significant difference between the methods can be observed. However, a significant difference is observed between the two patient populations. The P (r) of the patients with a disorganized and/or weak AA signal is lower than the ones characterized by short R-R intervals and varies in a shorter range. These factors are indicative of fine AF.  Fig. 7 (bottom) shows the SC values of the AA source estimated by the Löwner-BTDs computed by the NLS method and LCAGL, for each patient whose segments are characterized by short R-R intervals. In 6 out of the 10 patients of this population, LCAGL provided an AA estimate with significantly higher SC than NLS, while for the other 4 patients the SCs of both methods are practically the same. More precisely, for those 4 patients LCAGL still produces a slightly higher SC value (around 1%). In addition, the mean SC over these 10 patients is 45.72% for LCAGL and 39.25% for the NLS method. Fig. 9 shows the VA and AA estimates on lead V1 of the processed segment of patient 20. One can see that the several QRS complexes are well estimated, resulting in a clear AA signal when subtracted from the original ECG. The SC provided by LCAGL is 52.48%, practically equal to the 52.45% provided using the NLS method, manually choosing the best match (R, L r ) for the same segment. LCAGL provides a slightly stronger P (r), than the NLS method, which yields 3.11 × 10 −3 mV 2 for this segment. Fig. 10 shows the observed segment of Patient P18 on lead V1 and compares AA source estimations by the tensor-based methods. The AA estimate by CAGL (blue line) is deformed, probably due to the presence of 2 QRS complexes and the shape of the AA that may not be well approximated by an all-pole model, hampering the AA extraction. This loss of physiological shape also occurs in Fig. 5 (bottom). Although the AA source estimated by NLS is better than the one produced by imposing the Hankel structure, it contains some ventricular residuals, due to an inaccurate VA source estimation. Finally, it can be seen that for this segment, LCAGL is the method that provides the best VA cancellation from the original recording, yielding the best AA estimate (green line).

D. Segments With Short R-R Intervals
We note that, for all but patients P8 and P17, both Löwnerbased methods (NLS and LCAGL) provide the same DF for the AA source estimate. Before concluding this experimental assessment, it is important to highlight that classical techiques that estimate the VA and subtract it from the ECG to retrieve the AA, such as average beat subtraction and adaptive singular value cancellation [27] will certainly fail in theses short recordings, as they need sufficiently long ECGs to perform an accurate VA estimation.

VIII. DISCUSSION
Extracting the AA from AF ECGs is a particularly challenging task when the AA presents a very low amplitude and short R-R intervals appear often along the recording, as typically occurs in persistent forms of AF. In these cases, methods that focus on f-waves estimation do not provide a clean AA extraction or even fail. The proposed LCAGL method overcomes such limitations by taking into account that QRS complexes can be well modeled by rational functions. Estimating the VA can be simpler in such scenarios as the QRS morphology is not significantly changed during AF episodes as compared to NSR.
When compared to other tensor-based methods, LCAGL is more robust to the initialization of model parameters, as observed in the experiments on real and synthetic data. Also, the method guarantees the Löwner structure of BTD matrix factors. As compared to classical techniques for AA extraction based on QRS cancellation, LCAGL can operate successfully in very short ECG records, about one heartbeat long, whereas competing techniques require recordings typically composed of tens of seconds to estimate the QRS template well. This interesting advantage is a direct consequence of the deterministic signal model assumed by the BTD tensor approach, and renders LCAGL very attractive for real-time AF analysis. Another advantage lies in the fact that, after source separation, no additional processing is needed to automatically identify the AA source of interest, since the AA contribution results from the subtraction of the estimated VA from the original ECG. In consequence, potential errors resulting from AA source misidentification are avoided.
Limitations of the present work include the relative short sample size (20 patients) of the persistent AF population used in the experimental assessment, from which statistically solid conclusions are difficult to draw. The instability of LCAGL in estimating the T-wave should also be mentioned. Indeed, Fig. 6 shows an example where the T-wave is also captured with the QRS complex, whereas in Fig. 9 the T-waves are missing in the VA signal estimated by the same method. A possible explanation is that it is quite hard to distinguish the T-wave from the AA signal in segments with short R-R intervals, as one can see in Fig. 9. It should be remarked, however, that T-waves pose difficulties to most AA extraction methods in the literature.

IX. CONCLUSION
Persistent AF is an advanced stage of this arrhythmia where the ECG typically presents a fast heart rate, resulting in short R-R intervals and an AA signal that is difficult to identify, making the AA extraction a challenging task. The present work has put forward and investigated a tensor-based approach to model the VA and separate it from the AA in persistent AF ECGs, showing a number of benefits in these difficult scenarios. This approach consists in explicitly modeling the VA with rational functions and mapping them into Löwner matrices in order to built a BTD tensor model. The factors of this third-order tensor model are computed using the LCAGL algorithm introduced in this work. The Löwner structure of the matrices that compose the tensor is guaranteed by a linear projection onto the Löwner subspace derived here.
After an evaluation on synthetic signals, the Löwner-BTD computed by LCAGL was applied to a population of 20 patients suffering from persistent AF, whose segments are characterized by short R-R intervals and a disorganized and/or weak AA signal contribution. Experimental results with real AF ECGs in very short observation windows (around 1-1.5 seconds) showed LCAGL's better AA extraction performance than two other tensor-based methods in this challenging clinical scenario, while offering the possibility of beat-to-beat on-line processing.
More robust tensor-based techniques should be studied in future works, in order to derive methods for the computation of an approximate BTD model that requires less parameter tuning and is able to better capture the T-wave when hard to discern from the f-waves. Also, experiments in a larger population of AF patients are envisaged in order to confirm the clinical relevance of these results.