ECG denoising algorithm based on group sparsity and singular spectrum analysis

Biomedical Signal Processing and Control 50 (2019) 62–71

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

ECG denoising algorithm based on group sparsity and singular spectrum analysis Nasser Mourad a,b,∗ a b

Department of Electrical Engineering, Aswan Faculty of Engineering, Aswan University, Aswan 81542, Egypt Department of Electrical Engineering, College of Engineering, Shaqra University, Dawadmi, Riyadh, Saudi Arabia

a r t i c l e

i n f o

Article history: Received 9 September 2018 Received in revised form 1 December 2018 Accepted 7 January 2019 Keywords: ECG denoising Group sparsity SSA Majorization minimization

a b s t r a c t A new algorithm for denoising ECG signals contaminated by additive white Gaussian noise is proposed here. In the proposed algorithm, a clean ECG signal is modeled as a combination of two morphologically different signals. The first signal is a spike-like signal representing the QRS-complex, and, hence, it has the statistical characteristics of a group-sparse (GS) signal, i.e., a sparse signal its non-zero entries tend to concentrate in groups. The second signal, on the other hand, is a smooth signal representing the lowfrequency components of the ECG signal. Accordingly, the proposed algorithm consists of two stages, where in the first stage a new algorithm based on the majorization-minimization (MM) technique is developed to extract the GS signal representing the QRS-complex, while a modified version of the singular spectrum analysis (SSA) technique is utilized in the second stage to smooth the remaining signal. Simulation results on real and simulated ECG data show that the proposed algorithm can be successfully utilised to denoise ECG data. In addition, the proposed algorithm is also shown to produce significantly improved results compared to existing techniques used for performing similar tasks. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction The electrocardiogram (ECG) is the electrical signal reflecting the cardiac activity. ECG is an important tool for the diagnosis of cardiac abnormalities, and is extensively utilized by physicians for the interpretation and identification of physiological and pathological phenomenons. Unfortunately, in real situations, ECG recordings are usually corrupted by various kinds of artifacts and noise. Examples of such artifacts and noise are power line interference, instrumentation noise, ambient electromagnetic (EM) signals picked by the cables, electrode contact noise, baseline drift, motion artifacts, muscle contraction, and ECG signal amplitude modulation with respiration [1]. Unfortunately, the contaminating artifacts and noise severely limits the utility of the measured ECG signal, and, hence, they have to be removed for better clinical evaluation. Fortunately, some of these artifacts, like power line interference and baseline drift, are narrow band, and their bands of frequencies do not overlap with that of the ECG signal. Therefore, these artifacts can be easily corrected using appropriate filters. However, other artifacts and noise, like instrumentation noise and muscle contraction, occupy a

∗ Correspondence to: Department of Electrical Engineering, Aswan Faculty of Engineering, Aswan University, Aswan 81542, Egypt. E-mail address: [email protected] https://doi.org/10.1016/j.bspc.2019.01.018 1746-8094/© 2019 Elsevier Ltd. All rights reserved.

wide band of frequencies, which overlaps with the band of the ECG signal. Therefore, filtering is not appropriate to correct these kinds of artifacts. In this paper, these artifacts and noise components are modeled as white Gaussian noise (WGN), and a new algorithm for denoising ECG signal contaminated by additive WGN (AWGN) is developed. Several approaches have been developed in the literature to perform ECG denoising. These approaches include discrete wavelet transform (DWT) coefficient shrinkage [2,3], empirical mode decomposition (EMD) [4,5] and its extension called ensemble EMD (EEMD) [6,8,7], hybrid EMD-DWT approach [9,10], adaptive Fourier decomposition (AFD) [11], variational mode decomposition (VMD) [12], nonlocal mean (NLM) [13], non-local wavelet transform (NLWT) [14], sparsity-assisted signal smoothing (SASS) [15,16], and many other approaches. However, most of these approaches can be considered as general purpose denoising techniques, and they do not take into consideration the inherent characteristics of the ECG signal. Therefore, each one of these approaches has its own limitations in denoising ECG signals. A brief summary of some of these approaches with their limitations is presented in what follows. The DWT approach is widely used in denoising applications due to its ability to remove AWGN. However, the performance of the DWT approach depends on the selected mother wavelet, as well as on the decomposition level. On the other hand, EMD-based

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

approaches decompose a noisy ECG signal into a set of intrinsic mode functions (IMFs) in adaptive fashion. In these approaches, a denoised ECG signal is usually obtained by discarding the first IMFs which contain the largest amount of noise. Unfortunately, since the QRS-complex occupies a wide band of frequencies, 8–50 Hz or higher [17], removing the lower order IMFs usually introduces severe QRS-complex distortion in the form of R-wave amplitude attenuation [10]. Several approaches have been developed in the literature to reduce the distortion introduced by the conventional EMD-based denoising approach to the QRS-complex [9,10,7]. The main idea of these approaches is to perform signal denoising on the noise-dominant IMFs, rather than discarding them completely. For instance, it was suggested in [9,10] to denoise the IMFs using the DWT thresholding approach. It was shown in [18] that the hybrid EMD-DWT denoising can reduce noise from ECG signals more accurately and consistently in comparison to algorithms in EMD or wavelet domain alone. However, the main shortcoming associated with this approach is that the thresholding step inherent in the DWT approach cannot distinguish between the QRS information and high frequency noise. Therefore, this approach may also result in a distorted QRS-complex in the denoised ECG signal. The authors of [11] proposed to decompose the ECG signal using the adaptive Fourier decomposition (AFD) technique, which has the property of decomposing a signal according to its energy distribution. It was found in [11] that this approach has better performance than both the DWT and the EMD-based methods. On the other hand, the authors of [12] suggested to replace the EMD decomposition by the variational mode decomposition (VMD) approach. Combining VMD with the DWT, the authors of [12] proposed a VMD-DWT algorithm, which was found to outperform the conventional EMD-DWT approach. However, it was also reported in [12] that the performance of the VMD-DWT is inferior to that of the nonlocal mean (NLM) algorithm developed in [13]. In the NLM algorithm proposed in [13], estimates of the samples of a denoised ECG signal are obtained by weighted averaging of the samples of the noisy ECG signal. The applied weights are proportional to the similarity in the neighbourhoods, and they are independent of the temporal locations of the samples. As a result, higher weights are assigned to samples with quite similar neighbourhoods, whereas samples with dissimilar neighbourhoods are given lower weights. However, the main drawback of the NLM algorithm is that it exploits only the correlation among the non-local samples of the ECG signal, while ignoring the correlation among the local samples [14]. To overcome this difficulty, the authors of [14] proposed an algorithm called non-local wavelet transform (NLWT), which exploits the correlation among both local and non-local samples of the ECG signal. Unfortunately, the computation cost of this algorithm is high, and it depends on many parameters that have to be selected wisely to obtain good performance. In this paper, a new ‘two-stage’ approach for denoising ECG data is proposed. In the proposed approach, the ECG signal is modeled as a combination of two morphologically different components. The first component is a spike-like signal representing the QRS-complex, while the second component is a smooth signal representing the remaining low-frequency part, and including the T-wave and the P-wave. Therefore, the QRS-complex is modeled in this paper as a group-sparse (GS) signal, i.e., a sparse signal its nonzero entries tend to concentrate in groups. Accordingly, in the first stage of the proposed approach, an algorithm based on the majorization-minimization (MM) technique [19] is developed to extract the GS signal representing the QRS-complex. After extracting the QRS-complex from the measured ECG signal, a modified version of the singular spectrum analysis (SSA) technique [20–22] is then utilized in the second stage to denoise the remaining signal. The SSA is a nonparametric technique that makes no prior assumptions about the linearity, stationarity or Gaussianity of the

63

Fig. 1. A clean ECG signal u(t) (top) can be modeled as a combination of a groupsparse signal q(t) (middle) and a smooth signal s(t) (bottom).

analyzed data, and, hence, it works with arbitrary statistical processes [21]. The SSA technique decomposes a given signal into a sum of signals, so that each signal in this sum can be identified as either a trend, periodic or quasi-periodic component (perhaps, amplitude-modulated), or noise. Therefore, the SSA technique has useful applications in detrending, smoothing, signal denoising, extraction of periodicities with varying amplitudes and other tasks. See [20–22] and the references therein for more details about the SSA technique. The remaining part of this paper is organized as follows. The proposed algorithm is derived in Section 2. Section 3 presents some experiments to assess the performance of the proposed algorithm. A discussion about the limitations associated with the proposed algorithm, as well as the future study directions to overcome these limitations is presented in Section 4. Finally, concluding remarks are given in Section 5. 2. Methods 2.1. Notations and problem formulation The following notations are used throughout the full text. A bold upper case symbol, e.g., A refers to a matrix; similarly a refers to a column vector; at refers to the tth element of a; {at } is the set of elements in a; diag(a) is a diagonal matrix its main diagonal elements are the entries of a; a discrete-time series z(t), t = 1, . . ., T, is represented (using vector notations) by the vector z ∈ RT , where T is the number of time points; and zp =



T |z |p t=1 t

 1p

is the

p -norm of z. Referring to Fig. 1, it can be observed that a clean ECG signal u(t) (presented by the upper waveform) can be modeled as u(t) = q(t) + s(t), where q(t) and s(t) are two morphologically different signals presented, respectively, by the middle and the lower waveforms in Fig 1 . The signal q(t) is a spike-like signal, and, hence, it can be modeled as a GS signal. This signal occupies a wide band of frequencies (from 8 Hz to 50 Hz or higher [17]) and it represents the QRS-complex of the ECG signal. On the other hand, the signal s(t) is smooth, and it represents the low frequency components (less than 10 Hz) of the ECG signal, i.e., the T-wave and the P-wave. Accordingly, a recorded ECG data x(t) is modeled in this paper using the following equation x(t) = u(t) + n(t) = q(t) + s(t) + n(t),

t = 1, . . ., T,

(1)

where n(t) is a zero mean AWGN, and T is the number of time points (samples). Based on this model, the denoising algorithm proposed

64

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

the concentration of the large values of q(t) within the tth window. In addition, since q(t) is a GS signal, the vector f is sparse, and its nonzero entries are confined within the support of q(t). Accordingly, it is suggested in this paper to get approximate estimate of the QRScomplex using the following optimization problem 1 q˜ = argmin ˜x − q22 + (f ), q 2

Fig. 2. The block diagram of the proposed algorithm.

in this paper for recovering the clean ECG signal u(t) consists of two stages. In the first stage, presented in Section 2.2, an iterative algorithm is developed to extract the GS signal q(t) from the measured signal x(t). While in the second stage, presented in Section 2.3, an algorithm based on the SSA technique is developed to recover the low-frequency component s(t) by denoising the remaining signal xˆ (t):=x(t) − qˆ (t), where qˆ (t) is the estimated QRS-complex obtained in the first stage. The block diagram of the proposed algorithm is presented in Fig. 2. 2.2. Stage 1: Extracting the QRS-complex q(t) As stated before, the QRS-complex q(t) has the statistical properties of a GS signal. Therefore, it can be accurately extracted by an algorithm that maximizes the group-sparsity, rather than the sparsity, of the solution vector. Recently, group-sparsity has found many useful applications in the field of biomedical signal processing [24,25,23]. In this subsection, a new algorithm based on the MM technique is developed to extract a GS signal, representing the QRS-complex, from a measured ECG data. Let’s define the vector f = [f1 . . . fT ]T , where its tth entry ft is defined as the energy of the group of samples of q(t) confined within a short rectangular window (referred to as the “tth window”) extending from ti to tf and centered at the tth sample, that is ft =

tf 

q2j

(2)

j=ti

Normally, the size Lt of the tth window equals (2m + 1) samples, where the value of m > 0 is to be determined. However, To cope with the samples at the beginning and the end of the analyzed signal, the window size in (2) is adjusted by neglecting its portion that does not overlap with the signal. Therefore, ti and tf are calculated as ti = max(1, t − m), tf = min(t + m, T), and, hence, the window size Lt has the following bounds (m + 1) ≤ Lt ≤ (2m + 1). To determine the value of m, it is obvious that m should be selected based on the size of the QRS-complex. Since the duration of a normal QRS–complex approximately equals 100 ms [26], it is recommended to set m = 50 ms. This value of m is utilized in all the experiments presented in this paper. It can be observed from (2) that, with high probability, ft takes large (small) value if many samples within the tth-window have large (small) values. Therefore, ft can be regarded as a measure of

(3)

where (f ) : RT → R is a sparsity promoting regularization function, and  > 0 is a regularization parameter, which is used to manage the trade-off between data fit (measured by ˜x − q22 ) and the group sparsity of the solution vector (measured by (f)). Note that x˜ is a bandpass filtered version of the measured signal x, where the cutoff frequencies of the bandpass filter are selected to maximize the probability of detecting the support of the QRS-complex. See the first remark at the end of this subsection. Therefore, q˜ is not the desired QRS-complex signal q(t). Rather, it represents only an estimate of the QRS-complex q(t) within the band of x˜ . The estimation of the QRS-complex will be described by the end of this subsection. approach to promote the sparsity of f is to select A popular (f) =  f  1 = t |ft |. However, it was shown  in [27–29] that there is a class of nonconvex functions (f) = t=1 (ft ) that can perform better than the convex 1 -norm in encouraging the sparsity of its argument. Examples of such functions that have been extensively utilized in the literature are log (f ; ˛) = log(1 + |f|/˛) and p (f) = |f|p , where ˛ > 0 and 0 < p < 1. However, since (f) is nonconvex, the objective function in (3) has many local minima. Fortunately, this problem can be greatly alleviated by following the majorizationminimization technique (MM) [19,28,29]. Starting from an initial point f0 , the main idea of the MM approach is to locally replace the nonconvex function (f) by a surrogate smooth convex function g(f), which can be easily minimized. The surrogate function g(f) must satisfy the following two constraints: (1) g(f0 ) = (f0 ), and / f0 . Therefore, following the rules of Theorem (2) g(f) > (f), ∀ f = 2 in [28], it is suggested to select g(f) having the following linear expression g(f ) = (f0 ) +  (f0 )(f − f0 ) =  +  (f0 )f,

(4)

where  (f0 ) is the gradient of (f) at f0 , and  := (f0 ) −  (f0 )f0 is a constant that does not depend on f. Starting from an initial solution vector q˜ 0 , and following the MM approach, a new solution vector that reduces the objective function in (3) can be obtained by replacing the nonconvex function (f) in (3) by its convex surrogate function t g(ft ), where g(ft ) is as defined in (4). Following this procedure, and neglecting the constant term  in (4), a solution vector at the (k + 1)th iteration can be obtained by solving the following optimization problem

 1  (fk,t )ft , q˜ k+1 = argmin ˜x − q22 + k q 2 T

(5)

t=1

where k = 0, 1, . . . is the iteration index, and fk,t is the tth sample of fk , which is calculated from q˜ k using (2). In this paper, it is suggested to select (f) = log(1 + f/˛) (Note that f is nonnegative). Therefore, we have  (fk,t ) = 1/(fk,t + ˛k ), and the objective function (5) is reduced to the following simplified form

 ft 1 q˜ k+1 = argmin ˜x − q22 + k , fk,t + ˛k q 2 T

t=1

(6)

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

65

Table 1 Proposed algorithm. ˆ xˆ ] = ExtQRS(x, m, fH1 , fH2 ) [q, 1. Construct x˜ and x¯ using bandpass filters with cutoff frequencies [5, fH1 ] Hz, and [5, fH2 ] Hz, respectively. 2. Initialization: Set q˜ 0 = x˜ and ˛0 = 1. 3. for k = 0 : Kmax (a) Utilize q˜ k and ˛k to calculate  k using (2) and (8). (b) Utilize  k to calculate the value of k using the procedure described in [23]. (c) Utilize  k and k to calculate the weighting vector w k using (11). (d) Smooth out the sharp transitions in w k . (e) Update the vector q˜ k+1 = w k x˜ (f) if q˜ k+1 − q˜ k 2 /q˜ k 2 ≤ q , break, end.





(g) Update ˛k+1 = max 10−4 , ˛k /2 . ¯ and xˆ = x − q. ˆ 4. Calculate: qˆ = w k x,

Substituting the expression of ft from (2) into (6), and rearranging the terms, it is straight forward to show that (6) can be expressed as

 1 k,t q2t , q˜ k+1 = argmin ˜x − q22 + k q 2 T

(7)

t=1

while the second output xˆ = x − qˆ is the remaining part of the ECG data. In addition, the following remarks are in order.

where k,t =

tf  t=ti

1 , fk,t + ˛k

(8)

and ti and tf are as defined in (2). The optimization problem (7) can be written in the following compact form 1 1 q˜ k+1 = argmin ˜x − q22 + k ˙ k2 q22 , q 2

(9)

where ˙ k := diag( k ) is a diagonal matrix, and the the entry of its main diagonal  k is given by  k,t defined in (8). The optimization problem (9) is convex, and has the following closed form solution −1 q˜ k+1 = (I + k ˙ k ) x˜ = W k x˜ = w k x˜ ,

(10)

where I is the identity matrix, denotes the Hadamard prod−1 uct; i.e., element-by-element multiplication, W k = (I + k ˙ k ) = diag(w k ) is a weighting diagonal matrix, with the tth entry of w k given by wk,t =

1 , 1 + k k,t

t = 1, . . ., T.

(11)

and  k,t is calculated from q˜ k using (2) and (8). In order to calculate the value of the regularization parameter k , let’s arrange (11) in the following form wk,t =

(1/k,t ) (1/k,t ) + k

,

Fig. 3. Demonstration of the capability of the proposed ExtQRS algorithm in extracting the QRS-complex q(t) from a noisy ECG signal x(t). The signals from top to bottom represent, respectively, x(t), w k , qˆ (t), and xˆ (t).

t = 1, . . ., T.

which has similar structure as (12) in [23]. Therefore, following the same reasoning as in [23], it is suggested to calculate the value of k by first sorting the values of 1/ k,t in descending order, and then selecting the value of k as the value of 1/ k,t at the corner point of the sorted values. See [23] for more details about this approach. The proposed ExtQRS algorithm (which is used for extracting the QRS-complex) is summarized in Table 1. The algorithm has four inputs and two outputs. The first input is the recorded ECG data x, the second input is the value of m used for calculating fk,t and  k,t , and the last two inputs are the upper cutoff frequencies of the ¯ respectively. band pass filters (BPF) used for constructing x˜ and x, ¯ On the See the first remark below for more details about x˜ and x. other hand, the first output qˆ is an estimate of the QRS-complex,

2.2.1. Remarks 1. Normal QRS-complex signals are approximately confined within the frequency band 8–50 Hz [17]. However, abnormal ventricular conduction introduces notches on the QRS-complex, and, hence, increases its frequency band up to 70 Hz, or even higher [17]. On the other hand, it was reported in [30] that the optimum band of frequencies for detecting the QRS-complex is 8–20 Hz. Therefore, in Step 1 of Table 1 two bandpass signals are constructed. The first signal x˜ is utilized for estimating the support (which is marked by the weighting vector w k ) of the QRS-complex. Therefore, the cutoff frequencies of the BPF used for constructing x˜ are [5, fH1 ] Hz, where fH1 = 20 Hz or slightly higher. The second bandpass signal x¯ is utilized for estimating the QRS-complex. Therefore, this signal is constructed using a BPF with cutoff frequencies [5, fH2 ] Hz, where the value of fH2 depends on whether the QRS-complex is normal (fH2 = 50 Hz) or abnormal (fH2 = 70 Hz). Note that better performance of the proposed algorithm was obtained when the lower cutoff frequencies of the BPFs are set equal 5 Hz, rather than 8 Hz. 2. Usually, the weighting vector w k calculated in Step 3c has sharp transitions (approximately between 0 and 1) at the boundaries of the support of the QRS-complex (see the second waveform in Fig. 3). This sharp transition can introduce discontinuities in the denoised ECG signal at the boundary of the QRS-complex. To overcome this difficulty, it is suggested in Step 3d in Table 1 to smooth the sharp transitions in w k . This step is performed in the proposed algorithm using a 3-tap zero-phase moving average filter, which is implemented using the Matlab function filtfilt . m. 3. It is worth mentioning that the rule played by ˛ in the expression of log (f ; ˛) utilized in the proposed algorithm is similar to that of the exponent p in p (f). Accordingly, for a relatively large value of ˛, log (f ; ˛) becomes very close to the convex 1 -norm, while log (f ; ˛) converges asymptotically to the 0 -norm as the value of ˛ decreases to zero. Therefore, selecting a very small value for ˛ increases the probability of converging to local minima. To reduce this probability, it is suggested in Table 1 to initialize ˛ with a relatively large value ˛0 = 1, and then to reduce its value according to Step 3g.

66

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

¯ the 4. Since the QRS-complex q(t) is confined within the band of x, QRS-complex is estimated in Step 4 by multiplying the estimated ¯ rather than x˜ . weighting vector w k by x, 5. The proposed algorithm terminates if either the number of iterations exceeds Kmax (e.g., Kmax = 100), or the condition at Step 3f is satisfied, where q in Step 3f is a small constant, e.g., q = 10−3 . 6. All signals in this paper are filtered using the Matlab function “eegfilt . m” included in EEGLAB.1 The performance of the proposed ExtQRS algorithm is demonstrated in Fig. 3. The upper waveform in Fig. 3 represents a 5 s simulated ECG signal contaminated by AWGN at 10 dB SNR. The second waveform represents the estimated weighting vector w k . As shown in this figure, the entries of w k almost equal one at the support of the QRS-complex and zero elsewhere. The third and fourth waveforms in Fig. 3 represent the outputs from the ExtQRS algorithm, i.e., qˆ and xˆ , respectively. 2.3. Stage 2: Restoring the low-frequency signal s(t) by denoising xˆ (t) After estimating the QRS-complex qˆ (t), the remaining signal xˆ (t):=x(t) − qˆ (t) is basically the smooth low-frequency signal s(t) contaminated by AWGN. See the lower waveform in Fig. 3. Therefore, in the second stage of the proposed algorithm, an estimate of s(t) is obtained by smoothing xˆ (t). Since s(t) is a low-pass signal, one can use a low-pass filter (LPF) with appropriate cutoff frequency, e.g., fc = 10 Hz. However, a LPF will not remove the noise components within the band of s(t). Another approach is to denoise xˆ (t) using the DWT thresholding technique. Even though the DWT thresholding technique is efficient in denoising smooth signals contaminated by AWGN, this approach depends on some parameters (like the mother wavelet, the decomposition level, the threshold level, and the thresholding approach), which have to be selected wisely in order to get the desired results. Unfortunately, different combination of these parameters produce different results. In addition, a try-and-error procedure has to be followed in order to get the optimum combination of these parameters. In this paper, it is suggested to denoise xˆ (t) using the singular spectrum analysis (SSA) technique summarized in the next subsection. The superiority of the SSA technique, over both the LPF and the DWT thresholding approaches, in smoothing xˆ (t) is demonstrated in Section 3. 2.3.1. Summary of the SSA procedure Denoising xˆ (t) using SSA consists of the following five steps. Step 1: Construct a trajectory matrix X ∈ RL×K from the noisy signal xˆ (t), where L ≤ T/2 is called the ‘embedding’ parameter, and K = T − L + 1. The dependency of the denoising algorithm on the value of L is demonstrated in Section 3. The trajectory matrix has the following structure

⎡

xˆ 1

⎢ xˆ ⎢ 2 X=⎢ ⎢ . ⎣ .. xˆ L

xˆ K

xˆ 2

xˆ 3

...

xˆ 3

xˆ 4

...

.. .

.. .

..

.

.. .

xˆ L+1

xˆ L+2

...

xˆ T

⎤

⎥ xˆ K+1 ⎥ ⎥ ⎥ ⎦

(12)

Note that the trajectory matrix X is a Hankel matrix, i.e., it has equal elements on the ‘antidiagonals’ i + j = const. Note also that X can be expressed as X = S + N, where S and N are the trajectory matrices of s(t) and n(t), respectively.

1

This is a well known toolbox used for analyzing EEG data, and it is freely available at https://sccn.ucsd.edu/eeglab/.

Step 2: Construct the matrix R = XXT , and calculate its eigenvalue decomposition R = UEUT , where U is the matrix of eigenvectors, and E = diag(e) is a diagonal matrix containing the corresponding eigenvalues. It is assumed that e1 ≥ e2 ≥ . . . ≥ eL ≥ 0. Step 3: Split the eigenvectors U into two submatrices Us and Un , that is U = [Us Un ], where Us and Un represent the signal and the noise subspaces, respectively. The splitting step is the most critical step in the denoising process, and is presented shortly. Step 4: An estimate of the signal matrix S can then be obtained by projecting X into the column space of Us , that is Sˆ = U s U Ts X. Note that, generally, Sˆ is not a Hankel matrix. Step 5: Construct the one-dimensional smooth signal sˆ(t) by ˆ averaging the antidiagonals of the matrix S. 2.3.2. Estimating the signal subspace As can be observed from the above procedure, the key step in denoising a signal using SSA is the identification of the signal subspace Us . Conventionally, the signal subspace Us is constructed from the r–leading eigenvectors of U, i.e., Us = [u1 . . . ur ]. This approach relies on the separability assumption, i.e., the assumption that the smallest eigenvalue associated with the signal component is larger than the largest eigenvalue related to the noise component. Therefore, the value of r is usually determined as the index of the eigenvalue at which the noise floor starts [21,22]. However, this approach may fail when the separability assumption is violated. To overcome this difficulty, it is suggested to identify the signal subspace Us using the smoothness of the signal s(t). To accomplish this goal, let’s define the singular value decomposition (SVD) of the trajectory matrix X as X = UCV T =

L  j=1

cj uj v Tj =

L 

(13)

Xj,

j=1

where the columns uj and v j of U and V are respectively called the left and right singular vectors of X, C = diag(c) is a diagonal matrix containing the singular values along its main diagonal, and   X j :=cj uj v Tj

is a set of rank–one matrices. Since X = S + N, it is sug-

 

gested to construct S using only the set of matrices X j which have smooth rows. Note that only the rows are considered because we set of matrices. usually have K  L. Let Js refer to the indices of this  Therefore, an estimate of S can be expressed as Sˆ = X j . Since j∈J s

the rows of Xj are proportional to v Tj , it is suggested to utilize the smoothness of the columns of V to identify Js . The smoothness of a vector v can be quantified using the quantity z:=Dv2 , where D is the first-difference matrix. The smoother the vector v, the smaller the value of z. Let’s define the vector z ∈ RL , where its jth entry is defined as zj = Dv j 2 , and v j is the jth column of V. Assume without loss of generality that z1 ≤ z2 ≤ . . . ≤ zL . By plotting the entries of z, it is expected to get an abrupt change in the slope of the plotted curve at the boundary between the values of z associated with the singular vectors related to the signal and those related to the noise. See the upper panel in Fig. 4. Let jb refer to the index of this boundary point. Accordingly, Js is estimated as Js = {1, . . ., jb }. To estimate the value of jb , let a ∈ RL−2 be a vector containing the slopes calculated at each point of z, where the slope at the point zj is defined as aj = (zj+1 − zj−1 )/2 ∀ 1 < j < L. The value of jb is then estimated as the index of the largest value in a. See the lower panel in Fig. 4. It is worth mentioning that, in order to compute zj , the full SVD of the trajectory matrix X does not have to be computed. Rather, the relation between the EVD and the SVD is utilized, i.e., v j is calculated using the following equation vj =

1

 X T uj , j = 1, . . ., L ej

(14)

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

67

Fig. 5. Estimating the optimum value Lopt of the imbedding parameter L.

Fig. 4. Identification of the indices of V associated with the signal component.

where ej and uj are, respectively, the jth eigenvalue and the associated jth eigenvector of the matrix R calculated in Step 2. 3. Simulations and results In this section, four different examples are presented to examine the performance of the proposed algorithm. In all the presented examples, either simulated or real ECG signals are considered. The simulated ECG signals are generated using the algorithm from [31], and are sampled at 256 Hz. The real ECG signals, on the other hand, are taken from the MIT-BIH arrhythmia data base [32]. This database contains 48 different ECG records, and each record contains two channel signals with 30 min duration, and they are sampled at 360 Hz. Only the signals recorded at the first channel are considered here. To quantify the performance of various algorithms, the improvement in the signal-to-noise-ratio (SNRimp ) in dB is utilized as a metric. The SNRimp is computed using the following equation [10,7]

T

− u(t))2 t=1 (x(t) , 2 T ˆ u(t) − u(t) t=1

SNRimp [dB] = 10log10 





(15)

where x(t) is the measured ECG signal, T is the number of samˆ is its estimate. In addition, ples, u(t) is the clean ECG signal and u(t) the various parameters used with the proposed algorithm are as follows: m = 50 ms, fH1 = 20 Hz, fH2 = 35 Hz (70 Hz) for the simulated (real) ECG data, and, according to the results of 3.1Example 1 presented below, L = 2.2 s. 3.1. Example 1: Estimating the optimum value of L The performance of the SSA procedure depends on the value of the embedding parameter L, which controls the size of the trajectory matrix X ∈ RL×K . In order to find the optimum value of L (denoted as Lopt ) as a function of the input SNR, the following experiment is conducted. The values of the input SNR considered in this experiment are 5, 10, 15, and 20 dB. For each value of the input SNR, a simulated 10 s ECG data segment u(t) is first generated, and a noisy ECG signal x(t) is then constructed by adding a WGN at the specified SNR. The constructed noisy signal x(t) is then applied to the ExtQRS algorithm to extract the QRS-complex qˆ (t). After extracting the QRS-complex, the SSA procedure is followed to denoise the remaining signal xˆ (t). As described in Section 2.3.1, the trajectory

matrix X ∈ RL×K is constructed in the first step of the SSA procedure. The values of L considered in this experiment varies from 0.2 sec to 4 s in steps of 0.2 s. For each value of L, a denoised signal sˆ(t) is obtained using the SSA procedure proposed in Section 2.3, and the ˆ = qˆ (t) + sˆ(t). The SNRimp denoised ECG signal is constructed as u(t) parameter is then calculated using (15). This procedure is repeated 100 different times for each combination of L and the input SNR, where in each time a new ECG data segment and a new WGN are generated. The average of the 100 runs is then calculated, and the result is presented in Fig. 5. The value of Lopt is defined as the value of L at which the SNRimp metric is maximum. The value of Lopt associated with each value of the input SNR is marked by a black star in Fig. 5. As shown in Fig. 5 Lopt = 2.2 s for all the considered values of the input SNR. This value of L is utilized in all the remaining experiments. Note that the estimated value of Lopt is associated with data segments having duration 10 s. Longer data sets can be divided into contiguous segments of durations 10 s. For shorter segments, the value of Lopt can be scaled proportionally. 3.2. Example 2: Estimating the signal subspace As stated in Section 2.3.1, the most critical step in the SSA procedure is the splitting step (Step 3), in which the signal subspace Us is identified. In this example, a comparison between the conventional approach and the approach suggested in Section 2.3.2 for performing this step is performed. The comparison is made in terms of the SNRimp metric as a function of the input SNR. As in 3.1Example 1, the considered values of the input SNR are 5, 10, 15, and 20 dB. For each value of the input SNR, a noisy 10 s ECG data segment x(t) is constructed as in 3.1Example 1, and an estimate qˆ (t) of the QRS-complex is obtained using the ExtQRS algorithm. The remaining signal xˆ (t) is then denoised using the SSA procedure. The splitting step (Step 3 in the SSA procedure) is then performed using either the conventional approach or the proposed approach. In each case, a denoised signal sˆ(t) is obtained, and an estimate of ˆ the ECG signal is constructed as u(t) = qˆ (t) + sˆ(t). For each signal ˆ obtained from each method, the SNRimp metric is then calcuu(t) lated using (15). This procedure is repeated 100 times for each value of the input SNR, where in each time a new ECG data segment and a new WGN are generated. The average of the 100 runs is then calculated, and the result is presented in Fig. 6a. As shown in this figure, the proposed approach outperforms the conventional approach for all the considered values of the input SNR. The performance degradation of the conventional approach as the value of the input SNR increases can be explained as follows. As the value of the input SNR increases, some eigenvectors that represent fine details of the clean

68

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71 Table 2 Average convergence time (s) for the results in Fig. 7(a).

GSSSA GSWD GSLPF DWT EMD-DWT NLM AFD SASS

Fig. 6. Comparison between the conventional approach and the proposed approach in identifying the signal subspace.

ECG signal s(t) start to have small eigenvalues. Therefore, these vectors are erroneously assigned by the conventional approach to the noise subspace, rather than the signal subspace. A demonstrating example is presented in Fig. 6b. The upper waveform in Fig. 6b represents a 1 s noisy ECG data segment. The SNR of this signal equals 15 dB. The middle and lower waveforms in this figure represent the denoised ECG signals obtained by the proposed algorithm when the proposed and the conventional approaches are utilized, respectively. The two signals are superimposed on the original ECG signal s(t), which is presented by the black waveforms. As shown in this figure, there is a mismatch between the lower waveform and the original signal, especially over the interval of the QRS-complex. No such distortion is observed in the middle waveform. 3.3. Example 3: Quantitative comparison with other methods In this example, a quantitative comparison between the proposed GSSSA algorithm and the DWT, the hybrid EMD-DWT, the AFD, the NLM, and the SASS algorithms is presented. The comparison is made in terms of the average SNRimp (dB) metric calculated at different values (5, 10, 15, and 20 dB) of the input SNR. In addition, to demonstrate the efficacy of utilizing the SSA procedure in the second stage of the proposed algorithm, two other versions (GSWD and GSLPF) of the proposed algorithm are also considered. In these two versions, the second stage in the proposed algorithm is performed

5 dB

10 dB

15 dB

20 dB

3.627 1.087 1.083 0.238 1.604 5.467 37.624 13.749

3.480 0.946 0.944 0.451 2.993 6.286 46.985 12.370

3.449 0.890 0.888 0.509 3.399 5.742 45.094 13.278

3.416 0.873 0.870 0.510 3.339 5.592 45.069 12.502

using either WD (for the case of GSWD) or low pass filtering with cutoff frequency fc = 10 Hz (for the case of GSLPF). The value of fc is selected based on the fact that the frequency contents of the T-wave and the P-wave extend respectively to 5 Hz and 10 Hz [33]. On the other hand, the WD step in the GSWD algorithm is performed using the Matlab function “wden.m” with the following parameters. The utilized mother wavelet is “sym7”, and the signal is decomposed to the third level. “Soft” thresholding is used with the “wden.m” function, with the threshold value selected using the principle of Stein’s Unbiased Risk. In addition, rescaling of the wavelet coefficients is performed using level dependent estimation of level noise. These parameters are also used with the DWT and the EMD-DWT approaches. Moreover, simulated and real ECG data were also considered in this example to demonstrate the dependency of various algorithms on the type of the analyzed data. For the case of simulated ECG data, a noisy 10 sec ECG data segment at a specified value of the input SNR is constructed as in 3.1Example 1. The constructed noisy ECG data segment is then applied to various algorithms, and the SNRimp (dB) metric associated with each algorithm is then calculated. This procedure is repeated 100 different times for each value of the considered SNR, where in each run a new ECG data segment and a new Gaussian noise are generated. On the other hand, for the case of denoising real ECG data, the real ECG data is taken from the first channel in record 103 of the MIT-BIH arrhythmia data base. The data is first filtered using a high-pass filter with cutoff frequency fc = 0.5 Hz to remove the baseline wandering. The data is 30 min long, and is divided into contiguous segments of duration 10 s, which results in a total of 180 data segments. Unfortunately, the original ECG data is noisy, and, hence, the value of the input SNR cannot be adjusted accurately. Therefore, to construct a 10 s noisy ECG data segment with accurate value of the input SNR, the data segment is first denoised using the proposed algorithm before adding the WGN at the specified value of the input SNR. The remaining steps are similar to the ones followed with the simulated ECG data, where the SNRimp (dB) metric in this case is calculated as the average over the 180 data segments. The results are shown in Fig. 7. As shown in Fig. 7, and for both types of the ECG data and at all values of the considered input SNR, the general-purpose denoising algorithms (the DWT and the EMD-DWT) have the worst performance among all the considered algorithms, while the proposed GSSSA algorithm, on the other hand, has the best performance. In addition, it is also obvious from Fig. 7 that the SSA technique is more efficient in performing the second stage of the proposed algorithm than either the low pass filtering or the WD approaches. It can also be observed from Fig. 7 that the GSLPF performs better than the GSWD for the cases of relatively low SNR (SNR ≤ 10 dB), while the opposite is true for the case of relatively high SNR (SNR ≥ 15 dB). It can also be observed from Fig. 7 that all the considered algorithms, except the NLM, perform better for the case of simulated ECG data shown in Fig. 7(a) than the case of real ECG data shown in Fig. 7(b). The average convergence time associated with various algorithms, for the case of denoising simulated ECG data shown in Fig. 7(a), is presented in Table 2. Similar results were also obtained,

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

69

Fig. 8. Performance of the proposed algorithm for various data sets from the MITBIH arrhythmia data base.

SNR, the performance of the proposed algorithm is consistent with various data sets. 3.4. Example 4: Qualitative comparison with other methods

Fig. 7. Comparison between the proposed GSSSA algorithm and three other algorithms in terms of the mean SNRimp (dB) at different values of the input SNR using (a) Simulated ECG data, and (b) record 103 of the MIT-BIH arrhythmia database.

but not presented, for the case of real ECG data shown in Fig. 7(b). As can be observed from Table 2, the DWT has the shortest convergence time. However, referring to Fig. 7(a), it can also be observed that the DWT has the worst performance among all the considered algorithms. For the proposed GSSSA algorithm, it can be observed that it has the shortest convergence time compared to the competing AFD, NLM, and SASS algorithms. In addition, comparing the convergence time of the GSSSA with that of its two modified versions GSWD and GSLPF it can be inferred that most of the computation cost of the GSSSA algorithm resides in its second stage performed using the SSA technique. It is worth mentioning that the convergence time of the SASS algorithm can be substantially reduced by fixing the value of the regularization parameter used in the algorithm, rather than iteratively estimating its value using the procedure suggested in [15]. However, this may result in performance degradation. Finally, to check the dependency of the performance of the proposed GSSSA algorithm on the analyzed ECG data, the above procedure is repeated with records 100, 113, 115, 122, 124, and 231 of the MIT-BIH arrhythmia data base. The result is shown in Fig. 8. As shown in this figure, and for all the considered values of the input

In this example, a qualitative comparison between the proposed GSSSA algorithm and the AFD, the NLM, and the SASS algorithms considered in 3.3Example 3 is presented. Three different real ECG data segments taken from the MIT-BIH arrhythmia database are considered in this example. The three data segments are contaminated by real EMG artifacts at different levels of the input SNR. The first data segment has a duration of 5 s, and is taken from record 101. This data segment is presented by the upper waveform in Fig. 9(a), and it represents the case of relatively high input SNR. The cases of moderate and low input SNR are presented, respectively, by the upper waveforms in Fig. 9(b) and (c). These two 10 s data segments were taken, respectively, from records 228 and 104. As can be observed from Fig. 9, and for all the considered cases, the proposed GSSSA algorithm outperforms the three other algorithms considered in this example. More specifically, the ECG data segments were denoised successfully by the proposed algorithm, while residual noise can be observed in the data segments corrected by the other methods. On the other hand, the high performance of the proposed algorithm in correcting the noise is associated with preserving all features of the ECG data. The other methods, on the contrary, introduce distortions to the main features of the corrected ECG data. For instance, and generally speaking, the SASS algorithm introduces R-wave amplitude attenuation, while both the AFD and the NLM algorithms introduce over smoothness to the corrected data. 4. Discussion As shown in Section 3, the proposed GSSSA algorithm outperforms all the competing algorithms in terms of the accuracy and the convergence time. However, in this section some limitations associated with the proposed algorithm, as well as future study directions to overcome these limitations, are presented. The first limitation is associated with the sequential estimation of the two components (s(t) and q(t)) of the ECG signal u(t). It has been observed that, if an error occurred in the first stage, in which q(t) is estimated, this error may not be corrected in the second stage. As a demonstrating example, assume that we have one QRS-complex much shorter than the other QRS-complex signals within the analyzed data segment. In such a situation, the shortest QRS-complex may be missed by the

70

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71

the impact of this limitation on the overall response of the proposed algorithm is small. One approach to overcome this limitation is to estimate s(t) and q(t) simultaneously, rather than sequentially. The second limitation associated with the proposed GSSA algorithm is using a fixed value for m, which should be selected base on the duration of the QRS-complex. The impact of this limitation can be observed when the analyzed ECG data segment is contaminated by a strong nonstationary muscle artifact. In this case, if the value of m was not selected properly, then part of the strong muscle artifact may be captured by the ExtQRS algorithm, and it will appear in the estimated qˆ (t) signal. In this paper, this problem was greatly alleviated by first utilizing the filtered signal x˜ to estimate the weighting vector w k , and then estimating q(t) by multiplying the estimated ¯ The performance can be further improved if weighting vector by x. the value of m was selected adaptively, and/or the quasi periodicity of the QRS-complex is utilized in the suggested algorithm. As has been observed from Table 2, even though the proposed GSSSA algorithm has a relatively short convergence time compared to other algorithms, most of the computation cost of the proposed GSSSA algorithm resides in its second stage performed using the SSA technique. Therefore, the computation cost of the proposed GSSSA algorithm can be further improved by replacing the SSA technique by another technique that converges in a relatively short time without degrading the overall performance. 5. Conclusion In this paper, a new two-stage algorithm for denoising ECG signals contaminated by AWGN was proposed. In the proposed algorithm, a clean ECG signal was modeled as a combination of two morphologically different signals. The first signal is the QRScomplex, which has the characteristics of a GS signal, while the second signal is a smooth signal representing the low-frequency components of the ECG signal. Accordingly, in the first stage of the proposed algorithm a new algorithm was developed to extract the GS signal representing the QRS-complex. After extracting the QRScomplex, the remaining signal was smoothed in the second stage using a modified version of the SSA technique. Simulation results on real and simulated ECG data showed that the proposed algorithm can be successfully utilised to denoise ECG data. In addition, the proposed algorithm was also shown to produce significantly improved results compared to existing techniques used for performing similar tasks. References

Fig. 9. Qualitative comparison between the proposed GSSSA algorithm and three other algorithms in correcting real ECG data contaminated by real EMG noise at (a) high SNR, (b) moderate SNR, and (c) low SNR.

ExtQRS algorithm, and the SSA algorithm in the second stage may capture only the low-frequency components of the missed QRScomplex. This will result in a reduction in the R-peak of the shortest QRS-complex in the denoised ECG signal. Even though, this scenario was rarely observed, the results presented in Section 3 show that

[1] G.M. Friesen, T.C. Jannett, M.A. Jadallah, A comparison of the noise sensitivity of nine QRS detection algorithms, IEEE Trans. Biomed. Eng. 37 (1) (1990) 85–98. [2] P.M. Agante, J. de Sa, ECG noise filtering using wavelets with soft-thresholding methods, Comput. Cardiol. 26 (1999) 535–538. [3] M. Alfaouri, K. Daqrouq, ECG signal denoising by wavelet transform thresholding, Am. J. Appl. Sci. 5 (3) (2008) 276–281. [4] M. Blanco-Velasco, B. Weng, K.E. Barner, ECG signal denoising and baseline wander correction based on the empirical mode decomposition, Comput. Biol. Med. 38 (2008) 1–13. [5] S. Pal, M. Mitra, Empirical mode decomposition based ECG enhancement and QRS detection, Comput. Biol. Med. 42 (1) (2012) 83–92. [6] K.M. Chang, Arrhythmia ECG noise reduction by ensemble empirical mode decomposition, Sensors 10 (2010) 6063–6080. [7] P. Nguyen, J. Kim, Adaptive ECG denoising using genetic algorithm-based thresholding and ensemble empirical mode decomposition, Inf. Sci. 373 (2016) 499–511. [8] A. Komaty, A.-O. Boudraa, B. Augier, D. Daré-Emzivat, EMD-based filtering using similarity measure between probability density functions of IMFs, IEEE Trans. Instrum. Meas. 63 (2014) 27–34. [9] N. Li, P. Li, An improved algorithm based on EMD-wavelet for ECG signal de-noising, Proc. Int. Joint Conf. on Computational Sciences and Optimization (2009) 825–827. [10] M.A. Kabir, C. Shahnaz, Denoising of ECG signals based on noise reduction algorithms in EMD and wavelet domains, Biomed. Signal Process. Control 7 (2012) 481–489.

N. Mourad / Biomedical Signal Processing and Control 50 (2019) 62–71 [11] Z. Wang, F. Wan, C. Wong, L. Zhang, Adaptive Fourier decomposition based ECG denoising, Comput. Biol. Med. 77 (1) (2016) 195–205. [12] S. Lahmiri, Comparative study of ECG signal denoising by wavelet thresholding in empirical and variational mode decomposition domains, Healthc. Technol. Lett. 1 (3) (2014) 104–109. [13] B. Tracey, E. Miller, Nonlocal means denoising of ECG signals, IEEE Trans. Biomed. Eng. 59 (9) (2012) 2383–2386. [14] S.K. Yadav, R. Sinha, P.K. Bora, Electrocardiogram signal denoising using non-local wavelet transform domain filtering, IET Signal Process. 9 (1) (2015) 88–96. [15] I.W. Selesnick, et al., Sparsity-assisted signal smoothing, in: R. Balan (Ed.), Excursions in Harmonic Analysis 4, Birkhauser Basel, 2015, pp. 149–176. [16] I.W. Selesnick, Sparsity-Assisted Signal Smoothing (Revisited), ICASSP2017, March, New Orleans, LA, USA, 2017. [17] L.G. Tereshchenk, M.E. Josephson, Frequency content and characteristics of ventricular conduction, J. Electrocardiol. 48 (6) (2015) 933–937. [18] M.A. Kabir, C. Shahnaz, Comparison of ECG signal denoising algorithms in EMD and wavelet domains, IJRRAS 11 (3) (2012) 499–516. [19] D.R. Hunter, K. Lange, A tutorial on MM algorithms, Am. Stat. 58 (1) (2004) 30–37. [20] N. Golyandina, A. Zhigljavsky, Singular Spectrum Analysis for Time Series, Springer Briefs in Statistics, 2013, http://dx.doi.org/10.1007/978-3-64234913-3 2. [21] H. Hassani, R. Mahmoudvand, Singular Spectrum Analysis Using R, Palgrave Advanced Texts in Econometrics, 2018, http://dx.doi.org/10.1057/978-1-13740951-5. [22] T.J. Harris, H. Yuan, Filtering and frequency interpretations of Singular Spectrum Analysis, Physica D 239 (2010) 1958–1967.

71

[23] N. Mourad, Automatic correction of short-duration artifacts in single-channel EEG recording: a group-sparse signal denoising algorithm, IET Signal Process. 12 (4) (2018) 549–557. [24] Y. Jiao, Y. Zhang, X. Chen, E. Yin, J. Jin, X. Wang, A. Cichocki, Sparse group representation model for motor imagery EEG classification, IEEE J. Biomed. Health Inform. (2018), http://dx.doi.org/10.1109/JBHI.2018.2832538. [25] Y. Zhang, C.S. Nam, G. Zhou, J. Jin, X. Wang, A. Cichocki, Temporally constrained sparse group spatial patterns for motor imagery BCI, IEEE Trans. Cybern. (2019) (accepted for publication). [26] G.D. Clifford, F. Azuaje, P.E. McSharry, Advanced Methods and Tools for ECG Data Analysis, Enggineering in Medicine & Biology, London, 2006. [27] R. Chartrand, Exact reconstruction of sparse signals via nonconvex minimization, IEEE Signal Process. Lett. 14 (October (10)) (2007). [28] N. Mourad, J.P. Reilly, Minimizing nonconvex functions for sparse vector reconstruction, IEEE Trans. Signal Process. 58 (7) (2010). [29] N. Mourad, Advances in Sparse Signal Analysis With Applications to Blind Source Separation and EEG/MEG Signal Processing, Ph.D. dissertation, McMaster University, Hamilton, ON, Canada, 2009. [30] M. Elgendi, M. Jonkman, F. DeBoer, Frequency bands effects on QRS detection, in: 3rd Int. Conf. on Bio-inspired Sys. and Sig. Proc, BIOSIGNALS 2010 – Valencia, Spain, 20–23 January, 2010. [31] P.E. McSharry, G.D. Clifford, L. Tarassenko, L.A. Smith, A dynamical model for generating synthetic electrocardiogram signals, IEEE Trans. Biomed. Eng. 50 (3) (2003) 289–294. [32] The MIT-BIH arrhythmias database [Online]. Available from: http://physionet. org/physiobank/database/mitdb/. [33] D.E. Upasani, R.D. Kharadkar, Automated ECG diagnosis, IOSR J. Eng. 2 (5) (2012) 1265–1269.