Block Term Decomposition of ECG Recordings for Atrial Fibrillation Analysis: Temporal and Inter-Patient Variability

—Responsible for 25% of strokes and 1/3 of hospitalizations due to cardiac related disturbances, atrial ﬁbrilla-tion (AF) is the most common sustained cardiac arrhythmia in clinical practice, considered as the last great frontier of cardiac electrophysiology. Its mechanisms are not completely understood, and a precise analysis of the atrial activity (AA) signal in electrocardiogram (ECG) recordings is necessary to better understand this challenging cardiac condition. Recently, the block term decomposition (BTD) has been proposed as a powerful tool to noninvasively extract AA in AF ECG signals. However, this tensor factorization technique was performed only in short ECG recordings and illustrated in single patients. To assess its performance and variability through different subjects, BTD is applied to a population of 10 AF patients in this paper. Also, its time variability is evaluated by means of long segments of AF ECG with varying observation window size. Experimental results show the consistency of BTD as an AA extraction tool, outperforming two well-known matrix-based methods in most of the observed cases for long and short AF ECG recordings.


I. INTRODUCTION
C ARDIAC arrhythmias are heart diseases characterized by abnormal electrocardiography waveforms.Reducing life expectancy and increasing the healthcare bill, atrial fibrillation (AF) is the most common sustained cardiac arrhythmia encountered in clinical practice, especially affecting the eldery.Persistent AF represents a particularly complex case of this arrhythmia, where extensive atrial remodelling has taken place due to sustained AF, significantly affecting atrial activity (AA) and AF perpetuation itself.The importance given to this challenging cardiac condition has increased in the last years and is expected to increase further, becoming a major public health and economical concern, as its mechanisms are complex and not completely understood.About 160 000 new patients are diagnosed with AF every year only in USA and similar numbers are reported in European countries.By 2050, it is estimated that AF will increase from 2.3 to about 12.1-15.9 million patients only in USA, becoming then a new epidemic [1].Moreover, a patient with AF spends, anually, Pedro Marinho R. de Oliveira and Vicente Zarzoso are with the University Côte d'Azur, CNRS, I3S Laboratory, France (e-mail: {marinho,zarzoso}@i3s.unice.fr).
Digital Object Identifier: approximately $8 700 more in healthcare than a patient without AF.The treatment of this desease is estimated to add $26 billion per year to the USA healthcare costs [2].In a normal electrocardiogram (ECG), the P wave corresponds to atrial depolarization.During AF, the P wave is replaced by fibrillatory waves, called f waves, which are present throughout the whole ECG recording.However, they are masked by the QRST complex that corresponds to the ventricular depolarization and repolarization, i.e., the ventricular activity (VA), in each heartbeat.Signal processing techniques are useful and important tools in the analysis of cardiac signals, especially AF signals.In particular, they are necessary to noninvasively extract the AA from the standard 12-lead ECG for an accurate analysis of the f waves, in order to better understand the complex mechanisms behind this disease.The extraction of AA from multi-lead ECG recordings accepts a blind source separation (BSS) formulation [5].Matrix-based techniques to solve BSS problems, such as principal component analysis (PCA) [3] and independent component analysis (ICA) [4] have proven useful in noninvasive AA extraction [5]- [7].However, matrix decompositions are known to have some limitations, since constraints need to be imposed to guarantee uniqueness, e.g., orthogonality or statistical independence between the sources.Such constraints are mathematically convenient, but they may lack physiological grounds, making difficult the results interpretation.
In order to overcome the limitations imposed by matrix decompositions, a tensor factorization technique has been recently proposed to noninvasively extract the AA signal from AF ECG recordings [8], [9].If compared to matrix-based techniques, tensor decompositions present some remarkable features like, for example, their essential uniqueness with minimal or no constraints.Another example is the fact that the rank of the tensor can exceed its largest dimension, whereas in matrices the rank is limited by its lowest dimension.The block term decomposition (BTD) proposed as a technique to solve BSS problems in [10] suits the characteristics of AA in AF episodes, since atrial signals can be approximated by all-pole (or exponential) models and mapped onto Hankel matrices with rank equal to the number of poles [12].The Hankel matrices containing the ECG data are stacked along the third dimension of a third-order tensor, and then processed by BTD.Previous experiments using BTD in synthetic and real AF ECG data showed that this tensor decomposition can outperform the matrix-based techniques for noninvasive AA extraction in short segments of AF ECG recordings [8], [9], [11]- [13].
The performance of BTD strongly depends on its computation method.The most used for this application is the non-linear least squares (NLS) method, whose results, in turn, depend considerably on its initialization and model parameters, i.e., the multilinear rank (the rank of the matrix factors) and the number of block terms of the tensor.In [12] an estimation method of the multilinear rank of BTD in an AF ECG tensor is proposed.This estimation is based on the application of three popular methods for the estimation of auto-regressive (AR) model order in the TQ segment of a heartbeat, where only AA is present.In addition, a sucessful AA extraction depends on the precise selection of the AA source among the estimated sources.In [13] the BTD shows that the classical automated method of atrial source selection [5], [6] is not optimal, and then two automated methods are proposed that, though still suboptimal, provide better accuracy.
BTD has proven useful in extracting the AA signal from the AF ECG, being able to outperform well-known matrixbased techniques.However, the performance of this tensor approach was only assessed in fixed short segments of a single AF patient.Its performance in a population of patients remains an open question.Aiming answer this open question and to provide results with more clinical relevance, this paper presents the performance of BTD in long AF ECG segments, as well as the analysis of the observation window size on its performance [14].Furthermore, the present work assesses BTD for the first time in a population of patients with persistent AF and compares it to two popular matrix-based methods for AA extraction: RobustICA-f [15] and PCA.The correlation of two AA quality indices is another outcome of this work.Experimental results using Monte Carlo simulations evaluate the AA extraction performance of BTD on long and short segments of AF ECG recordings.
The rest of this paper is organized as follows.Section II introduces the notation used in the present work, while Section III recalls the BTD as a tensor approach to solve BSS problems.Section IV discusses quantitative measures of AA content.Section V presents the database and the setup used in the present experiments, whose results are reported in Section VI.Finally, Section VII gives an end to this paper with conclusions of the work and prospects of future research.
The matrix transpose is represented by (•)

III. BLOCK TERM DECOMPOSITION
The BTD of an arbitrary third-order tensor T ∈ R I1×I2×I3 is written as where c r ∈ R I3 is a nonzero vector and E r ∈ R I1×I2 has rank L r and admits a factorization E r = A r B T r , where A r ∈ R I1×Lr and B r ∈ R I2×Lr have full column rank L r .Equation (1) may then be rewritten as The visual representation of BTD as a sum of the outer product of its matrix and vector factors is shown in Figure 1.One can see that the BTD is a decomposition of T in multilinear rank-(L r ,L r ,1) terms.Several conditions guarantee the essential uniqueness of this decomposition.For example, in [10, 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.The first condition requires that R r=1 L r ≤ I 1 , I 2 .Milder uniqueness conditions are also presented in [10].Now, ECG recordings from K leads composed of N time samples can be stored in a matrix: 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, ventricular and other sources (noise, respiration, muscular activity, etc.), and R is the number of sources [7].Since the goal is to estimate M and S from the matrix Y (the only observed data), it is clear that AA extraction in an AF ECG recording is a BSS problem.
As previously described in [10], the BTD is proposed as a solution of a BSS problem like (3).The idea to obtain a tensor from Y is to map each of its k th row onto a Hankel matrix where i = 1, ..., I, j = 1, ..., J, and k = 1, ..., K. Next, the tensor is built by stacking each Hankel matrix along the third dimension (as frontal slices) of a third-order tensor Y ∈ R I×J×K , that is In scalar form, the third-order tensor Y can be written as The k th matrix slice of this tensor can be represented as where ∈ R I×J is a Hankel matrix built from the r th row of S. One can see that for each r, the outer product between matrix H (r) S and the r th column of M, i.e., m .r, is being performed to obtain the contribution of each source to the observed tensor.This way, the third-order tensor Y can be written as Comparing Equations ( 1) and ( 8), it can be concluded that the ECG data tensor follows a BTD tensor model with the following correspondence Due to the quasi-periodic nature of AF, atrial sources can be represented by the all-pole model: where n = 1, ..., N represents the discrete-time index, r = 1, .., R, L r is the number of exponential terms, z ,r is the th pole of the r th source, and λ ,r is the scaling coefficient.This way, their associated Hankel matrices accept the Vandermonde decomposition [16]: with and In the case of different poles, the Vandermonde matrix with L r ≤ I, J will have full-column rank L r , so if M does not have proportional columns, the BTD in (8) will be essentially unique.In the case of equal poles, milder conditions can assure the uniqueness of (8) [10].
As previously described, BTD can outperform the matrix approaches regarding AA extraction in fixed short AF ECG recordings [8]- [13].However, its performance in a population of AF patients remains an open question, since, to the best of the authors' knowledge, the BTD performance has only been assessed in fixed short AF ECG segments of single patients.

IV. ATRIAL ACTIVITY CONTENT MEASUREMENT
Signal processing techniques used to solve BSS problems separate the observed signal in several sources.In AF ECGs, typically at least one of these sources contains the AA.It is unknown if the AA is concentrated only in a single source.However, in the present work, as in previous works, it is considered that the AA is concentrated only in one source, and the source with the most significant AA content is called the atrial source.
Measuring the quality of the estimation or the AA content of real signals is a challenging task.Since there is no ground truth for comparison, one needs to take advantage of some features present in AA during AF.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).As in [13], it is defined as potential atrial source any source with DF in such an interval.In this section, three parameters used to measure AA extraction quality are presented.In this work, the atrial source is selected using the automated method that provided the highest accuracy proposed in [13], helped by visual inspection.The first parameter for AA content measurement is the spectral concentration (SC), that is, the relative amount of energy around the DF.The SC is computed as in [6]: 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 [6].The second parameter used to better analyze the potential atrial sources is the kurtosis, denoted K, of the signal in the frequency domain acquired by a 4096-point FFT.As in [15], the general expression of kurtosis valid for noncircular complex data is used, given by where S r is the FFT of the r th source.Kurtosis is a measure of peakedness and sparsity of a distribution and, when computed in the frequency domain, it naturally provides a quantitative measure of harmonicity of the signal.A high kurtosis in the frequency domain means that the power spectral density is sparse, which is suggestive of a quasi-harmonic signal like AA in AF.In addition, kurtosis is parameter free, whereas SC depends on the DF and the definition of a suitable interval for interpretation [13].The third parameter, used to discard sources with irrelevant content that can be mistaken as the atrial source, is the power contribution to lead V1, which is given by in mV 2 , where m is the contribution of the r th source to lead V1 (given by the corresponding element of the estimated mixing matrix) and s r. is the r th source in time domain, corresponding to the r th row of matrix S in Equation (3).The power contribution to lead V1 by an AA source is expected to be relatively strong (> 10 −4 mV 2 ), since this lead is the one that typically best reflects AA in AF ECGs.

V. DATABASE AND EXPERIMENTAL SETUP
Experiments on real AF ECG data evaluate the performance of BTD regarding AA extraction in long segments varying the observation window size and short segments in a 10-patient population.

A. Real AF ECG Data and Preprocessing
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.

B. Data for Time Variability
To assess time variability in long segments (≥ 2.5 seconds), experiments are performed in 4 segments varying the observation window size of a standard 12-lead ECG recording from a single patient suffering from persistent AF.All 12 leads are considered.A 15-second segment in lead II from this patient is shown on Figure 2. The segments assessed for this patient have about 2.5, 5, 10 and 15 seconds and they all have the same starting sample.They are downsampled by a factor of 4, 8, 16 and 24, respectively, since the originally built third-order tensors pose some difficulties to Tensorlab MATLAB toolbox [17].Downsampling includes a low-pass filtering with cutoff frequency f c = f s /2, where f s is the new varying frequency.For the matrix-based techniques PCA and RobustICA-f, no downsampling is needed.

C. Data for Inter-Patient Variability
To assess the performance in a population of patients, and providing more relevant clinical results, experiments consider a randomly selected heartbeat (QRS-T complex + TQ segment) of ten real standard 12-lead AF ECG recordings from ten different patients suffering from persistent AF.A singlebeat segment in lead II from one of the patients in the 10patient population is shown on Figure 3, where the TQ interval can be seen just after the QRS-T complex.All the beats (one per patient) have between 1000 and 1400 samples (1.02 and 1.43 seconds).They are downsampled by a factor of 2, for the same reason previously explained.Fig. 3.A single heartbeat segment of an AF ECG recording, shown in lead II, from one patient of the observed population.A heartbeat consists in the QRS-T complex, followed by the TQ segment, where only AA is observed.

D. BTD Setup
BTD is implemented using the NLS method available in Tensorlab choosing R = 12 and L r = 95, for all r = 1, 2, ..., R.This choice is made based on the work [12], that shows that such values provide satisfactory results for the heartbeat with the largest TQ segment of the patient considered in the present work for the long segment performance assessment.The tolerance threshold for BTD convergence is set to 10 −9 and the maximum number of iterations is set to 1000.Ten Monte Carlo runs, with Gaussian random initialization for the matrix and vector factors at each run, are used to analyze the performance of BTD in AF patients regarding AA signal extraction.Monte Carlo runs are needed since the NLS method to compute BTD has a strong dependence on initialization.This tensorial technique is compared to the matrix-based methods PCA and RobustICA-f, which have already proven their efficiency solving BSS problems.

A. Temporal Variability Assessment
Table I shows the values of SC in % for PCA and RobustICA-f.The mean of the ten independent runs and the maximum value is shown for BTD.It can be seen in Table I that the mean of the SC for BTD is very close to the matrix-based techniques, outperforming PCA in all the observed segments and RobustICA-f in 3 out of the 4 segments with different lengths.Also, BTD is superior for all observed segments in 5 to 8 out of the 10 independent runs, showing that it can easily have superior performance if the right model parameters (R, L r ) and initialization are chosen.It can also be seen that the maximum SC presents a high value (≥ 75%) for all the observed segments, and in all cases, superior to that of matrix techniques, meaning that, with the right initialization, BTD outperforms the matrix-based techniques for the considered window lengths.
In the observed segments with different lengths, BTD finds more potential atrial sources than PCA and RobustICA-f.Discovering more than one potential atrial sources may be an interesting outcome, since it increases the possibility of finding some features that, although weakly contributing to the overall AA, may provide useful physiological and clinical information about the arrhythmia.This possibility, however, will not be explored in the present work.
Figures 4 and 5 illustrate the estimated atrial source by the three BSS techniques compared in this paper, in the time and frequency domain, respectively, for the best performance of BTD in the 15-second segment.These two figures show the satisfactory performance of BTD in estimating the AA in long segments of an AF ECG, as well as its superiority compared to the matrix-based methods, as quantified by the higher SC and kurtosis values.
The observation window size used in the experiments extends from 2.5 to 15 seconds in order to analyze its influence on BTD performance.Figures 6 and 7 show the variation of SC and kurtosis, respectively, over the 10 independent runs for all the window observation sizes analyzed in this work.This variation is illustrated by a box-and-whisker plot, where the red lines represent the median and the edges of the blue boxes the 25th and 75th percentiles.The black whiskers represent the extreme data values that are not considered as outliers, which are represented by red dots.
It can be seen that there exists a certain variation of SC and kurtosis over the runs, for each observation window size.This is expected, since BTD performance depends considerably on the initialization of its model factors.However, there is no clear trend in the computed parameters, which seems to indicate that the influence of the observation window size on the performance of BTD is not very significant or critical.The only drawback observed when processing long segments of ECG using BTD is the fact that the original segment must be downsampled by high factors, since Tensorlab has some difficulties in processing a tensor with large dimensions.Downsampling by a high factor could cause some loss of information in the signal, due to the frequency filter used to avoid spectral aliasing.However, previous experiments in short segments have shown that the impact of the downsampling factor on the SC is minimal.Actually, it was observed that the downsampled segments provides a slightly lower SC comparing to the original segment (without any downsampling), which makes the results reported in the present paper a "lower  bound" to the performance of BTD in the observed scenario.

B. Inter-Patient Variability Assessment
Table II shows the values of SC in % for PCA and RobustICA-f in the 10-patient population.The mean of the ten independent runs and the maximum value is shown for BTD.We can see in Table II that the mean of the SC for BTD is very close to the matrix-based techniques, outperforming PCA in almost all the observed patients (except patient 8) and RobustICA-f in 5 patients.Also, BTD is superior for all patients in 3 to 7 out of the 10 independent runs, showing again that it can easily have superior performance if the right parameters (R, L r ) and the right initialization are chosen.We can see also that the maximum SC is satisfactory (≥ 65%) for all but one patient (P10), and in all cases, superior to that of matrix techniques, meaning that, with the right initialization, BTD was superior in all cases.In the observed population of patients, BTD found more potential atrial sources than PCA and RobustICA-f.As previously explained, finding many potential atrial sources may be an interesting result, increasing the possibility of finding features that weakly contribute to the AA, while providing important physiological and clinical information about a complex arrhythmia like persistent AF.
As explained in Section IV, this work assumes, as in previous works, that the AA is concentrated on a single source, which is selected as the source presenting the most significant AA content, measured by the parameters previously presented.However, in the present experiments, in some runs of 5 from the 10 patients, AA appears in more than one source estimated by BTD.This can be seen in Figure 8(a), where two estimated sources by BTD presents significant AA (signals 2 and 3).This could mean that BTD is able to extract more information about the AA than the other methods compared here.To keep the focus on the subject of this work, a deep analysis of the cases where BTD provided more than one source with significant AA is not discussed in this paper.
Figure 8 shows the 3 most relevant potential atrial sources of one of the observed patients (Patient 6).For clarity, all the potential atrial sources with power contribution to lead V1 less than 10 −4 are not showed here, since they do not present significant features as they are very weak.In the time domain, shown in Figure 8(a), we can see that the two last sources have the atrial signature and present a high power contribution to lead V1, while looking at the frequency domain, in Figure 8(b), we can see that those sources present high kurtosis and SC, represented by indices K and SC in the figure legend, respectively.
Figure 9 shows how the SC of the atrial source is distributed over the 10 independent runs for the population of 10 observed AF patients.We can see that the median (red line) and the percentiles (blue box) present significant variations over BTD initialization and over patients.A fixed parameter match and Gaussian random initializations for each run are used to compute the BTD in 10 AF patients, providing satisfactory results.However, it can be observed that the chosen parameters work better in some patients than in others, recalling that BTD performance depends strongly on its parameters and initialization, and opening the challenge of finding the best parameter choice for a particular patient.It can be pointed out that in patients 9 and 10 the matrixbased methods could not successfully extract the AA signal from the considered segment, while the BTD could extract it in 6 out of the 10 runs for patient 9, and 3 out of 10 runs for patient 10.Despite its higher computational complexity, this outcome shows that BTD is capable of extracting the AA signal from segments where the matrix-based techniques here compared could not, pointing out the superiority of BTD as an AA extraction tool.
We can see in Figure 10 that the DF of the estimated atrial  source by BTD does not change considerably over the runs for a given patient, which means that BTD seems to be able to target the AA source always, although with varying accuracy depending on the algorithm convergence.Over the population of observed patients, the DF of the atrial sources estimated by BTD are in the interval of 4.77 to 6.67 Hz, very close to the DF interval of the atrial sources estimated by PCA and RobustICA-f, which lie in 4.77 to 6.44 Hz.

C. Correlation of Kurtosis and SC
In Figure 11 it can be seen how kurtosis and SC of the estimated atrial source are correlated across the observed population of AF patients.It is observed that these parameters are positively correlated, since in general, as long as the values of kurtosis increase, the values of SC also increase.Moreover, it is valid to point that, in the observed population, for high values of kurtosis we will always have high values of SC.However, for some patients it is observed a high value of SC but a low value of kurtosis, as shown in Figure 11.It can also be seen in Figure 11 the difficulties of BTD in extracting the AA from the AF ECG of some patients, as for example, Patient 9, whose points are concentrated in low values of kurtosis and SC.On the other hand, patients with successful AA extractions are observed, as Patient 4, for example, which has its points scattered in high values of kurtosis and SC.Note that some patients have less points than others.This is due to the fact that, as explained before, BTD could not successfully extract the AA in some Monte Carlo runs for some patients.Only the cases where the extraction was performed are analyzed in this final experiment.

VII. CONCLUSION
The present work has evaluated the temporal and interpatient variability of BTD for noninvasive AA extraction in AF ECGs.The first contribution is the analysis of BTD in AF ECG segments with varying time length, showing that this tensor technique presents a satisfactory performance not only in short segments (2.5 s), but also in long segments (15 s).The choice of the size of the observation window is shown not to be very significant regarding the quality of the atrial signal extracted.The only drawback observed is the need of downsampling the long segment by a high factor, which may cause some loss of information.The second contribution of the present work is the performance assessment of BTD in a population of 10 AF patients, confirming the strong variability to initialization of this technique in this particular biomedical application, not only for a fixed patient, but for every patient in an observed population.Nevertheless, with a suitable initialization, BTD can provide a better estimation of the atrial source (as measured by kurtosis and SC) and more information about the AA than matrix-based techniques, since in some cases more than one source with significant AA is recovered.The assessment of the inter-patient variability of BTD in AF ECGs demonstrates that BTD can provide a better AA extraction performance in a population of AF patients when compared to the matrix-based methods PCA and RobustICA-f.Also, it can extract AA in some patients for whom the matrix techniques fail.The correlation of the AA content quality parameters is also presented, showing another way to classify whether a patient poses difficulties in noninvasive AA extraction.
Future works will aim to develop a more robust method, with respect to initialization and the choice of model parameters, to compute BTD and provide better indexes of AA content measurement.Analyzing the cases in which the AA is spread across more than one source is also an objective to be presented in following works, as well as studying the clinical relevance of these additional AA sources.Finally, experiments should be performed in a larger database of AF patients to fully validate this tensor technique and provide results with increased statistical significance.

Fig. 1 .
Fig. 1.Visual representation of the block term decomposition of an arbitrary third-order tensor.

Fig. 2 .
Fig.2.A 15-second segment of an AF ECG recording from the patient used for evaluate the time variability.For concision, only lead II is shown; all 12 leads are processed.

Fig. 4 .Fig. 5 .
Fig.4.Atrial source contribution to lead V1 estimated by BTD, RobustICA-f and PCA, for the 15-second segment, showed in the time domain (in mV).AA signal estimates are vertically shifted for clarity.The power contribution to lead V1 for each technique is also showed.

Fig. 6 .
Fig.6.SC (%) of the AA signals estimated by BTD from the observed AF ECG segments over ten independent runs.

Fig. 7 .
Fig. 7. Kurtosis of the AA signals estimated by BTD from the observed AF ECG segments over ten independent runs.

Fig. 8 .
Fig. 8. Potential atrial sources contribution to lead V1 estimated by BTD for Patient 6. (Top) In the time domain, measured in mV.(Bottom) In the frequency domain, measured in mV/ √ Hz.

Fig. 9 .
Fig. 9. Variation of SC (%) of the atrial source estimated by BTD over 10 tensor factor initializations for the observed population of AF patients.

Fig. 10 .
Fig. 10.Variation of DF (Hz) of the atrial source estimated by BTD over independent runs for the observed population of AF patients.

Fig. 11 .
Fig. 11.Scatter plot showing the correlation between kurtosis and SC of the selected atrial source by BTD, for the observed population of AF patients.A positive correlation can be seen.

TABLE I VALUES
OF SC (%) FOR PCA AND ROBUSTICA-F.FOR BTD, THE MAXIMUM (BTDmax) AND THE MEAN (BTDmean) VALUES OF SC (%) OF TEN INDEPENDENT RUNS IS SHOWN.

TABLE II VALUES
OF SC (%) FOR PCA AND ROBUSTICA-F.FOR BTD, THE MAXIMUM (BTDmax) AND THE MEAN (BTDmean) VALUES OF SC (%) OF TEN INDEPENDENT RUNS ARE SHOWN.