Sunday, September 11, 2016

QRS detection in 12-lead ECG based Multiresolution Singular Value Decomposition

This is a summary of work "Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução".

Multiresolution Singular Value Decomposition (MRSVD) 

The Singular Value Decomposition (SVD) is a method to decompose a matrix in three other matrices: two orthogonal and one diagonal. Let $X$ $m×n$ be a real-valued matrix, with $m ≤ n$, then the SVD of $X$ may be written [4]:

$X = U Γ V^{T}$ (1)

where $U$ $m×m$ and $V$ $n×m$ are orthogonal matrix and $Γ$ $m×m$ a diagonal matrix.

The columns of $U$, called left singular vectors, are eigenvectors of $XX^{T}$ and the columns of $V$, called right singular vectors, are eigenvectors of $X^{T}X$. The diagonal entries of $Γ$ are the singular values, square root of the eigenvalues of $XX^{T}$ or $X^{T}X$ and are ordered so that $γ 1 ≥ γ 2 ≥ · · ·≥0 $. The orthogonality of $U$ also allows to write equation (1) as $X̂ = ΓV T$ , where

$X̂ = U^{T}X$ (2)

The Principal Component Analysis (PCA) is a method related to SVD which decomposes a data matrix in uncorrelated principal components. The columns of the decomposition matrix are formed by the eigenvec-tors of the covariance matrix. Let $X$ be a data matrix and $X̄$ is the same matrix with the mean removed. The eigendecomposition of the scatter matrix, estimation of the covariance matrix, may be written $X̄X̄^{T} = UΛU^{T}$ [4], where the columns of $U$ are eigenvectors associated with eigenvalues of $Λ$. The principal components are given by equation $Y = U^{T}X̄$, which is the same as (2), except for the mean removed. This implies that, when the mean of $X$ is removed before the SVD, the principal components of $X̂$ are uncorrelated. This way, the first principal component $φ = X̂(1, ·)$ is the approximation signal, and the other $ψ = X̂(2, ·)$, the detail component, where $X̂(n, ·)$ is n-th row of matrix $X̂$. Relating to the structure of a two-channel filter bank, $U$ represent the filtering process and $X̂$ a sub-band decomposition of $X$ [4].

In order to implement a multiresolution SVD, for the one-dimensional case, [4] propose the construction of a signal $X = [x(1) · · · x(n)]$, being n divisible by $2^L$ for any integer $L ≥ 1$, as a data matrix with the first row containing the odd-numbered entries, and the other ones containing the even-numbered entries. Therefore, the first resolution level the data matrix becomes

$X_{1} = (\table x(1), ... , x(n-1); x(2), ... , x(n))$ (3)

The next step is the removing of the mean, multiplying $X_1$ by $H_1$ , be $H_l = I_l − ( 1/l )e_l e_{l}^{T}$ for the $l$th level, where $e_l$ is the $l × 1$ vector containing all unit entries. For the decorrelation filter $U$ the eigen-decomposition of the scatter matrix is computed. The components in subband are obtained by multiplying of this filter by the data matrix with the mean removed (equation 7).

The following equations define the MRSVD for $l = 1, · · · , L − 1$ levels, $n_l = n/2^l$ and $φ_0 = X$:

$X_{l} = (\table φ_{l-1} (1), ... , φ_{l-1} (n-1); φ_{l-1} (2), ... , φ_{l-1} (n))$ (4)
$ X̄ _{l}= X _{l} H_{n_l} $ (5)
$ X̄ _{l} X̄ _{l}^{T} = U_l Γ_{l}^{2} U_{l}^{T}$ (6)
$X̂_{l} = U_{l}^{T} X̄ _{l}$ (7)
$ φ_{l} = X̂_{l}(1, cdot)$ (8)
$ ψ_{l} = X̂_{l}(2, cdot)$ (9)

Proposed Method 

The sum this detail components provides a signal whose QRS complexes are in highlighted, however due the oversampling operation the location of the start point and end in that sum is not same of the QRS complex in the original signal, therefore an oversampling operation is need. Besides it is appropriate to calculate the absolute value of the oversampled detail components, to avoid the calculation of two thresholds.

The ECG signal after the pre-processing stage is:

$ψ_{sum} = |(↑2^{2})ψ_2| + |(↑2^{3})ψ_3|$ (11)

In the decision stage, the threshold is calculated by equation:

$thr = 0.3 max(ψ_{sum})$ (12)

and checked if the amplitude of $ψ_sum$, for each sample, is higher than $thr$. In the affirmative case the beginning of the QRS complex is found. After detecting the starting point, if the amplitude of $ψ_sum$ is less than $thr$, then is set the end of the QRS complex. However due the oversampling operation (equation (11)) have zeros between samples and this causes the algorithm finds the end of one QRS complex when in reality it is not the end.

The solution is using the physiological information of the heart to avoid the premature detection of end of the QRS complex, based on the fact which depolarization of the ventricles lasts from 0.06 to 0.10 (or 0.12) seconds [11]. Thus it cannot be found the end of the QRS complex before $0.06 Fs$ samples, where $Fs$ is the sampling rate. Let $qrs_ s$ and $qrs_e$ be the location of the starting samples and end, respectively, of each QRS complex, so the difference between them must satisfy the relation:

 $| qrs_s - qrs_e| >= 0.12 Fs$ (13)

It is clear that the limitation in equation (13) can returns an interval greater than the QRS interval, perhaps until ST segment. Nevertheless the maximum point it will still be a peak R, because in $ψ_sum$ the waves of low frequencies, waves T and P, have little energy, soon their amplitudes are attenuated. After found the intervals of the QRS complexes, a search is made, in the original signal, by the maximum point corresponding to the R peak.

The decision stage is executed to check if the R peak detected is authentic or not. The initial verification is relating the distance of the first R peak, which must occur after the depo-larization of the atria, i.e, after the P wave. Let $r(t)$ be the vector which stores the R peaks found, where $t = (t_1 , · · · , t_n )$ is the vector with the localization of all R peaks. For first R peak, the following equation must be satisfied:

$t_1 > 0.2 Fs$ (14)

When the algorithm has detected at least one R peak, it executes the next verifying based in the distance between two R peaks that is related to rate heart $RR = 60/HR$ [6], where $HR$ is a rate heart, measured in bpm (beats per minute) and $RR$ seconds. The maximum rate heart was obtained by experiments of [12] and is expressed by equation $HR_{max} =208−0.7a$, where $a$ is the patient age. Taking ${HR_max} = 208$, since there is no information about the age of the patients in selected database, the minimum samples amount between two R peaks is

$RR = 60/208 Fs$ (15)

Therefore, if RR, than the minimum, it should be chosen the one with the largest amplitude, removing the other which does not satisfy the criterium. Otherwise, both peaks are stored as correct. After detected all R peaks, its analyzed, in the correction stage, the need to search other R peaks. When exists a R peak 30% larger than the most of the peaks the detection fails because the threshold, equation (12), will be above the outhers. To circumvent this inconvenience a second threshold should be computed, removing the highest R peak and computing a new threshold for the resulting signal.

To decide when this calculation should be performed, it is considered the minimum amount of R peaks for any ECG signal, which is expressed by equation

$R_min = τ HR_min / 60$ (16)
where $τ$ is the signal duration, in seconds, and $HR_min = 35$ [5] is a minimum rate heart. The manner as this is done in practice, is by copying the original signal and replacing the of QRS interval by a null vector, and compute a new threshold.

Python scripts

 file MSVD.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import numpy as np

class MSVD:
    
    def __init__(self, signal):
        
        self.signal = signal
    
    '''
    Decomposes the signal at a certain level
    If level_dec = 0, decompose in max level
    '''    
    def dec(self, level_dec, mean_removal = False):
    
        phi = signal
        
        if level_dec == 0:
            level_dec = self.levelmax(self.signal)
    
        for l in range(1, level_dec + 1):
            N = phi.size
            even = phi[0:N:2]
            odd = phi[1:N:2]
            X = np.vstack((even , odd))
            
            if (mean_removal):
                Xbar = X - X.mean()
            else:
                Xbar = X
            
            U, S, Vt = np.linalg.svd(Xbar)
            Xhat = np.dot(U.T , Xbar)
            phi = Xhat[0]
            psi = Xhat[1]
    
        return phi, psi, U
        
    '''
    Returns max level decomposition
    '''
    def levelmax(self, signal):
        return int(np.log(signal.size) / np.log(2)) - 1

------------------------------------------------------------------------------------------------------------------

file QRSdetect.py

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import numpy as np
from MSVD import MSVD
from multirate import multirate 


class QRSdetect:
    
    def __init__(self, ecg, Fs, HR_min = 30):
        
        #taxa de amostragem
        self.Fs = Fs
        #sinal a ser analisadao
        self.ecg = ecg
        #batimentos minimo
        self.HR_min = HR_min
            

    def detect(self, level_dec = 2):


        msvd = MSVD(self.ecg)
        phi, psi, U = msvd.dec(level_dec)
        mt = multirate(psi)
        signal0 = mt.upsample(2**level_dec, -2**level_dec)
        
        msvd = MSVD(self.ecg)
        phi, psi, U = msvd.dec(level_dec + 1)
        mt = multirate(psi)
        signal1 = mt.upsample(2**(level_dec + 1), -2**(level_dec + 1) -2**level_dec)
        
        self.signal = np.abs(signal0) + np.abs(signal1)
                        
        signal_copy = 0
        R_min = (self.HR_min/60.0) * (self.signal.size/self.F       
        self.thr = 0.3*np.max(self.signal)
        
        while True:
            
            start = 0
            end = 0
            Rpeak = []
            
            for x in range(0, self.signal.size):
                
                if self.signal[x] >= self.thr:
                    if start == 0:
                        start = x
                else:
                    if start != 0:
                  
                        if abs(start - x) >= self.Fs * 0.12:
                            end = x - 1
                
                if (start != 0) and (end != 0):
            
                    
                    v_max = max(self.ecg[start:end])
                    v = v_max    
                    
                    t = np.where(self.ecg == v)
                    
                    
                    if t[0].size > 1:
                        idx = t[0]
                        for ii in range(0, t[0].size):
                            if (idx[ii] >= start) and (idx[ii] <= end):
                                t = idx[ii]
                    else:
                        t = t[0]
                    
                    
                    if type(t).__name__ == 'ndarray':
                        t = np.copy(t[0])                    
                    
                    
                    l = len(Rpeak)
                    attr = 0
                    if l > 1:
                        tp, vp = Rpeak[l - 1]
                        if (t - tp) < self.Fs * 0.29:   
                            
                            if v > vp:
                                Rpeak.pop()
                                attr = 1
                        else:
                            attr = 1
                    else:
                        if t > self.Fs * 0.20:    
                            attr = 1
                    
                    if attr == 1:
                        Rpeak.append((t, v))
                        
                    start = 0
                    end = 0
            
            if (len(Rpeak) > R_min):
                break;
            else:
                if (signal_copy == 0):
                    signal_copy = np.copy(self.signal)
                
                lw = int(self.Fs * 0.12)
                wz = np.zeros((1, lw))
                
                l = np.where(self.signal == np.max(self.signal))
                
                signal_copy[abs(int(lw/2) - l[0]):l[0] + int(lw/2)] = wz
                
                self.thr = 0.3*np.max(signal_copy)
           
        self.Rpeak = np.array(Rpeak)
        
        return self.Rpeak
       
    def serie(self):
        self.rr = []
        for i, peak in enumerate(self.Rpeak[:,0]):
            if (i + 1) < self.Rpeak[:,0].size:
                d = self.Rpeak[i + 1, 0] - peak
                self.rr.append(d)
        
        return self.rr
       
  
file multirate.py (link)

Example of detection

#!/usr/bin/python
# -*- coding: UTF-8 -*-
import sys
sys.path.append('/home/.../Scripts/Python')

import numpy as np
import matplotlib
import matplotlib.pylab as pl
from QRSdetect import QRSdetect
import scipy.io
import time

#Configurações do gráficos
matplotlib.rc('xtick', labelsize=8)
matplotlib.rc('ytick', labelsize=8)
   
''' 
ECG
'''
Fs = 500.0
mat = scipy.io.loadmat('/home/.../ECG/C2011/set-a/1002867.mat', mat_dtype=True)
dados = mat['val']

t = np.arange(0, 8.192, 1/Fs)

V1 = dados[6]
V2 = dados[7]
V3 = dados[8]
V4 = dados[9]
V5 = dados[10]
V6 = dados[11]

X_Grupo2 = (np.c_[V1[0:4096], V2[0:4096], V3[0:4096], V4[0:4096], V5[0:4096], V6[0:4096]].T)/100.0

I = dados[0]
II = dados[1]
III = dados[2]
aVR = (-1)*dados[3]
aVF = dados[4]
aVL = dados[5]

X_Grupo1 = (np.c_[I[0:4096], II[0:4096], III[0:4096], aVR[0:4096], aVF[0:4096], aVL[0:4096]].T)/100.0


derivacoes = ['I','II','III','aVR','aVF','aVL','V1','V2','V3','V4','V5','V6']
X = np.c_[X_Grupo1.T, X_Grupo2.T].T

fig1 = pl.figure(1)
fig1.suptitle('Sinais ECG originais', fontsize=14)
fig1.text(0.46,0.03, 'Tempo (s)')
fig1.text(0.07,0.55, 'Amplitude', rotation='vertical')
for i in range(0, 12):
    pl.subplot(4,3,i+1)    
    pl.plot(t, X[i]/np.max(np.abs(X[i])))
    pl.title(derivacoes[i])
    pl.xlim((0, 8.192))
fig0 = pl.figure(2)
fig0.suptitle('Espectrograms', fontsize=14)
fig0.text(0.46,0.03, 'Time (s)')
fig0.text(0.07,0.60, u'Frequency (Hz)', rotation='vertical')
for i in range(0, 12):
    pl.subplot(4,3,i+1)    
    pl.specgram(X[i], Fs = Fs)
    pl.title(derivacoes[i])
    pl.xlim((0, 8))      
       
fig2 = pl.figure(3)
fig2.suptitle(u'R peaks', fontsize=14)
fig2.text(0.46,0.03, 'Sample')
fig2.text(0.07,0.55, 'Amplitude', rotation='vertical')
sigma = 0
media = 0
tempo_medio = 0
for i in range(0, 12):
    pl.subplot(4,3,i+1)

    inicio = time.time()
    
    qrs = QRSdetect(X[i]/np.max(np.abs(X[i])), Fs)
    Rpeak = qrs.detect()   
    fim = time.time()
    serie = qrs.serie()
    print Rpeak[:,0]
    
    tempo_medio += fim - inicio
        
    pl.stem(Rpeak[:,0], Rpeak[:,1])
    pl.title(derivacoes[i])
    pl.xlim((0, 4096))
    pl.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
    med = np.mean(serie)
    media = med + media
    sig = np.std(Rpeak[:,0])/med
    sigma = sig + sigma
    print derivacoes[i]
    print u'Time execution: ', fim - inicio
    print 'Beats: ', Rpeak[:,1].size

figs = [fig1, fig2]

for fig in figs:   
    for ax in fig.get_axes():
        ax.titlesize = 10
        for line in ax.get_lines():
            line.set_color('black')
            
pl.show()



References 

1. Bashir M. E. A., Lee D. G., Aksha M., Yi G. M., Cha E., Bae J., Cho M. C., Ryu K. H, Highlighting the current is- sues with pride suggestions for improving the performance of real time cardiac health monitoring, Information Tech- nology in Bio-and Medical Informatics ITBAM 2010, 6266, 226-233 (2010).
2. Goldberger A. L., Amaral L. A. N., Glass L., Hausdorff J. M., Ivanov P. Ch., Mark R. G., Mietus J. E., Moody G. B., Peng C. K., Stanley H. E., PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals.Circulation101(23): e215-e220 [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215], (2000).
3. Gussak, I., Antzelevitch, C. and Wilde, A. A. M. and Friedman, P. A. and Ackerman, M. J. and Shen, W.- K., Electrical Diseases of the Heart: Genetics, Mecha- nisms, Treatment, Prevention, 600. Springer-Verlag, Lon- don (2008).
4. Kakarala R., Ogunbona P. O., Signal Analysis Using a Multiresolution Form of the Singular Value Decomposition. IEEE Transaction on Image Processing, 10, 724-735 (2001).
5. Kokkinos, P., Physical Activity and Cardiovascular Dis- ease Prevention, 418. Jones and Bartlett Publishers, On- tario (2010).
6. Lee V. S., Cardiovascular MRI: Physical Principles to Practical Protocols, Lippincott Willians & Wilkins, Philadelphia (1966).
7. Lepage R., Boucher J. M., Blanc J. J., Cornilly, J. C., ECG Segmentation and P-wave Feature Extraction: Application to Patients Prone To Atrial Fibrillation, 23rd Annual Inter- national Conference of the IEEE Engineering in Medicine and Biology Society, 1, 298-301 (2001).
8. Mehta S. S., Lingayat N. S., Comparative Study of QRS detection in Single Lead and 12-Lead ECG Based on En- tropy and Combined Entropy Criteria Using Support Vec- tor Machine, Journal of Theoretical and Applied Informa- tion Technology, 8-18 (2007).
9. Rastogi, N. Mehra, R., Analysis of Savitzky-Golay Filter for Baseline Wander Cancellation in ECG using Wavelets, International Journal of Engineering Sciences & Emerging Technologies, 6, 15-23 (2013).
10. Pan J., Tompkins W. J., A Real-time QRS detection al- gorithm, IEEE Transaction on Biomedical Engineering, 32, 230-236 (1985).
11. Randall D. C., El-Wazir Y. M., ECG Interpretation, Hayes Barton Press, Raleigh (2004).
12. Tanaka H., Monahan K. D., Seals D. R., Age-predicted maximal heart rate revisited. Journal of the American Col- lege of Cardiology, 37, 153-156 (2001).
13. Zidelman Z., Amirou A., Adnane M., Belochrani A., QRS detection base wavelet coefficients. Computer Methods and Programs and Biomedicine, 07, 490-496 (2012).

Links to paper: 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - UNESP 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - IX ENAMA 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - Academia.edu 

QRS detection in 12-lead ECG based Multiresolution Singular Value Decomposition - Academia.edu (english version) 

Detecção de Complexos QRS do ECG pela Decomposição em Valores Singulares em Multirresolução - GoogleScholar

No comments:

Post a Comment