This is an open-access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work, first published in JMIR Biomedical Engineering, is properly cited. The complete bibliographic information, a link to the original publication on http://biomedeng.jmir.org/, as well as this copyright and license information must be included.
Background
Heart rate variability (HRV) is derived from the series of R-R intervals extracted from an electrocardiographic (ECG) measurement. Ideally all components of the R-R series are the result of sinoatrial node depolarization. However, the actual R-R series are contaminated by outliers due to heart rhythm disturbances such as ectopic beats, which ought to be detected and corrected appropriately before HRV analysis.
Objective
We have introduced a novel, lightweight, and near real-time method to detect and correct anomalies in the R-R series based on the singular spectrum analysis (SSA). This study aimed to assess the performance of the proposed method in terms of (1) detection performance (sensitivity, specificity, and accuracy); (2) root mean square error (RMSE) between the actual N-N series and the approximated outlier-cleaned R-R series; and (3) how it benchmarks against a competitor in terms of the relative RMSE.
Methods
A lightweight SSA-based change-point detection procedure, improved through the use of a cumulative sum control chart with adaptive thresholds to reduce detection delays, monitored the series of R-R intervals in real time. Upon detection of an anomaly, the corrupted segment was substituted with the respective outlier-cleaned approximation obtained using recurrent SSA forecasting. Next, N-N intervals from a 5-minute ECG segment were extracted from each of the 18 records in the MIT-BIH Normal Sinus Rhythm Database. Then, for each such series, a number (randomly drawn integer between 1 and 6) of simulated ectopic beats were inserted at random positions within the series and results were averaged over 1000 Monte Carlo runs. Accordingly, 18,000 R-R records corresponding to 5-minute ECG segments were used to assess the detection performance whereas another 180,000 (10,000 for each record) were used to assess the error introduced in the correction step. Overall 198,000 R-R series were used in this study.
Results
The proposed SSA-based algorithm reliably detected outliers in the R-R series and achieved an overall sensitivity of 96.6%, specificity of 98.4% and accuracy of 98.4%. Furthermore, it compared favorably in terms of discrepancies of the cleaned R-R series compared with the actual N-N series, outperforming an established correction method on average by almost 30%.
Conclusions
The proposed algorithm, which leverages the power and versatility of the SSA to both automatically detect and correct artifacts in the R-R series, provides an effective and efficient complementary method and a potential alternative to the current manual-editing gold standard. Other important characteristics of the proposed method include the ability to operate in near real-time, the almost entirely model-free nature of the framework which does not require historical training data, and its overall low computational complexity.
Oscillations in the time interval between successive heart beats, referred to as heart rate variability (HRV), have long been known to allow for insight into the intricate control mechanisms of the autonomic nervous system (ANS) [1-4]. Accordingly, research into HRV has attracted considerable attention over the past 4 decades, as evidenced by an exponential growth of published work [5,6].
In a nutshell, heart rhythm and rate can be said to be governed by the 2 competing divisions of the ANS, that is, the sympathetic nervous system and the parasympathetic nervous system. An increased sympathetic input to the sinoatrial (SA) node yields an increase in heart rate (HR) mediated by the release of adrenaline and noradrenaline and the subsequent activation of β-1-adrenoceptors, resulting in an accelerated diastolic depolarization. On the other hand, an increase in the parasympathetic outflow decreases the HR through the release of acetylcholine by the vagus nerve, which binds to and activates M2 muscarinic acetylcholine receptors in the SA node and eventually slows down the diastolic depolarization rate [4,7-9].
Various HRV measures have been established and are usually categorized as either time domain, frequency domain, or nonlinear analysis techniques [2,4]. To a large extent, the popularity of HRV is because of its easy acquisition and seemingly straightforward interpretation.
The Necessity of Preprocessing the R-R Series
It is crucial to realize that by virtue of its very definition, the HRV encompasses only oscillatory phenomena between heart beats resulting from the SA node depolarization [4,10]. While, ideally, all components of an R-R series originate from the SA node, the actual R-R series in both healthy and diseased subjects are contaminated by outliers due to artifacts and heart rhythm disturbances such as ectopic beats (ie, heart beats whose origin is outside of the SA node). Thus arises the necessity to ensure that the HRV analysis is performed on a series representing only the actual normal sinus rhythm (NSR) interbeat intervals, commonly referred to as N-N series (as in normal-to-normal).
The detrimental impact of ectopics on HRV measures is pronounced and well-documented [6,11-17]. In a recent study Stapelberg et al [18] examined the sensitivity of 38 time domain, frequency domain and nonlinear HRV measures to the addition of artifacts in real and artificial 24-hour recordings. In accordance with previous findings, they concluded short-term time domain HRV measures to be more sensitive to the presence of artifacts than their long-term counterparts. Furthermore, frequency domain measures were found to generally be more sensitive than time domain measures, whereas the less commonly used nonlinear measures were shown to exhibit some inherent robustness properties.
Ectopic heart beats are ubiquitous phenomena and not limited to patients with cardiac disorders and diseases. Hingorani et al [19] retrospectively examined the prevalence of cardiac arrhythmias by scrutinizing 24-hour Holter recordings of 1273 healthy volunteers from 22 phase I clinical trials; note that this specific population is not a representative sample of the overall population. The sample was heavily biased in that, consistent with the requirements of phase I clinical trials, subjects were extensively screened by history, physical examination, and laboratory tests and were therefore significantly healthier than the overall population. Crucially, all subjects underwent 12-lead pretrial electrocardiography (ECG) and were excluded if they had any cardiac conduction disorder or disease, including a personal or familiar history thereof. Those exhibiting ≥2 consecutive ectopics as well as bigeminy, trigeminy and quadrigeminy in their 12-lead ECG examination were also excluded. Despite these rigorous exclusion criteria, Hingorani et al [19], among other findings, found premature atrial complexes in 60.8% of the examined Holter records, followed by premature ventricular complexes (PVCs), which were observed in 43.4% of healthy volunteers. While multifocal PVCs occurred in 5.3%, 3.3% exhibited >200 PVCs per 24-hour ECG recordings. A relatively high prevalence of occasional ectopic beats has been reported by numerous other authors as well both in healthy and diseased subjects [20-23].
In summary, the detection and correction of ectopic beats in tachograms are not merely imperative for formal consistency with the HRV definition but arises from the impact of ectopics on HRV measures and the general prevalence of ectopics.
Prior Work and State-of-the-Art
Compared with the overall research interest directed toward HRV applications, the crucial issue of the discrepancy between most real-world R-R series and their respective N-N series has arguably not received as much attention. The actual gold standard of manual R-R series editing, advocated for in the field’s relevant guideline by the European Society of Cardiology and the North American Society of Pacing and Electrophysiology Task Force [4], remains unaltered to this day [11,24].
It should be noted that, aiding and abetting the aforementioned exponential growth of the HRV-related (applied) research, a number of well-designed and widely used software packages for the HRV analysis have been developed and made available to the general public; these include, but are not limited to, Kubios [25,26], Nevrokard’s aHRV [27], and others [28-30]. More often than not, these commercial software packages tend to be closed source and therefore provide only a limited benefit to algorithm developers. On the other hand, it ought to be acknowledged that they are well suited for the needs of practitioners, as they combine advanced preprocessing and analysis features, intuitive graphical interfaces, and support for various file formats.
In automated approaches the task of the detection and correction of outliers is usually tackled as a two-step process.
First, outliers in the R-R tachogram are detected, usually by straightforward thresholding based on a fixed percentage difference to previous intervals or their SD [31], although some more sophisticated methods have also been proposed [12].
Second, the detected (regions of) artifacts from step 1 are either discarded (thereby reducing the effective length of the tachogram), substituted by preceding uncorrupted intervals or by intervals obtained through interpolation, with linear and cubic spline interpolation being the most popular techniques [11].
Objective of This Work
The objective of this study was to introduce a general framework of low computational complexity that allows for an automatic near real-time detection and correction of outliers in the R-R series based on the singular spectrum analysis (SSA). The main novel contributions of this work include (1) the application of a recently proposed model-free lightweight SSA change-point detection (l-SSA-CPD) algorithm [32]; (2) a modification of l-SSA-CPD through the use of an adaptive control limit sequential ranks (AC-SRC) control chart [33] to drastically reduce detection delays; and (3) upon detection of an anomaly, the substitution of the corrupted tachogram segment by an approximation obtained through recurrent SSA forecasting based on a small outlier-free tachogram segment.
An extensive simulation study comprising 198,000 5-minute long R-R series, obtained by randomly inserting varying amounts of simulated ectopic beats in records taken from the MIT-BIH Normal Sinus Rhythm Database (NSRDB), is conducted to assess the artifact detection and correction performance of the proposed algorithm.
Fundamentals of Basic Univariate Singular Spectrum Analysis and Singular Spectrum Analysis Based Change-Point Detection
Recently, this author proposed a low-complexity model-free approach based on SSA and nonparametric cumulative sum (CUSUM) control charts for real-time cardiac anomaly detection, referred to as l-SSA-CPD [32]. It was shown that l-SSA-CPD reliably detects anomalies even when directly applied to unprocessed (ie, no preprocessing was performed) raw ECG and photoplethysmographic records from common databases publicly available through Physiobank [34].
In this study, modifications to the original l-SSA-CPD algorithm will be introduced, and its capabilities will be expanded by adding a recurrent forecasting feature to generate appropriate substitutes of the corrupted tachogram entries. First, however, a brief restatement of SSA and l-SSA-CPD [32] is required.
Basic Singular Spectrum Analysis Algorithm
SSA is a powerful technique of time series analysis owing much of its appeal to an inherent model-free concept combining elements of conventional time series analysis, multivariate geometry and statistics, and signal processing [35]. The univariate basic SSA can be considered as the application of the principal component analysis to the Hankel matrix (obtained through an embedding of the original univariate time series) with the subsequent attempt to reconstruct the original series.
Consider N observations = (x1,..., xN) of a univariate time series and a lag-integer or embedding dimension M (1 < M < < N) also commonly referred to as the window length. The basic SSA algorithm then encompasses the following 4 steps:
1. Embedding = (x1, ..., xN → Χ ∈ ℝM × K)
A trajectory matrix Χ is constructed by mapping into a sequence of K = N − M + 1 lagged column vectors Xj = (xj, ..., xj+M-1)T, j = 1, ..., K of size M, yielding
Χ is called a Hankel matrix because of its characteristic of having equal elements on the antidiagonals. One can think of Χ as multivariate data with M characteristics and K observations and accordingly Xj of Χ as column vectors in the M-dimensional space ℝM.
2. Singular Value Decomposition of Χ
Taking the singular value decomposition (SVD) of Χ decomposes the trajectory matrix into its orthogonal bases and yields a set of M eigenvalues and eigenvectors. Let the eigenvalues of ΧΧT be λ1 ≥ ··· ≥ λM ≥ 0 and U1,..., UM be the respective eigenvectors. Then, with d = max (i, such that λi>0) denoting the rank of Χ, the SVD of Χ can be expressed as the sum of d elementary matrices
with rank 1 matrices and being the eigenvectors of ΧTΧ. Accordingly, the decomposition in Equation (2) is completely characterized by the set of so-called “eigentriples” .
Note that owing to the symmetry of left and right singular vectors, the SVDs of Χ with lag M and K = N − M + 1 are equivalent, entailing the lack of any additional benefit from the use of larger embedding dimensions M > N / 2 [32,35-39].
3. Eigentriple Grouping
Consider the task of extracting a particular component or signal of interest from observed data contaminated by various artifacts such as noise. To do so, during the third stage of the basic SSA, disjoint subsets of indices {1, ..., d} are determined such that the respective systems of eigenvectors span the subspaces associated with those signal components.
In the example of aiming at the separation of a signal from unwanted noise and other disturbances, this entails determining an appropriate subset of indices I = {i1, ..., il}, l < d ≤ M that span an l-dimensional subspace in ℝM, denoted as ℒI ⊂ ℝM = span {UI} = span {Ui1, ..., Uil}, representing the signal whereas the remaining eigentriples with Ī = {i1, ..., id} / I are said to span the noise subspace ℒĪ ⊂ ℝM = span {UĪ}.
The trajectory matrix component ΧI corresponding to the subset I of eigentriples associated with the signal of interest is then
and the component ΧĪ corresponding to the subset Ī = {i1,..., id} I associated with the remainder of the observed signal is
such that the overall SVD of Χ can be rewritten as
In the case of separability (see, eg, [36] at 17), the contribution of ΧI to the entire observed signal Χ is represented by the respective share of eigenvalues .
4. Diagonal Averaging
For perfectly separable signal components, all matrices in the expansion of Equation (5) are Hankel matrices, which require that they all have equal entries on their antidiagonals. In real-world problems, however, perfect separability is rarely achievable. While an approximate separability usually suffices, the final fourth step in the SSA algorithm is required, as the matrices Χi from step 3 would have unequal entries on their antidiagonals. Thus, the Hankelization of all Χi is performed, that is, the compliance with the Hankel structure is enforced by taking the average of each antidiagonal and then replacing each element of that antidiagonal with the determined average. This yields
with being the diagonally averaged matrices.
One can then, for example, easily reconstruct the approximation of the signal of interest through the eigentriples with indices I through the one-to-one correspondence between and the respective time series which provides an approximation of either the entire time series or some components of it, depending on the particular choice of indices I.
Singular Spectrum Analysis Based Change-Point Detection Rationale
The sequential application of SSA for change detection is based on Moskvina and Zhigljavsky [40] (see [39] for an exhaustive discussion; see also [35] for an earlier account of the basic idea). The rationale is to slide 2 (possibly intersecting) moving windows over the time series and to apply the first 3 stages of the basic SSA each time. Think of the first window as adaptively generating a low-dimensional base subspace representation (capturing the main structure in the series) and the second window, which contains at least M new data samples, as generating a test subspace. As both windows are slid over the series, some sort of distance between the base and the test subspace is monitored. If no change-point occurs, the distance remains small, whereas it spikes in the presence of a sudden significant change in the main structure of the series.
The algorithm by Moskvina and Zhigljavsky has been successfully applied to various change detection problems and some improved variations of it have been presented as well (see [32] for a detailed background). Computational complexity, however, ought to be acknowledged as a potentially significant limitation of the said algorithm. Recently, l-SSA-CPD [32] has been proposed as a lightweight alternative. Contrary to the algorithm by Moskvina and Zhigljavsky, which recalculates the SVD of 2 sliding windows each time new observations become available, l-SSA-CPD relies on a nominal low-dimensional signal subspace representation computed only once at startup. As new observations become available, l-SSA-CPD simply allocates them into a vector and calculates a test statistic based on the squared Euclidean distance and the angle between this new data vector and the previously determined nominal reference subspace.
Singular Spectrum Analysis Parameter Selection
The behavior of SSA is largely determined by 2 tuning parameters, namely the window length M and the number l of eigenvalues to be used in the embedding and grouping steps, respectively. While it is undisputed that the successful application of SSA hinges on an appropriate choice of M and l, practitioners are faced with the nuisance that no clear-cut rules as to how they ought to be determined exist.
For a series of length N and window length (or embedding dimension) 1 < M < N, as previously discussed, there is no additional benefit from using M > N / 2 (see, eg, [35] at 69 and[36] at 47). It is widely acknowledged that as the appropriate window length M depends on both the underlying data and the particular application, no universal rule for determining it exists [35-39]. General recommendations vary between M ≃ N / 4 [37] and M ≃ N / 2 [35,36]. In the case of a periodic signal with period T, it is important for M to be, at least, equal to T for the SSA to capture the main structure of the series. Furthermore, M should be proportional to T [35,36,38]. Of note, small windows act like a smoothing filter of width 2M − 1 [36].
Similarly, it is well known that selecting a proper l is also crucial, for issues are known to arise when eigentriple grouping comprises either an insufficient or an excessive number of eigenvalues. More specifically, if l is too small, the SSA is unable to capture the entire signal whereas if it is too large noise components are unwittingly captured. Since the contribution of a particular series component (associated with the corresponding eigentriple) to the entire series is proportional to the respective eigenvalue, it is common practice to select l by identifying the leading eigenvalues through visual inspection of the eigenvalue spectra or by setting a predefined eigenvalue share (see, eg, [35,36,38] for a more exhaustive discussion).
In addition to the above-referred monographs, interested readers are also referred to some fairly recent contributions by Hassani et al [41] and Khan and Poskitt [42,43], which approach the issue through the concept of separability and information theory, respectively.
In the light of the pragmatic suboptimal approach pursued here, a detailed discussion of the SSA parameter selection is, however, beyond the scope of this study and therefore omitted.
MethodsThe Proposed Outlier Detection and Correction Framework
The proposed framework consisted of 2 separate stages, namely the detection and the correction of outlier corrupted segments of the tachogram.
Outlier Detection Using Lightweight Singular Spectrum Analysis Change-Point Detection and Adaptive Control Limit Sequential Ranks Control Charts
Let M, N, l, p, q be fixed integers such that l < M < N / 2 and 0 ≤ p < q. Then l-SSA-CPD proceeds as follows:
1. Initialization at n = 0
Akin to MZ, the first 3 steps of the basic SSA are applied on the interval [n + 1, n + N] to obtain a low-dimensional nominal subspace ℒI = ℒI(n=0) which captures the main structure of the series. This involves embedding as in Equation (1) to obtain a trajectory matrix, which, hereafter, we refer to as the base matrix ΧB = ΧB(0) = ΧB(n=0), taking the SVD of ΧB and obtaining ℒI through an appropriate choice of I = {i1,...,il}, l < d ≤ M with d = max (i, such that λi > 0).
2. Then, for each n = 0,1,... l-SSA-CPD proceeds as follows:
Construction of a test matrix ΧT(n) on the interval [n + p + 1, n + q + M − 1] as
Computation of the detection statistic Dn,I,p,q as
with
and with the angle ∠ (ΧT(n), ℒI) taking values in [0, π / 2] and accordingly . Χj(n) = [xn+j, ..., xn+j+M −1]T are the column vectors of ΧT(n) whereas UI = [Ui1, ..., Uil] denotes the M × l matrix of eigenvectors spanning ℒI and º denotes the Hadamard (or element-wise) product. Note that throughout this study, p = N, q = N + 1 is used, that is, the width Q = q × p of the test matrices ΧT(n) used always equals 1. In other words, each test matrix is actually a single column vector containing M new observations. Q = 1 is primarily chosen to allow for an agile response time and to minimize the computational burden [32].
Monitoring of Dn,I,p,q using the AC-SRC Control Chart
To reduce detection delays, instead of the sequential ranks CUSUM (SRC) by McDonald [44], as used in the conventional l-SSA-CPD [32], a sequential ranks-based CUSUM with adaptive control limits referred to as AC-SRC [33] is used. AC-SRC is inspired by a distribution-free CUSUM proposed by Chatterjee and Qiu [45] where instead of a fixed threshold a sequence of control limits, determined by the bootstrap estimate of the conditional distribution of the CUSUM test statistic given the last time it was zero, is used. AC-SRC carries the approach of Chatterjee and Qiu over to McDonald’s SRC; the main advantage of the SRC (which is that it does not require training data and control limits can thus easily be determined in advance, eg, through Monte Carlo simulations) is retained while detection delays are significantly reduced because of the use of an adaptive sequence of control limits. Furthermore, owing to a trade-off favoring the lower computational burden and simplicity over the optimal performance, AC-SRC lacks some of the refined fine-tuning routines used in the Chatterjee and Qiu method. The author refrains from a detailed discussion of AC-SRC and refers the interested reader to previous studies [33,45].
Let the sequential rank of Dn,I,p,q be denoted as
The AC-SRC then starts with the same recursive equation as the SRC, ie,
with C0 = 0 and k a reference constant. Let Tn denote the so-called “sprint length” expressing the time elapsed since Cn was last 0, ie,
Furthermore, let Yj be a random variable with distribution
Chatterjee and Qiu showed that the conditional distributions of Cn |Tn in Equation (11) are easier to handle than the unconditional distribution of Cn and under some regularity conditions (see [45] and references therein for details), depend only on j and the underlying process generating distribution, but not on n. Then, for any positive integer jmax ≤ n, the distribution of Cn can be approximated by means of the conditional distributions in Equation (11) as
with Y * ∼ Cn |Tn > jmax and I being the indicator function.
It can be shown that [44], given the observed process is in-control, the quantities Rn / (n + 1) in Equation (10) are independent and discrete uniform on . Herein lies the key advantage of AC-SRC, in that for some fixed jmax and k one can determine a sequence of control limits from
with representing a simulation of Equation (10) for up to some NAC-SRC wherein the quantities Rn / (n + 1) are substituted by the set of random variables as discrete uniform on . The sequence of control limits can then easily and without the need for training data samples be determined by means of Monte Carlo simulations as outlined elsewhere [33].
The AC-SRC signals a change to have occurred if
or if
To limit the computational burden, it is recommended to calculate control limits only up to a reasonably small jmax after which the control limit is kept fixed at hAC-SRCjmax. In addition, following recommendations by Chatterjee and Qiu, the reference constant k is set proportionate to the expected average sprint length E {Tn}. Throughout this work, k was empirically calibrated such that , whereas the average run length (ARL) was set to a nominal ARL0=500. A detailed discussion as to the empirical calibration of control charts would, however, be beyond the scope of this work [33]. Note that as k increases, accordingly jmax decreases E {Tn}, that is, the likelihood of Cn bouncing back to 0 increases [33,45]. For the particular use case examined here, the use of a small jmax to allow for agile change detection is recommended (see Tables 1-4 and the respective discussion in the next section).
Finally, as with the SRC in the conventional l-SSA-CPD algorithm, the control chart is reinitialized after a change has been detected (by setting Cn=0 and Tn=0) to allow for the detection of multiple and potentially nearby change-points [32].
Monitoring of the l-SSA-CPD test statistic Dn,I,p,q with N=20, M=10 by means of AC-SRC with jmax=8, E {Tn}=6 is showcased on a 5-minute tachogram excerpt of record 16,273 MIT-BIH NSRDB artificially contaminated with 2 and 8 ectopics in Figures 1 and 2, respectively. Note the remarkably low detection delay due to the use of AC-SRC.
Outlier Correction Using Recurrent Singular Spectrum Analysis Forecasting
Once an outlier has been detected at say n=τ, accounting for the inherent detection delay and accordingly an uncertainty as to the exact outlier location, the interval [τ − Δ1, τ], is designated as corrupted tachogram segment. Note that Δ1 depends on the (average) detection delay and should be chosen carefully. However, it need not be exact or otherwise optimal. Then, a second interval of length Δ2 immediately preceding the corrupted one is used to obtain an anomaly-free forecast to substitute the corrupted segment with.
The rationale of the proposed framework is depicted in Figure 3, with green and yellow rectangles representing the designated uncorrupted and corrupted segments and the detected change highlighted in red.
As the proposed method operates sequentially as new observations arrive, an appropriate uncontaminated data segment (either uncorrupted to begin with or previously cleaned) to serve as a basis for forecasting is always available. While it may be beneficial to use large values for Δ2 (assuming sufficient past observations were available), keeping Δ1, Δ2 as small as possible is recommended; as doing so allows for near real-time operation and significantly reduces the computational burden. For the purpose of this study, simplicity was prioritized over (more often than not unnecessary) optimality.
Detection performance of adaptive control limit sequential ranks with jmax=6 and lightweight singular spectrum analysis change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples.
Record
Sensitivity
Specificity
Accuracy
16,265
0.9712
0.9827
0.9825
16,272
0.9415
0.9802
0.9794
16,273
0.9905
0.9910
0.9909
16,420
0.9771
0.9778
0.9777
16,483
0.9953
0.9846
0.9846
16,539
0.9326
0.9878
0.9870
16,773
0.9692
0.9832
0.9829
16,786
0.9722
0.9758
0.9755
16,795
0.8800
0.9922
0.9905
17,052
0.9521
0.9830
0.9825
17,453
0.9768
0.9904
0.9902
18,177
0.9549
0.9757
0.9754
18,184
0.9719
0.9805
0.9803
19,088
0.9711
0.9917
0.9914
19,090
0.9900
0.9745
0.9745
19,093
0.9744
0.9837
0.9834
19,140
0.9689
0.9920
0.9917
19,830
0.9948
0.9909
0.9909
Detection performance of adaptive control limit sequential ranks with jmax=8 and lightweight singular spectrum analysis change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples.
Record
Sensitivity
Specificity
Accuracy
16,265
0.9335
0.9860
0.9854
16,272
0.8595
0.9850
0.9828
16,273
0.9199
0.9922
0.9913
16,420
0.9141
0.9810
0.9803
16,483
0.9569
0.9893
0.9888
16,539
0.7494
0.9889
0.9863
16,773
0.8781
0.9858
0.9844
16,786
0.9038
0.9779
0.9767
16,795
0.7210
0.9936
0.9897
17,052
0.8259
0.9851
0.9829
17,453
0.9013
0.9933
0.9921
18,177
0.8101
0.9785
0.9771
18,184
0.8508
0.9833
0.9819
19,088
0.8667
0.9924
0.9910
19,090
0.9403
0.9793
0.9789
19,093
0.8920
0.9869
0.9854
19,140
0.8796
0.9946
0.9933
19,830
0.9510
0.9945
0.9941
Detection performance of adaptive control limit sequential ranks with jmax=12 and lightweight singular spectrum analysis change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples.
Record
Sensitivity
Specificity
Accuracy
16,265
0.8816
0.9875
0.9864
16,272
0.8240
0.9871
0.9844
16,273
0.8710
0.9940
0.9925
16,420
0.8741
0.9850
0.9839
16,483
0.8990
0.9900
0.9889
16,539
0.7218
0.9908
0.9879
16,773
0.8202
0.9871
0.9850
16,786
0.8421
0.9814
0.9794
16,795
0.6619
0.9943
0.9898
17,052
0.7753
0.9893
0.9865
17,453
0.8450
0.9939
0.9921
18,177
0.7539
0.9824
0.9805
18,184
0.7952
0.9853
0.9832
19,088
0.8030
0.9930
0.9909
19,090
0.8869
0.9824
0.9814
19,093
0.8379
0.9892
0.9869
19,140
0.8004
0.9951
0.9931
19,830
0.9063
0.9948
0.9939
Detection performance of adaptive control limit sequential ranks with jmax=16 and lightweight singular spectrum analysis change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples.
Record
Sensitivity
Specificity
Accuracy
16,265
0.8623
0.9878
0.9866
16,272
0.7947
0.9885
0.9854
16,273
0.8774
0.9942
0.9929
16,420
0.8585
0.9851
0.9838
16,483
0.8950
0.9904
0.9893
16,539
0.7260
0.9920
0.9890
16,773
0.7951
0.9880
0.9857
16,786
0.8022
0.9814
0.9790
16,795
0.6913
0.9941
0.9900
17,052
0.7866
0.9905
0.9881
17,453
0.8176
0.9951
0.9932
18,177
0.7644
0.9841
0.9823
18,184
0.7846
0.9866
0.9847
19,088
0.8119
0.9936
0.9917
19,090
0.8733
0.9837
0.9827
19,093
0.8299
0.9898
0.9875
19,140
0.7949
0.9953
0.9934
19,830
0.8979
0.9952
0.9944
Adaptive control limits sequential ranks cumulative sum (AC-SRC) monitoring with jmax=8, E {Tn}=6, ARL0=500 of Dn,I,p,q with N=20, M=10, Q=1, 75% of leading eigentriples (top) for a 5-minute tachogram excerpt of record 16273 MIT-BIH Normal Sinus Rhythm Database (NSRDB) with 2 artificially inserted ectopics (bottom).
Adaptive control limits sequential ranks cumulative sum (AC-SRC) monitoring with jmax=8, E {Tn}=6, ARL0=500 of Dn,I,p,q with N=20, M=10, Q=1, 75% of leading eigentriples (top) for a 5-minute tachogram excerpt of record 16273 MIT-BIH Normal Sinus Rhythm Database (NSRDB) with 8 artificially inserted ectopics (bottom).
Schematic depiction of the proposed R-R series outlier detection and correction framework. Adaptive control limits sequential ranks cumulative sum with low detection delay is used to monitor lightweight singular spectrum analysis (SSA) change-point detection’s test statistic; when it signals (red vertical bar), a segment of length Δ1 immediately preceding the detection is labeled as being corrupted (yellow rectangle) and a larger, yet still reasonably small segment of length Δ2 (green rectangle) is used to construct a trajectory matrix for the recurrent SSA forecasting which eventually yields a short-term SSA forecast of length Δ1 which the corrupted segment is then substituted with.
Time series forecasting plays a crucial role in many scientific fields and arguably also represents one of the most popular (real-world) application areas of the SSA. Note that there exist different approaches for SSA forecasting [46], most notably recurrent and vector forecasting. Refraining from any detailed discussion (see [35-38] for excellent monographs on the subject matter), the recurrent SSA forecasting approach used in this study will briefly be outlined. As with other more advanced capabilities that go beyond the realm of basic SSA, forecasting and missing value imputation (note that the 2 are closely related with the latter being the more general problem incorporating forecasting as a special case) in a way sacrifice some of the model-free beauty and appeal of SSA by invoking certain model assumptions. On the other hand, this appears to be outweighed by the successful application of SSA forecasting in numerous areas [47-52]. Furthermore, Golyandina and Zhigljavsky [36] point out that for short-term forecasts, there is actually only very little use of the imposed model, that is, the series approximately satisfy some linear recurrent formulas.
Consider again a univariate time series = (x1, ..., xN) of length N, a fixed embedding dimension M (1 < M < < N) and the respective trajectory matrix Χ = [Χ1, ..., ΧK with K = N − M + 1 as in Equation (1). Let ℒr ⊂ ℝM be a linear space of dimension r < M and U1, ..., Ur (ie, the eigenvectors obtained through the SVD of Χ) be an orthonormal basis in ℒr. Furthermore, let , that is, is the orthogonal projection of the column vector Χi of Χ onto ℒr, and its Hankelized counterpart be the trajectory matrix corresponding to the time series . In addition, given a vector Y = (y1, ..., yM)T ∈ ℝM, let Y∇ ∈ ℝM −1 denote the vector consisting of the first M − 1 components of Y and similarly Y Δ ∈ ℝM −1 the vector consisting of the last M −1 components of Y. Set v2 = π12 + ··· + πr2 where πi|i=1, ..., r denotes the last component of the eigenvector Ui. Assuming ℯM = (0, 0, ..., 0, 1)T ∉ ℒr (ie, ℒr not to be a vertical space), v2 represents the squared cosine of the angle between ℯM and ℒr and v2 < 1 holds. The last component yM of any vector Y = (y1, ..., yM)T ∈ ℒr can then be shown to be a linear combination of its first M − 1 components y1, ..., yM −1 (ie, of Y∇)
with vector A = (a1, ..., αM −1)T being
Note that the representation in Equation (15) does not depend on the choice of the basis U1, ..., Ur in ℒr.
With the required notation now having been introduced, the time series can be defined as
with being the reconstructed series as in basic SSA and yN +1, ..., yN+NF the recurrent SSA forecast of length NF.
Two examples of the automatic detection and correction of outliers using the method proposed in this study are shown in Figures 4 and 5 and compared with the benchmark approach of replacing the corrupted segment with an immediately preceding block of R-R data, referred to as block replacement [28]. Note how in both cases, while all outliers have been detected and corrected, there are ectopic-free segments of tachogram which have also been corrected because of false alarms.
Ectopic beats detected on a 5-minute tachogram excerpt of record 19093 MIT-BIH Normal Sinus Rhythm Database (NSRDB) with 6 artificially inserted ectopics using lightweight SSA change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples and adaptive control limits sequential ranks cumulative sum with jmax=6, E {Tn}=4, ARL0=500 are corrected by means of the proposed method with Δ1=M, Δ2=N (top) or by replacing the corrupted segment of length Δ1 with the immediately preceding block of equal length (bottom).
Ectopic beats detected on a 5-minute tachogram excerpt of record 16265 MIT-BIH Normal Sinus Rhythm Database (NSRDB) with 5 artificially inserted ectopics using lightweight SSA change-point detection with N=20, M=10, Q=1, 75% of leading eigentriples and adaptive control limits sequential ranks cumulative sum with jmax=6, E{Tn}=4, ARL0=500 are corrected by means of the proposed method with Δ1=M, Δ2=N (top) or by replacing the corrupted segment of length Δ1 with the immediately preceding block of equal length (bottom).
Performance Evaluation
The performance and utility of the proposed method are evaluated using records that are publicly available through Physiobank [34]. More specifically, excerpts of all 18 records in the MIT-BIH NSRDB are used. The European Society of Cardiology/North American Society of Pacing and Electrophysiology Task Force guideline [4] recommends to perform the HRV analysis on ECG recordings of at least 5 minutes, and the said length can be considered as the accepted standard for the short-term HRV analysis. Accordingly, for each of the 18 records in MIT-BIH NSRDB excerpts corresponding to the first 5-minute period of NSR are collected. Note that the records in MIT-BIH NSRDB are cardiologist-annotated, that is, they have manually been checked for errors in beat detection and labeling. Thus, one can circumvent the step of QRS complex-detection and use the R-R intervals as provided (and validated by experts) by Physiobank directly. For all 18 records the first 5 minutes were free of ectopic beats and therefore selected. However, it will be pointed out that 4 records contained a couple of extreme outliers (most likely measurement errors), which were removed manually.
Artificial premature ventricular complex obtained by modifying two consecutive tachogram values.
RR(tau)=RR(tau)-RR(tau)/3;
RR(tau+1)=RR(tau+1)+RR(tau+1)/3;
Specifically, the 4 affected records were records number 16,272, 18,177, 19,088, and 19,140, for which 2, 1, 2, and 1 data points were removed, respectively. Besides being readily discernible because they deviate by some order of magnitude from surrounding interbeat intervals, these outliers are also marked as artifacts in the respective annotation files.
Next, the effect of a PVC on the R-R interval was simulated by decreasing the respective entry to 2/3 of its actual value and by increasing the following entry to 4/3 of its value (Textbox 1). That is, for a PVC at say n = τ, one sets
which is consistent with the actual signature of a single PVC (see, eg, [53]) in an R-R series characterized by a premature beat (and accordingly a significantly smaller R-R interval) followed by a compensatory pause (longer R-R interval) followed by a return to baseline (NSR).
Extensive Monte Carlo simulations were then used to determine the performance characteristics of the proposed method, wherein a large number of artificial R-R series based on the 18 5-minute MIT-BIH NSRDB records were created by first randomly drawing an integer between 1 and 6 and then randomly selecting the position(s) within the series where the respective PVCs would be placed (with the constraint that PVCs be at least 5 samples apart).
Results
With the proposed method consisting of 2 stages (detection and correction), a separate analysis of each stage appears appropriate.
Detection Step
In reporting results pertaining to the detection performance of the first stage, established metrics commonly used in the literature are relied upon the sensitivity (Se), specificity (Sp), and accuracy (Acc), which are defined as
and
with TP, FP, TN, and FN representing the number of true positives, false positives, true negatives, and false negatives, respectively. Intuitively, sensitivity quantifies the ability to correctly detect (actual) anomalies, specificity quantifies the proportion of nonabnormal segments that are correctly identified as such, and accuracy combines both of the aforementioned aspects. More formally, Se expresses the empirical statistical power, 1 – Sp the type I (false positive), and 1 – Se the type II (missed detection/false negative) error rate.
As reported previously [32], for an event occurring at a time instance n = τ, one allows for a certain detection delay Δd and considers a signal from the AC-SRC control chart as true positive if it falls in the interval [τ, τ + Δd]. All results presented in this study were obtained using Δd = Δ1 = M.
Note that the design and calibration of a detection algorithm involves an inherent trade-off between false alarms and missed detections and that furthermore reducing the false alarm rate usually entails an increase in the detection delay. For this specific use case, high sensitivity and short detection delays are arguably of the highest priority; for one wants a high likelihood of detecting all actual outliers (to reliably “clean” the R-R series, thereby transforming it into the respective N-N series) with a manageable (ie, small) detection delay to be able to narrow down the corrupted segments in the series as much as possible. The latter is particularly relevant because the proposed method substitutes the (allegedly) corrupted segment by an approximation of the outlier-free segment obtained using recurrent SSA forecasting; one aims to keep the forecasting horizon small to obtain acceptable results without having to tweak the SSA parameters extensively. This comes at the expense of higher false alarm rates, as is evidenced by the results presented in Tables 1-4.
As in the considered design, jmax is linked to the average sprint length E {Tn}, better agility is generally achieved with rather small jmax in a use case, as the one considered here, were AC-SRC monitors for large deviations [33,45].
While AC-SRC achieved satisfactory performance over a wide range of values jmax (results not reported here), upon examination of the results reported in Tables 1-4, which are based on 18,000 R-R series, jmax=6 was chosen to proceed to the second stage of the proposed method (ie, the correction of corrupted segments of the R-R series).
Correction Step
As mentioned earlier this design choice entails higher false alarm rates, as discernible in the 2 examples shown in Figures 4 and 5, which both exhibit actually uncorrupted segments that have been “corrected” by the proposed algorithm because of false alarms, thereby introducing some distortions. More importantly though, both examples also clearly show that all artificially inserted outliers because of PVCs were detected and appropriately corrected.
Average root mean square errors of artificially contaminated R-R series.
Record
RMSESSAa,b
RMSEblock
RRMSEc
Percent increase or decrease
16,265
0.0225
0.0251
0.8964
10.36
16,272
0.0239
0.0354
0.6751
32,49
16,273
0.0136
0.0221
0.6154
38.46
16,420
0.0246
0.0326
0.7546
24.54
16,483
0.0162
0.0196
0.8265
17.35
16,539
0.0282
0.0342
0.8246
17.54
16,773
0.0299
0.0509
0.5874
41.26
16,786
0.0249
0.0366
0.6803
31.97
16,795
0.0323
0.0450
0.7178
28.22
17,052
0.0213
0.0414
0.5145
48.55
17,453
0.0256
0.0344
0.7442
25.58
18,177
0.0296
0.0357
0.8291
17.09
18,184
0.0298
0.0450
0.6622
33.78
19,088
0.0151
0.0198
0.7626
23.74
19,090
0.0187
0.0245
0.7633
23.67
19,093
0.0350
0.0439
0.7973
20.27
19,140
0.0123
0.0156
0.7885
21.15
19,830
0.0110
0.0125
0.8800
12.00
Average
0.0230
0.0319
0.7210
27.9
aRMSE: root mean square error.
bSSA: singular spectrum analysis.
cRRMSE: relative root mean square error.
Results pertaining to the second stage of the proposed method are reported in Table 5, wherein the root mean square error (RMSE) is used to quantify the discrepancy between the outlier-cleaned series and the actual (uncontaminated) series. Let be the uncontaminated N-N series (ie, the ground truth) and denote our allegedly cleaned R-R series.
Calculation of the RMSE is straightforward,
Furthermore, to better compare the RMSE of the proposed method to its competitor (the commonly used block replacement method is used here for benchmarking purposes), consider also the relative RMSE (RRMSE) given by
Note that if the RRMSE, as given in Equation (18) is <1, it implies that the novel approach presented in this paper outperforms the competing benchmark approach by 100·(1 − RRMSE). It shall be emphasized, however, that the block replacement method as it its applied here (to correct for outliers in the R-R series) uses the first stage of the proposed method as well, its performance should thus be seen in the context of being applied to pretty tight segments of the series that ought to be substituted, which is all because of the AC-SRC l-SSA-CPD algorithm of the first stage of the framework presented in this study.
That being said, the results reported in Table 5, which are based on 180,000 R-R series, show that the proposed method, which substitutes corrupted tachogram segments by an approximation of the respective outlier-free segment obtained by
means of recurrent SSA forecasting, clearly outperforms the competing block replacement method. In fact, it outperforms its competitor for all 18 records (on which the extensive simulations built upon), on average by almost 30%.
DiscussionBackground and Relevance
The HRV has long been known to allow for valuable insights into the intricate control mechanisms of the ANS and has a long track record of exponentially growing research interest and output. A crucial prerequisite for any HRV analysis is the exclusion of all artifacts and abnormalities, that is, a clean N-N series is required. Such a clean N-N series is virtually never available, mostly because of various interferences, beats not detected by the QRS complex-detector, and various highly prevalent heart rhythm disturbances.
Compared with the overall research effort directed toward the HRV, the issue of automatic R-R series (pre)processing has received rather little attention. In fact, manual editing still represents the gold standard. It should be noted that every automatic approach comes with risks and benefits; the former lie in the inherently nonzero type I and type II error rates, whereas the latter include, but are not limited to, considerable savings in both cost and time. Considering that error rates for manual approaches are nonzero as well and that they are most likely an increasing function of time (whereas in automatic approaches they are usually not), the choice of which one or which particular combination to use should always be made on an individual basis, considering all requirements and circumstances of the particular case. Accordingly, while the proposed method could be construed as a potential alternative to manual editing, in the light of the above considerations, this author would rather recommend it as a complementary method and or a preprocessing tool.
The main contribution of this work is the introduction of a general framework of low computational complexity that allows for an automatic near real-time detection and correction of outliers in R-R series based on the SSA. While related work pertaining to the use of the SSA in ECG and R-R data processing exists [54], to the best of the author’s knowledge, this study is the first to propose a general SSA-based framework that handles both detection and correction of (unwanted) anomalies in the R-R series.
Principal Findings
An extensive simulation study comprising 198,000 5-minute R-R series, obtained by randomly inserting varying amounts of simulated ectopic beats in records taken from the MIT-BIH NSRDB was conducted to assess the artifact detection and correction performance of the proposed algorithm.
It was shown that the proposed algorithm reliably detects outliers in the R-R series and achieved an overall sensitivity of 96.6%, specificity of 98.4%, and accuracy of 98.4%. Furthermore, it compares favorably in terms of discrepancies of the cleaned R-R series compared with the actual N-N series, outperforming a block replacement approach on average by almost 30%. It should also be emphasized that a suboptimal pragmatic approach was pursued, deliberately refraining from an optimization of the various (SSA and AC-SRC) tuning parameters, which would most likely have resulted in further performance improvements at the expense of simplicity and computational complexity.
Other important characteristics of the proposed method include the ability to operate in near real-time, the almost entirely model-free nature of the framework, and the low computational complexity. Moreover, it should be pointed out that the proposed method is not limited to the R-R series but can be applied broadly, as all of its components are deliberately kept as general as possible. This entails the additional benefit of not being limited to the removal of PVCs (although only this has extensively been investigated and tested so far) but rather being able to deal with all kinds of anomalies that may contaminate the observed series of interest.
Limitations and Future Research
This work has several limitations, some of which open (new) avenues for future research.
A major limitation is because of the lack of established and generally recognized and validated procedures other than manual editing, which transforms the performance comparison among competing methods in a nontrivial task. Second, to increase the readability and decrease the length, the paper fails to provide general recommendations and guidelines as to how various tuning parameters should be chosen. Finally, the SSA represents a very active research field, as in the last couple of years, and several promising advances and new developments have been reported, but thus far have not been considered in the research reported in this work. This, however, does not infringe on the relevance of this work, as the objective was to introduce a simple, general framework of low computational complexity and to investigate the use of SSA in both the detection and correction step.
Several of the abovementioned recent developments, however, appear to open promising avenues for future research. In particular, the author intends to expand and improve the work presented here by focusing on the issue of the SSA parameter selection and by considering latest developments, especially pertaining to SSA forecasting. For the latter, the interested reader is referred to recent publications by Hassani and Kalantari [55-58].
AbbreviationsAC-SRC
adaptive control limits sequential ranks cumulative sum
ANS
autonomic nervous system
ARL
average run length
CUSUM
cumulative sum
ECG
electrocardiography
HR
heart rate
HRV
heart rate variability
l-SSA-CPD
lightweight SSA change-point detection
NSR
normal sinus rhythm
NSRB
Normal Sinus Rhythm Database
PVC
premature ventricular complex
RMSE
root mean square error
RRMSE
relative root mean square error
SA
sinoatrial
SRC
sequential ranks cumulative sum
SSA
singular spectrum analysis
SVD
singular value decomposition
ML was supported by the “Excellence Initiative” of the German Federal and State Governments and the Graduate School of Excellence Computational Engineering at Technische Universität Darmstadt. The views expressed in this article are solely those of the author in his private capacity and do not necessarily reflect the views of Technische Universität Darmstadt or any other organization.
None declared.
BillmanGEHeart rate variability - a historical perspectiveFront Physiol201128610.3389/fphys.2011.0008622144961PMC3225923BillmanGEHuikuriHVSachaJTrimmelKAn introduction to heart rate variability: methodological considerations and clinical applicationsFront Physiol201565510.3389/fphys.2015.0005525762937PMC4340167RajendraAUPaulJKKannathalNLimCMSuriJSHeart rate variability: a reviewMed Biol Eng Comput200612441210315110.1007/s11517-006-0119-017111118ESC/NASPEHeart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and ElectrophysiologyCirculation199603019351043658598068HuikuriHVMäkikallioTAiraksinenKEMitraniRCastellanosAMyerburgRJMeasurement of heart rate variability: a clinical tool or a research toy?J Am Coll Cardiol19991234718788310588197S0735-1097(99)00468-4SassiRCeruttiSLombardiFMalikMHuikuriHVPengCSchmidtGYamamotoYAdvances in heart rate variability signal analysis: joint position statement by the e-Cardiology ESC Working Group and the European Heart Rhythm Association co-endorsed by the Asia Pacific Heart Rhythm SocietyEuropace20150917913415310.1093/europace/euv01526177817euv015LaskeTGShrivastavMIaizzoPAIaizzoPAThe cardiac conduction systemHandbook of cardiac anatomy, physiology, and devices, 3rd ed2015Cham, SwitzerlandSpringer International Publishing215234IaizzoPAFitzgeraldKIaizzoPAAutonomic nervous systemHandbook of cardiac anatomy, physiology, and devices, 3rd ed2015Cham, SwitzerlandSpringer International Publishing235250KaremakerJMAn introduction into autonomic nervous functionPhysiol Meas201705385R89R11810.1088/1361-6579/aa678228304283BerntsonGGBiggerJTEckbergDLGrossmanPKaufmannPGMalikMNagarajaHNPorgesSWSaulJPStonePHvanDMMWHeart rate variability: origins, methods, and interpretive caveatsPsychophysiology199711346623489401419PeltolaMARole of editing of R-R intervals in the analysis of heart rate variabilityFront Physiol2012314810.3389/fphys.2012.0014822654764PMC3358711NabilDBereksi ReguigFEctopic beats detection and correction methods: A reviewBiomedical Signal Processing and Control2015041822824410.1016/j.bspc.2015.01.008SaloMAHuikuriHVSeppänenTEctopic beats in heart rate variability analysis: effects of editing on time and frequency domain measuresAnn Noninvasive Electrocardiol2001016151711174857KimKKLimYGKimJSParkKSEffect of missing RR-interval data on heart rate variability analysis in the time domainPhysiol Meas200712281214859410.1088/0967-3334/28/12/00318057513S0967-3334(07)50608-8KimKKKimJSLimYGParkKSThe effect of missing RR-interval data on heart rate variability analysis in the frequency domainPhysiol Meas200910301010395010.1088/0967-3334/30/10/00519713596S0967-3334(09)19239-0KimKKBaekHJLimYGParkKSEffect of missing RR-interval data on nonlinear heart rate variability analysisComput Methods Programs Biomed2012061063210810.1016/j.cmpb.2010.11.01121194782S0169-2607(10)00287-7LippmanNSteinKMLermanBBComparison of methods for removal of ectopy in measurement of heart rate variabilityAm J Physiol1994072671 Pt 2H411810.1152/ajpheart.1994.267.1.H4117519408StapelbergNJCNeumannDLShumDHKMcConnellHHamilton-CraigIThe sensitivity of 38 heart rate variability measures to the addition of artifact in human and artificial 24-hr cardiac recordingsAnn Noninvasive Electrocardiol20180123110.1111/anec.1248328670841HingoraniPKarnadDRRohekarPKerkarVLokhandwalaYYKothariSArrhythmias Seen in Baseline 24-Hour Holter ECG Recordings in Healthy Normal Volunteers During Phase 1 Clinical TrialsJ Clin Pharmacol2016125678859310.1002/jcph.67926626443ManolioTAFurbergCDRautaharjuPMSiscovickDNewmanABBorhaniNOGardinJMTabatznikBCardiac arrhythmias on 24-h ambulatory electrocardiography in older women and men: the Cardiovascular Health StudyJ Am Coll Cardiol199403152349162581066970735-1097(94)90638-6SimpsonRJCascioWESchreinerPJCrowRSRautaharjuPMHeissGPrevalence of premature ventricular contractions in a population of African American and white men and women: the Atherosclerosis Risk in Communities (ARIC) studyAm Heart J20020314335354011868062S0002870302597489MessineoFCVentricular ectopic activity: prevalence and riskAm J Cardiol19891205642053J56J24807100002-9149(89)91200-9KostisJBMcCroneKMoreyraAEGotzoyannisSAglitzNMNatarajanNKuoPTPremature ventricular complexes in the absence of identifiable heart diseaseCirculation198106636135167226480LabordeSMosleyEThayerJFHeart Rate Variability and Cardiac Vagal Tone in Psychophysiological Research - Recommendations for Experiment Planning, Data Analysis, and Data ReportingFront Psychol2017821310.3389/fpsyg.2017.0021328265249PMC5316555Kubios HRV2018-05-13http://www.kubios.com/TarvainenMPNiskanenJLipponenJARanta-AhoPOKarjalainenPAKubios HRV--heart rate variability analysis softwareComput Methods Programs Biomed201411312102010.1016/j.cmpb.2013.07.02424054542S0169-2607(13)00259-9nevrokard HRV2018-05-13http://www.nevrokard.eu/hrv.htmPichotVRocheFCelleSBarthélémyJChouchouFHRVanalysis: A Free Software for Analyzing Cardiac Autonomic ActivityFront Physiol201611755710.3389/fphys.2016.0055727920726PMC5118625PerakakisPJoffilyMTaylorMGuerraPVilaJKARDIA: a Matlab software for the analysis of cardiac interbeat intervalsComput Methods Programs Biomed20100498183910.1016/j.cmpb.2009.10.00219896748S0169-2607(09)00286-7KaufmannTSütterlinSSchulzSMVögeleCARTiiFACT: a tool for heart rate artifact processing and heart rate variability analysisBehav Res Methods20111243411617010.3758/s13428-011-0107-721573720CliffordGMcSharryPTarassenkoLCharacterizing artefact in the normal human 24-hour RR time series to aid identification and artificial replication of circadian variations in human beat to beat heart rate using a simple thresholdComputers in Cardiology20022912913210.1109/CIC.2002.1166724LangMA Low-Complexity Model-Free Approach for Real-Time Cardiac Anomaly Detection Based on Singular Spectrum Analysis and Nonparametric Control ChartsTechnologies20180215612610.3390/technologies6010026LangMZoubirAMA nonparametric cumulative sum scheme based on sequential ranks and adaptive control limitsProceedings of the 23rd European Signal Processing Conference (EUSIPCO)201523rd European Signal Processing Conference (EUSIPCO)2015 Aug 31Nice, France1984198810.1109/EUSIPCO.2015.7362731GoldbergerALAmaralLAGlassLHausdorffJMIvanovPCMarkRGMietusJEMoodyGBPengCKStanleyHEPhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signalsCirculation2000061310123E2152010851218GolyandinaNNekrutkinVZhigljavskyAAnalysis of Time Series Structure: SSA and Related Techniques2001Boca Raton, FL, USAChapman and Hall/CRCGolyandinaNZhigljavskyASingular Spectrum Analysis for Time Series (SpringerBriefs in Statistics)2013Berlin/Heidelberg, GermanySpringerElsnerJBTsonisAASingular Spectrum Analysis: A New Tool in Time Series Analysis1996New York, NY, USASpringeraSaneiSHassaniHSingular Spectrum Analysis of Biomedical Signals2015Boca Raton, FL, USACRC PressMoskvinaVApplication of the singular spectrum analysis for change-point detection in time seriesPh.D. Thesis, Cardiff University: Cardiff, UK2001MoskvinaVZhigljavskyAAn Algorithm Based on Singular Spectrum Analysis for Change-Point DetectionCommunications in Statistics - Simulation and Computation2003010632231935210.1081/SAC-120017494HassaniHMahmoudvandRZokaeiMSeparability and window length in singular spectrum analysisComptes Rendus Mathematique2011934917-1898799010.1016/j.crma.2011.07.012Atikur Rahman KhanMPoskittDSA Note on Window Length Selection in Singular Spectrum AnalysisAust. N. Z. J. Stat201306195528710810.1111/anzs.12027KhanARPoskittDSSignal Identification in Singular Spectrum AnalysisAust. N. Z. J. Stat20160323581719810.1111/anzs.12143McDonaldDA cusum procedure based on sequential ranksNaval Research Logistics19903764610.1002/1520-6750(199010)37:5<627::AID-NAV3220370504>3.0.CO;2-FChatterjeeSQiuPDistribution-free cumulative sum control charts using bootstrap-based control limitsAnn. Appl. Stat2009033134936910.1214/08-AOAS197GhodsiMHassaniHRahmaniDSilvaESVector and recurrent singular spectrum analysis: which is better at forecasting?Journal of Applied Statistics2017111545101872189910.1080/02664763.2017.1401050HassaniHGhodsiZGuptaRSegnonMForecasting Home Sales in the Four Census Regions and the Aggregate US Economy Using Singular Spectrum AnalysisComput Econ20151217491839710.1007/s10614-015-9548-xSilvaESHassaniHOn the use of singular spectrum analysis for forecasting U.S. trade before, during and after the 2008 recessionInternational Economics201505141344910.1016/j.inteco.2014.11.003HassaniHHeraviSBrownGAyoubkhaniDForecasting before, during, and after recession with singular spectrum analysisJournal of Applied Statistics20131040102290230210.1080/02664763.2013.810193HassaniHHeraviSZhigljavskyAForecasting UK Industrial Production with Multivariate Singular Spectrum AnalysisJ. Forecast2012041032539540810.1002/for.2244MiranianAAbdollahzadeMHassaniHDay-ahead electricity price analysis and forecasting by singular spectrum analysisIET Generation Transmission & Distribution201374346HassaniHThomakosDA review on singular spectrum analysis for economic and financial time seriesStatistics and Its Interface20103337739710.4310/SII.2010.v3.n3.a11CliffordGDTarassenkoLQuantifying errors in spectral estimates of HRV due to beat replacement and resamplingIEEE Trans Biomed Eng200504524630810.1109/TBME.2005.84402815825865ThuraisinghamRAUse of SSA and MCSSA in the Analysis of Cardiac RR Time SeriesJournal of Computational Medicine201320131710.1155/2013/231459KalantariMYarmohammadiMHassaniHSingular Spectrum Analysis Based on L1-NormFluct. Noise Lett2016031501165000910.1142/S0219477516500097KalantariMYarmohammadiMHassaniHSilvaESTime Series Imputation via L1 Norm-Based Singular Spectrum AnalysisFluct. Noise Lett2018061702185001710.1142/S0219477518500177HassaniHKalantariMYarmohammadiMAn improved SSA forecasting result based on a filtered recurrent forecasting algorithmComptes Rendus Mathematique20170935591026103610.1016/j.crma.2017.09.00425904163HassaniHKalantariMA novel signal extraction approach for filtering and forecasting noisy exponential seriesComptes Rendus Mathematique201805356556357010.1016/j.crma.2018.03.00625904163