Abstract
Fragmented QRS (fQRS) has been proven to be an efficient biomarker for several diseases, including remote and acute myocardial infarction, cardiac sarcoidosis, nonischaemic cardiomyopathy, etc. It has also been shown to have higher sensitivity and/or specificity values than the conventional markers (e.g. Qwave, STelevation, etc.) which may even regress or disappear with time. Patients with such diseases have to undergo expensive and sometimes invasive tests for diagnosis. Automated detection of fQRS followed by identification of its various morphologies in addition to the conventional ECG feature (e.g. P, QRS, T amplitude and duration, etc.) extraction will lead to a more reliable diagnosis, therapy and disease prognosis than the stateoftheart approaches and thereby will be of significant clinical importance for both hospitalbased and emerging remote health monitoring environments as well as for implanted ICD devices. An automated algorithm for detection of fQRS from the ECG and identification of its various morphologies is proposed in this work which, to the best of our knowledge, is the first work of its kind. Using our recently proposed time–domain morphology and gradientbased ECG feature extraction algorithm, the QRS complex is extracted and discrete wavelet transform (DWT) with one level of decomposition, using the ‘Haar’ wavelet, is applied on it to detect the presence of fragmentation. Detailed DWT coefficients were observed to hypothesize the postulates of detection of all types of morphologies as reported in the literature. To model and verify the algorithm, PhysioNet's PTB database was used. Forty patients were randomly selected from the database and their ECG were examined by two experienced cardiologists and the results were compared with those obtained from the algorithm. Out of 40 patients, 31 were considered appropriate for comparison by two cardiologists, and it is shown that 334 out of 372 (89.8%) leads from the chosen 31 patients complied favourably with our proposed algorithm. The sensitivity and specificity values obtained for the detection of fQRS were 0.897 and 0.899, respectively. Automation will speed up the detection of fragmentation, reducing the human error involved and will allow it to be implemented for hospitalbased remote monitoring and ICD devices.
1. Introduction
Recently in last 5 years, fragmented QRS (fQRS) has gained clinical significance in the diagnosis of various cardiologic disorders, including remote and acute myocardial infarction, cardiac sarcoidosis, nonSTelevation myocardial infarction, ventricular aneurysm, etc. [1–15]. These studies have shown that fQRS complexes can be an important biomarker for detection of several diseases and has resulted in higher sensitivity and/or specificity than other conventional markers e.g. Qwave, STelevation, etc. [1–15]. However, despite its enormous diagnosis significance cardiologists often disregard or do not report fragmentation in most of the cases except those of bundle branch block (BBB). Therefore, automatic detection and identification of morphologies of fQRS will precisely report on all the cases and may help in finding correlations of potential clinical significance by analysing hundreds of tracings, and hence will facilitate its widespread clinical acceptance, adaptation and application.
In general, the diagnostic information is obtained from standard 12lead ECG using the conventional biomarkers, e.g. P, Q, R, S waves, their duration and peaks, ST elevation/depression, PT interval and (non)inverted T wave [16]. However, there are certain diseases, e.g. cardiac sarcoidosis, myocardial infarction, etc. [2–5,9], which cannot be detected by these conventional biomarkers, e.g. Qwave, ST elevation/depression, PT interval, etc. Patients suffering from these diseases may have to undergo several invasive or noninvasive tests for reliable diagnosis which may be unaffordable [2–5,9]. fQRS in this context has been found to be a marker for aforementioned diseases and its occurrence in standard 12lead ECG makes it an inexpensive and easily available tool for diagnosis. Table 1 provides information on diseases which have been found to have high sensitivity and/or specificity values for diagnosis using fQRS [2–5,8,9] along with the conventional markers/tests that are used to diagnose these diseases and their limitations. Computerized ECG interpretation and feature extraction are being successfully used in a hospitalbased environment [17]. Furthermore, significant prevalence of cardiovascular diseases throughout the world resulted in emerging remote health monitoring systems, which also demand automated, lowcomplexity feature extraction algorithms ported into sophisticated mobile devices or within onboard sensors processing modules ([16,18], http://www.chironproject.eu/). An automated algorithm can find direct applications in remote health monitoring. It has been found that fQRS can be used as a selection criterion for implantation of ICD devices [19,20]. Longterm monitoring is required for this purpose and an automated system will considerably reduce the effort required for verification of the signals. However, to the best of our knowledge fQRS detection and identification of its various morphologies have not been automated and implemented in practice. Here, we propose an automated algorithm for fQRS detection and morphology identification.
2. Material and methods
Figure 1 presents the procedure followed in the proposed method. The raw ECG signal is passed through a preprocessing module comprising baseline wandering removal and denoising, which will be discussed in §2.2. Then, out of existing ECG feature extraction techniques [16,18,21–24] our recently proposed time–domain morphology and gradient (TDMG)based algorithm [16] is applied to extract the QRS complex which is interpolated next and fed as an input to the fragmentation detection and identification (FDMI) module. We skip here the detailed discussion on the TDMG algorithm. Interested readers may consult Mazomenos et al. [16] for the details. If fragmentation is detected, its morphology is then identified. This FDMI module is discussed next. Online QRS detection has been investigated by several researchers [21–24], and hence the availability of an accurately detected QRS can be safely assumed. In this work, DWT has been applied using the ‘Haar’ wavelet. As the Haar wavelet is discontinuous and antisymmetric [25], it is suitable for discontinuity and edge detection [25]. Moreover, owing to its low computational complexity and power consumption, it is suitable for remote health monitoring applications [18]. Hence, Haar is used in this investigation. Detailed discussion on the ‘Haar’ waveletbased transform beyond the scope of this paper. However, interested readers are requested to consult section III of reference [18] for more detail.
2.1. Fragmentation detection and identification module
This section has been further subdivided into two subsections viz. fragmentation detection and morphology identification.
2.1.1. Fragmentation detection
The detailed DWT coefficients behave in a particular manner when a discontinuity^{1} is encountered. A peak or nadir in any signal can be detected by the zero crossing of the wavelet transform. This technique has been previously used to detect QRS complexes of ECG [22]. In this section, we explain this behaviour with the help of figure 2. Consider a line segment with positive slope joining two points and the corresponding bar plot of the detailed DWT coefficient as shown in figure 2a, it can be seen that the coefficient is negative for an increasing segment. Figure 2b shows a line segment with negative slope and the corresponding bar plot of the detailed DWT coefficient which in this case is positive. In figure 2c, a set of points have drawn along with the bar plot of their detailed DWT coefficients. It can be inferred from figure 2a–c that the increasing part of the curve results in negative and the decreasing part of the curve results in positive detailed DWT coefficients, respectively. It should also be noted that the magnitude of the coefficients depends on the slope of the tangent at that point. The greater the slope, the greater the magnitude of the coefficient. Whenever a local extremum appears, there is a transition in the sign of the detailed coefficient depending on the presence of a maxima or minima. We confirm this hypothesis with figure 2d in which both the extrema are identified, i.e. first the maxima and then the minima. This phenomenon has been used in modelling and designing of the proposed fragmentation detection algorithm (FDA). The notches [1–15] that occur in the QRS complex are identified as frequent changes in the sign of detailed coefficients and a peak is identified as a sudden change in sign with a constant follow up of the coefficients with the same sign.
In figure 2, bar plots of detailed DWT coefficients of interpolated QRS complex are plotted and rules for identifying extrema and notches are formulated by observing and correlating the patterns occurring in the QRS complex and detailed coefficients. It is to be noted that the bars may not seem to align because the scales of the xaxis of the wave and the bar plot of detailed discrete wavelet transform (DWT) coefficients are different. When DWT is applied the number of detailed coefficients obtained is half the number of samples on which the DWT was applied. The main purpose of using figure 2 is to show the behaviour of detailed coefficients when an increasing or decreasing part of the wave is encountered. In figure 2c, the wave has seven samples and the total number of detailed coefficients is 4, however, the 4th coefficient is zero. Similarly, in figure 2d, the number of samples is 9 and coefficients are 5, however, the last coefficient is zero, which is indicated in figure 2. Detailed discussion of the DWT can be found in [18]. With the help of few patients, the criteria in table 2 were hypothesized and were then iteratively refined by applying them on leads of 40 subjects such that the mentioned criteria accurately and reliably captured all sorts of discontinuities occurring in the QRS complexes. Table 2 presents the rules framed for detection of local extrema and notches. The algorithm starts evaluating the QRS complex from the leftmost side of the bar plot and proceeds to the right. While traversing through the coefficients, patterns are recognized and if any pattern matches to those mentioned in table 2 the corresponding discontinuity^{2} is realized and noted. Algorithm pointer (k) pointing at a particular point of the QRS complex is incremented according to the presence and type of discontinuity encountered. If a discontinuity is spotted, it is identified and the pointer increments as per the mentioned rule (table 2) otherwise it increments by one. For example, in pattern A1 as shown in table 2 there are two consecutive sign changes in the values of detailed coefficients which can be attributed to the occurrence of a local extremum pair in close proximity, and hence can assumed to be a notch. Similarly, in pattern A2, there are three consecutive sign changes implying that there are three extrema and the one that is identified as a notch depends on the magnitude of the detailed coefficient, the other is identified as an extremum. Other patterns for the identification of notch and extremum can be interpreted in a similar fashion as mentioned in table 2.
2.1.2. Morphology identification
There exist six fundamental morphologies of fQRS and several other variations of RSR’ patterns which were exemplified in the literature [1–15]. Apart from those mentioned, for the sake of completion we have attempted to encompass all the possible variations in RSR’ patterns, for example, some morphologies were found to have a Qwave that was not reported in the six fundamental morphology as shown in B2, C1, etc., in tables 3 and 4. However, all the originally mentioned morphologies identified as fragmented QRS have been included. Tables 3 and 4 summarize and state the criteria for identification of the corresponding morphologies. Table 3 presents all the 10 morphologies that were quantified for QRS ≤ 120 ms, and table 4 presents the criterion for the identification of morphologies with QRS ≥ 120 ms which are generally encountered in practice and appeared in literature. The criterion of identification preferably starts once the number of maxima, minima and notches, the point of occurrence, i.e. positive or negative side of the reference axis, sequence of occurrence and height and depth of R and S waves, respectively, have been obtained from the fragmentation detection step. We have attempted to maintain a clear difference between an extremum pair and a notch. Morphologies with two R waves (R’) viz. A, C, E and I have been considered to dominate over the presence of notches so as to prevent such morphologies to be identified as fQRS. Notched S morphology where the notch is being identified as an extremum pair was encountered, this case has also been taken into account. Similar morphologies for R waves exist and has been identified as RsR’ patterns. These criteria have been formulated based on the appearance of the QRS complex and have been assigned to morphologies such that each one can be identified distinctly, i.e. criteria assigned to morphologies can distinctively identify them without any conflict.
2.2. Preprocessing and feature extraction
Waveletbased techniques have been implemented for lowpower and lowcomplexity applications [18]. Several wavelet transformbased artefact removal algorithms have been proposed in the literature [26–29]. However, there exists no clear demarcation upon the performance of these denoising techniques for comparison. Hence, all these techniques were employed and the denoised signals were visually observed. DWT has been used to remove the baseline wandering and have tested four different variants of wavelet transform to denoise the raw ECG signal. Baseline wandering removal using DWT [26] involved decomposition down to level 9 and the wavelet filter used was Symlet 10. The denoising approaches are as follows:

Approach 1. Denoising and artefact removal using DWT [27]. Decomposition down to level 3. Wavelet filter used was Symlet 4.

Approach 2. Denoising and artefact removal using stationary wavelet transform (SWT) [28]. Decomposition down to level 4. Wavelet filter used was Symlet 4. Empirical Bayes posterior median thresholding was used.

Approach 3. Denoising and artefact removal using undecimated wavelet transform (UWT) [29]. Wavelet filters were Daubechies 6. The level of decomposition was selfdetermined by the code [30]. We used classical standard deviation type of variance estimator and hard thresholding.

Approach 4. Denoising and artefact removal using translation invariant wavelet transform (TIWT) [26]. Wavelet filters used were Symlet 8. Hard thresholding was used and level of decomposition was selfdetermined by the code [31].
It was found that Approach 2 and in very rare cases Approach 1 were tampering with the QRS complex. Figure 3 shows two instances taken from two different subjects, where the tampering with the original signal was found while using Approach 2 (SWT), the denoised signal is plotted against the baseline wanderingremoved signal which was initially verified as not tampering with the signal. It can be seen that notches have been introduced even though the original complex did not have one; other denoising techniques have smoothly retraced the original QRS complex. This visual inspection was done for more than 40 subjects and in rare instances DWT was also found to tamper the complex but none such cases were encountered for Approaches 3 and 4. Approach 3 (UWT) and Approach 4 (TIWT) were satisfactorily denoising the signals. For designing and verification of the algorithm, Approach 4 was adopted as this denoising technique was applied and verified with signals obtained at higher sampling frequencies [32] and satisfactory results were obtained.
The proposed algorithm was implemented on MATLAB (v. 7.10.02010a). Appendix A provides the MATLAB code snippet for the implementation of the baseline wandering and denoising techniques.
3. Experiments and results
This section has been divided into the following subsections. Section 3.1 presents the experimental setup, §3.2 discusses the case studies to understand the working principle of the algorithm, §3.3 presents the evaluation methods used to measure the performance of the algorithm in terms of accuracy, and §3.4 presents the results.
3.1. Experimental setup
The PTB database (PTBDB) [33,34] from PhysioNet has been used for the designing and verification of the proposed algorithm. PTBDB is an unprocessed or raw 15lead database comprising conventional 12 leads and three orthogonal Frank leads digitized simultaneously at a sampling frequency of 1 kHz and captured at the standard speed of 25 mm s^{−1} and 10 mm mV^{−1} with grid intervals being 0.2 s and 0.5 mV. The database was categorized on the basis of cardiac disorders reported and ECGs of patients from various categories were used for designing and modelling the algorithm. PTBDB consists of patients belonging to various diagnostic classes viz. myocardial infarction, cardiomyopathy/heart failure, BBB, dysrthymia, myocardial hypertrophy, Valvular heart disease, myocarditis, healthy controls and other miscellaneous conditions [33,34]. A high sampling rate is desired to capture the occurrence of highfrequency notches in QRS complexes. A simple linear interpolation is then applied on the preprocessed signal to further increase the number of samples so that the detailed coefficients obtained after applying DWT can accurately detect all discontinuities (local extrema and notches), as DWT diminishes the time resolution by a factor of 2. Equation (3.1) shows the methodology adopted for interpolation 3.1
A mean was calculated for every two consecutive samples and was inserted in between them thus increasing the total number of samples. Interpolation in this manner is extremely simple and does not affect the points of extrema. Requirement of interpolation is based on observational interpretations during the designing of the algorithm. Upon interpolation, it was observed that the detailed coefficients were found to capture all the discontinuities without fail. The MATLAB code snippet for the implementation of linear interpolation and DWT can be found in appendix A.
3.2. Case studies
The outputs obtained from the algorithm are the number of extrema points (maxima and minima), notches, time instants at which they occur, magnitude of extrema, the occurrence on the lower or upper half of the magnitude–time plane and the morphology of fragmentation encountered. In this section, we discuss eight different cases and the output obtained from the algorithm for them which are intended to show the working principle of our proposed algorithm. These complexes were randomly selected from the large database avoiding the selections where similar morphologies were selected. Figure 4a–h shows the plot of these cases with part ‘(i)’ of the plot showing the interpolated QRS complex and part ‘(ii)’ of the plot showing the bar plot of the detailed coefficients obtained after applying DWT. Notches and extrema are denoted by circle and rectangle, respectively. The following cases correspond to the QRS complexes shown in figure 4. This case study provides insight into the working of the algorithm and will help in the reproduction of the work. Rules mentioned in table 2 have been used to detect the discontinuities, and criteria mentioned in tables 3 and 4 have been used to identify the morphology.

Case 1. This case is an example of notched R (rsR’) morphology. Number of maxima, minima and notches are 1, 2 and 1, respectively. On examining figure 4a(i), it can be seen that a minima (<0) is encountered first followed by a notch (>0), maxima (>0) and a second minima (<0).

Case 2. This case is an example of RsR’ without ST elevation. Number of maxima, minima and notches are 2, 2 and 0, respectively. In figure 4b(ii), a maxima (>0) is encountered first followed by minima (>0), maxima (>0) and minima (<0).

Case 3. This case is an example of Rsr’. Number of maxima, minima and notches are 1, 2 and 1, respectively. In figure 4c(ii), a minima (<0) is encountered first followed by maxima (>0), notch (>0) and minima (<0).

Case 4. This case is an example of notched S. Number of maxima, minima and notches are 2, 2 and 0, respectively. In figure 4d(ii), a maxima (>0) is encountered first followed by minima (<0), maxima (<0) and minima (<0).

Case 5. This case is an example of rSr’. Number of maxima, minima and notches are 2, 1 and 1, respectively. In figure 4e(ii), a maxima (>0) is encountered first followed by notch (>0), minima (<0) and maxima (<0).

Case 6. This case is an example of fQRS. Number of maxima, minima and notches are 2, 3 and 0, respectively. In figure 4f(ii), a minima (<0) is encountered first followed by maxima (>0), minima (<0), maxima (>0) and minima (<0).

Case 7. This case is an example of fQRS. Number of maxima, minima and notches are 3, 4 and 0, respectively. In figure 4g(ii), a minima (<0) is encountered first followed by maxima (>0), minima (>0), maxima (>0), minima (>0), maxima (>0) and minima (<0).

Case 8. The length of QRS is more than 120 ms (QRS ≥ 120 ms). This case resulted in fragmented BBB (fBBB) with 3 or more R waves (R’). Number of maxima, minima and notches are 3, 2 and 1, respectively. In figure 4h(ii), a maxima (>0) is encountered first followed by minima (<0), maxima (>0), notch (>0), minima (>0), maxima (>0).
3.3. Evaluation criteria
Forty patients were selected at random from the database, and the QRS complexes were extracted using the TDMG feature extraction algorithm and were examined by the cardiologists. Out of 40 patients, nine patients were removed from the study owing to the poor quality of their ECGs as per the suggestion of the cardiologists. The remaining QRS complexes were then independently examined by the two cardiologists in a blind fashion, and finally a consensus was reached to produce the final result which will be hereby referred as cardiologist's status (CS). The measurements obtained from the cardiologists are assumed to be the gold standard. The QRS complexes were then input to the FDA and the results obtained were compared with CS, and sensitivity and specificity values were calculated. The following will discuss the test results and test nature:

True positive. The cardiologists detected fragmentation in a particular lead of the patient and the algorithm reported correctly.

False positive. The cardiologists did not detect fragmentation but the algorithm reported the presence of fragmentation in a particular lead.

True negative. The cardiologists did not detect fragmentation and the algorithm reported correctly.

False negative. The cardiologists detected fragmentation in a particular lead but the algorithm could not detect.
The sensitivity and specificity were calculated using the following equations: 3.2and 3.3
3.4. Experimental results
The QRS complexes obtained from TDMG were evaluated by two experienced cardiologists and were simultaneously used to obtain results from the algorithm. From 372 leads (31 patients) selected, i.e. only 12lead ECG of 31 patients and not the Frank leads, the results of 89.8%, i.e. 334 leads from FDA complied with that of CS. The sensitivity and specificity values obtained were 0.897 and 0.899, respectively. As cardiologists are often interested in fragmentation in BBB subjects, out of randomly selected 31 patients five patients had BBB and the sensitivity and specificity were found to be 0.932 and 0.933, respectively. It is to be noted that in this paper we have evaluated the presence or the absence of fragmentation at a lead level rather than patient level which is usually the case in the literature. It is important that the algorithm detects fragmentation correctly in every lead of the patient as well identifies its various morphologies. The sensitivity and specificity values were obtained for the presence or the absence of fragmentation and not for the morphology obtained by the algorithm and cardiologists. We have not found any work that standardizes the morphologies of fragmentation for the common agreement, and hence we have avoided the evaluation of the morphological results obtained from the algorithm.
4. Discussion
From figure 4a–h, we can see that whenever the gradient of the wave is high the magnitude of the detailed coefficients is more. Sensitivity of the algorithm can be estimated from the fact that it not only captures the formation of an extremum pair but also sudden changes in gradient which did not result in an extremum pair, e.g. in figure 4f,h a sudden gradient change before the occurrence of final maxima (figure 4f(ii)) and after the first maxima (figure 4h(ii)) can be seen (denoted by a fourpoint orange star) with the magnitude of detailed coefficients suddenly decreasing and then increasing. This highlights the sensitivity of the detailed coefficients to sudden changes in the ECG, and hence it is important to denoise the ECG before applying the FDA algorithm otherwise a noisy part in the QRS complex may be detected as a notch, which upon denoising may not be present. Hence, baseline wandering removal and denoising is required for the correct detection of morphology of fQRS.
Sampling rate plays an important role in capturing fragmentation. We have observed that a sampling rate of 2 kHz will be appropriate for algorithm implementation. This can be deduced from the role played by interpolation to double the number of samples as the sampling rate of the PTB database is 1 kHz. The significance of number of samples can be seen from figure 5a–d (notches have been denoted by a red circle). Figure 5a–d(i)(iii) shows a plot of interpolated QRS along with a bar plot of detailed coefficients. Figure 5a–d(ii)(iv) of original QRS along with the bar plot of its detailed coefficients. All four figures show that on increasing the number of samples, undetectable notches can also be identified. This has been our main motivation behind interpolating. A notch that is detectable will only become elongated or may become converted to an extremum pair; this can be resolved with the help of criteria and postulates for morphology realization and identification, respectively, but an undetectable notch will result in an incorrect determination of the number of notches, and hence the morphology, which may lead to an incorrect diagnosis. Figure 5a was identified by the algorithm as notched S, 5B as Rsr’, 5C as fQRS and 5D as rSR’.
When patterns similar to table 2, A1 and A3 are encountered, it unfailingly denotes the presence of a notch but when patterns like A2 and A4 are encountered it becomes difficult to interpret whether such a pattern should be treated as a notch or an extremum pair. Figure 5b(i)(iii) shows a similar case. In the absence of magnitude criterion, the maxima would have been identified as a notch and the peak of the notch would have been identified as a maxima.
Properties of the detailed coefficients as explained in figure 2 can be used to explain and understand the occurrence of such patterns. Similarly, if the magnitude criteria stated in tables 3 and 4, B1 and B2 are spared then on implementation of the algorithm in figure 5c(ii)(iv) would have resulted in detection of the first maxima and second minima as two notches which is inaccurate. On iterative refinement of the algorithm and evaluating it on more than 40 patients, we deemed it necessary to incorporate the magnitude criterion in the postulates. In the designing and modelling phase, the method of QRS detection was adopted from the appendix of [25]. To obtain accurate results from the algorithm, it is necessary to input an accurately detected QRS complex. Any extra discontinuity will lead to the identification of some other morphology. In the worstcase scenarios, the algorithm will identify the morphology as fQRS, leaving no fQRS undetected.
5. Conclusion
We have introduced a novel approach for detection of discontinuities in the QRS complexes and have verified it using PTBDB. A less complex ‘Haar’ wavelet was used for lowpower consumption which can be used in batteryoperated devices viz. mobile phone/PDA/tablets [18] or batteryoperated ICD devices. This approach is not signal specific and the method can be applied to any other kind of biomedical signal for detection of its specific important aspects and features. For ECGspecific applications, we have formulated the postulates for detection of notches and extrema and have proposed criteria for identification of various morphologies. The significance of denoising techniques and all types of discrepancies encountered have been discussed.
Funding statement
This work is partly supported by the DIT, India under the ‘Cyber Physical Systems Innovation Hub’ under grant no.: 13(6)/2010CC&BT, dated 1 March 2011 and ‘IOT for Smarter Healthcare’ under grant no.: 13(7)/2012CC&BT, dated 25 February 2013.
Appendix A
The appendix provides a snippet of the Matlab codes which will be helpful in reproduction of the work. The part of the code considered important and necessary has been provided. Figure 6 provides the Matlab code for implementation of denoising techniques, and figure 7 provides the code for implementation of interpolation and DWT.
Footnotes
 Received August 16, 2013.
 Accepted September 25, 2013.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.