Unraveling cardiac arrhythmia frequency: comparative analysis using time and frequency domain algorithms

During cardiac arrhythmia, the heart frequency is an important physiological parameter that can be identified by analyzing electrocardiogram (ECG) signals. However, the accuracy of the frequency estimation becomes increasingly challenging as the ECG morphology becomes more complex, for example, during transitions from tachycardia to fibrillation. In this paper, the authors compare seven conventional and novel time- and frequency-domain methods for cardiac arrhythmia frequency analysis, including an algorithm used in implantable cardioverter defibrillators. The objective of this study is to identify the approaches that reveal the potential presence of a dominant frequency and its role in characterizing different arrhythmia types. By evaluating the strengths and weaknesses of each method, the authors aim to establish an informative framework for extracting meaningful insights from electrocardiogram data in the context of cardiac arrhythmia frequency. In order to ascertain the statistical relevance of the methods, a dataset comprising 112 ECGs from arrhythmic murine hearts was analyzed. Additionally, a dataset comprising human arrhythmia data was examined to validate the techniques presented. The R-library, which contains the frequency determination algorithms, as well as the murine data set, is made available to the reader for the purposes of further testing and supplementation.

1 Introduction

The electrocardiogram (ECG) has been a fundamental tool in the understanding of cardiac electrophysiology for over a century. However, despite significant technical advances, the mechanisms initiating and sustaining cardiac arrhythmia remain to be elucidated. A hallmark of the ECG that differentiates normal sinus rhythm from life-threatening ventricular tachyarrhythmias is the heart rate (or heart frequency), which inversely correlates with the cardiac oscillation period or cycle length (CL). As a critical clinical marker, heart rate monitoring has become essential for diagnosing and managing cardiac diseases, particularly in identifying pathological rhythms that require immediate therapeutic intervention. In this regard, the development of algorithms to automatically detect heart rate, along with their subsequent integration into implantable cardioverter-defibrillators (ICDs) and automated external defibrillators (AEDs), has significantly improved early detection and treatment of dangerous rhythms, such as ventricular tachycardia (VT) and ventricular fibrillation (VF) (Moss et al. (1996), Moss et al. (2002); Weisfeldt et al. (2010)).

The development of ECG mapping technologies (Cakulev et al., 2021) and optical mapping of intact hearts (George and Efimov, 2019) has significantly contributed to the understanding of ventricular tachyarrhythmia mechanisms. Optical mapping has demonstrated that spiral-like activation waves are a primary driver of ventricular arrhythmia (Gray et al., 1995b), whereas ECG mapping has revealed re-entrant circuits as a potential source of VT (Kastor et al., 1972). The combination of these techniques has enabled the correlation of ECG patterns with arrhythmia mechanisms, providing a more comprehensive understanding of the underlying processes. Specifically, studies have shown that the anchoring of a single spiral wave to a anatomical or functional heterogeneity can give rise to monomorphic ventricular tachycardia (mVT) as demonstrated by Annoni et al. (2017), while the deviation of spirals within the myocardial substrate can lead to the distinctive morphology of polymorphic ventricular tachycardia (pVT), which is characterized by a complex and irregular ECG pattern of activation (Gray et al., 1995a). Furthermore, ventricular fibrillation (VF) is marked by a chaotic and disorganized ECG pattern, originated from multiple, rapidly drifting spiral waves (Jalife and Gray, 1996).

Early research employed the Fast Fourier Transform (FFT) to analyze the frequency components of ECG signals during arrhythmia. This analysis revealed that different arrhythmia morphologies are associated with distinct frequency spectra (Nolle et al., 1980; Kinoshita et al., 1992; Clayton et al., 1993). In particular, multiple peaks in the frequency-domain representation of ECG signals have been identified, which are linked to the presence of multiple spiral waves with varying rotation periods (Pandit and Jalife, 2013). This finding suggests that the frequency components of ECG signals can provide valuable information about the underlying mechanisms of arrhythmia.

The determination of arrhythmia frequency also known as dominant frequency has gained significant importance for therapeutic approaches, as studies have shown that delivering multiple low-energy shocks within a specific time frame relative to the arrhythmia frequency increases the likelihood of successful defibrillation (Li et al., 2009; Ambrosi et al., 2011). One approach that leverages this concept and is already employed in clinical settings is Anti-tachycardia Pacing (ATP). ATP involves delivering low-voltage pulses at a frequency higher than the detected heart rate, a technique known as overdrive pacing. While ATP is often successful in terminating VT, its effectiveness is limited in the case of VF (Ricci et al., 2001; Maria et al., 2017). Conversely, other studies have demonstrated that far-field low-energy overdrive pacing is an effective strategy for the termination of atrial fibrillation (AF) and VF (Fenton et al., 2009; Luther et al., 2011; Ji et al., 2017; Hornung et al., 2017).

In the light of the importance of determining the cardiac frequency of arrhythmia, for both diagnostic and therapeutic purposes, we have developed an open-source platform for interspecies ECG arrhythmia frequency analysis. This platform, implemented in an R-library (Wickham and Bryan, 2023), provides a flexible and user-friendly tool for researchers to analyze ECG data. The R-library FibFreq includes three frequency-domain approaches: (i) the well-established maximum of the power spectrum, herein referred to as the maximum of the periodogram (MaxP); (ii) a method to estimate the spectrum via the correlation function, called Adaptive Variable Period (AVP), which is a faster version of the algorithm presented by Diaz-Maue and Nobach (2022); and (iii) a spectrum-fitting method using a Lorentz distribution (Lorentz), as presented in Diaz-Maue et al. (2022). In addition, time-domain approaches included in the R-Library encompass the following methods: (iv) a sinus fit to the time series (SinFit); (v) distance analysis between succeeding peaks (M2M); (vi) a zero-crossing algorithm (UZC), both of which can be considered as CL estimators; and (vii) the Biotronik heart rate analysis algorithm (ICD) as outlined in Brüggemann et al. (2016). The authors also provide the dataset utilized for the analysis of arrhythmia, allowing for further examination by interested researchers.

2 Experimental data set

2.1 Mouse data

The first dataset utilized in this study was obtained from a series of experiments conducted on intact hearts from transgenic mouse containing the Channel Rhodopsin-2 (ChR2) encoding gene. The light-sensitive ion channel, ChR2, is expressed exclusively in cardiac tissue and it enables the elicitation of action potentials through the use of light at an appropriate wavelength. In comparison to wildtype mice of the same background, this model exhibited normal electrical propagation throughout the cardiac tissue (Zaglia et al., 2015). The induction of arrhythmia was achieved through the delivery of a rapid series of optical stimuli by an external blue light source as demonstrated in Diaz-Maue et al. (2022). For this analysis, a subset of seven adult mice, aged between 91 and 113 weeks, was selected.

Single-Lead epicardial electrograms were recorded using a surface electrode positioned at the center of the left ventricle and a reference electrode placed in close proximity to the wall of the perfusion reservoir as shown in Figure 1. The acquisition system is equipped with a low-pass hardware filter with a corner frequency $f_{\text{LPhw}}=35\text{\hspace{0.17em}Hz}$ and a high-pass filter with a corner frequency $f_{\text{HPhw}}=0.5\text{\hspace{0.17em}Hz}$. The data set provided in this paper encompasses $N_{\text{arr}}=112$ ECG signals, each containing one second of arrhythmia, sampled at a frequency $f_{\text{s}}=1\text{\hspace{0.17em}kHz}$. In order to eliminate the power line hum (50 Hz) the ECG time series underwent digital low-pass filtering with a fifth-order Butterworth filter with a cutoff frequency of $f_{\text{LPsw}}=45\text{\hspace{0.17em}Hz}$ (van Boxtel et al., 2021).

Figure 1

Figure 1. The ECG signal from Langendorff-perfused intact mouse hearts was measured using a surface electrode positioned on the left ventricle of the heart and a reference electrode located near to the reservoir wall. The ECG signal was then amplified and digitized with an acquisition system.

In consideration of the potential loss of information during signal processing, the preprocessing steps were kept to a minimum. Hence, no windowing function was used, however, to ensure the comparability of ECGs, the signal mean was subtracted to eliminate direct current (DC) components. This step is necessary to account for any DC offset that may be present in the ECG signals. Four examples of the examined mice data are shown in Figure 2, these are ordered from very regular (top) to irregular (bottom).

Figure 2

Figure 2. ECG signals of arrhythmic mouse hearts are displayed in increasing complexity from the top to the bottom. Figure created using the ggplot2 library of R (Wickham et al. (2019).

2.2 Human data

Three cases of human ventricular arrhythmia were retrieved from record 420 (VT), 418 (Ventricular flutter (VFlu)) and 422 (VF) of the MIT-BIH-malignant-ventricular-ectopy-database-1.0.0. For detailed information, please refer to Moody and Mark (2001) and Greenwald (1986). The ECG signals were digitized with a 12-bit Analog-to-Digital Converter (ADC) with a voltage resolution of 200 adu/mV (analog-digital-unit (adu)) and a sampling frequency of $f_{\text{s}}=250\text{\hspace{0.17em}Hz}$. The subsequent analysis was conducted using the R-Library on recording ECG:0 (see Figure 3). To prepare the data for further processing, the mean value of the ECG was subtracted. Subsequently, possible power line interferences were filtered using a low-pass fifth-order Butterworth filter with a corner frequency of $f_{\text{LPhum}}=20\text{\hspace{0.17em}Hz}$. The removal of all frequency components with a frequency lower than $f_{\text{HPhum}}=0.5\text{\hspace{0.17em}Hz}$, including baseline wander, was achieved through the application of inverse Fourier filtering. The presence of movement artifacts was not detected in the extracted segments. However, these issues can be effectively mitigated by using Fourier filtering method as previously outlined.

Figure 3

Figure 3. The ECG signals of human origin presented in this study exemplify three instances of ventricular tachyarrhythmia retrieved from the MIT-BIH-malignant-ventricular-ectopy-database-1.0.0 (Moody and Mark, 2001; Greenwald, 1986), each exhibiting unique characteristics. Figure created using the ggplot2 library of R (Wickham et al. (2019).

3 Signal processing methods

We will examine seven techniques for determining arrhythmia frequency, distinguishing between frequency-domain and time-domain. In the frequency-domain, we will discuss the Maximum of the periodogram, the Lorentz fit to the amplitude spectrum and the Adaptive Variable Period. In addition, we will consider time-domain methods, including the Sinus Fit to the time series, the Average Maximum to Maximum Distance, the Average Distance between upward Zero-Crossings and an algorithm implemented in ICDs.

The following notations will be applied: the ECG time series $x_j$, where $j=1,\dots ,N$ with $x_j$ representing the ECG at time $t_j=j\cdot \mathrm{\Delta }$. The sampling interval is denoted by $\mathrm{\Delta }$ and the number of time series elements by $N$. For the mouse data, the number of time series elements is $N=1000$ and the sampling frequency $f_{\text{s}}$ is 1,000 samples per second, corresponding to a sampling interval of $\mathrm{\Delta }=1/f_{\text{s}}=1\text{\hspace{0.17em}ms}$; Each human ECG time series comprises $N=1250$ data points and is sampled at a frequency of $f_{\text{s}}=250\text{\hspace{0.17em}Hz}$ which corresponds to a sampling interval of $\mathrm{\Delta }=4\text{\hspace{0.17em}ms}$ and an observational interval of $5\text{\hspace{0.17em}s}$.

3.1 Frequency-domain methods for the estimation of the arrhythmia frequency

The conversion of a time series into the frequency-domain can provide a more comprehensive insight into the prominent frequencies present. Therefore, the analysis of the frequency content by estimating the power spectrum has become a wide-spread method. However, it is essential to note that there are multiple estimators of the power spectrum, which exhibit varying degrees of estimation error and frequency resolution. In this section, three methods for the estimation of the frequency spectrum and the identification of a dominant frequency are presented.

3.1.1 Maximum of periodogram (MaxP)

The simplest estimator of the power spectrum is based on the discrete Fourier transform (DFT)

$X_k={\displaystyle \sum _{j=1}^N}{x}_j{\mathrm{e}}^{-2\pi \mathrm{i}\frac{jk}{N}}(1)$

of the signal $x_j,j=1,\dots ,N$, which can be computed with fast algorithms, such as the Cooley–Tukey algorithm for FFT (Cooley and Tukey, 1965).

For real valued time series such as the ECG signals considered in this study, the entire power spectrum $P(f_k)$ is symmetric. Provided that the signal has no DC component and also no contributions at half the sampling frequency $f_{\text{s}}=1/\mathrm{\Delta }$, it is sufficient to consider the one-sided spectrum with the frequencies $f_k=k\cdot f_{\text{s}}/N,k=1,\dots ,N/2-1$ only. That leads to the factor two in the numerator of the full-scale periodogram

$P\left(f_k\right)=\frac{2}{N^2}|X_k{|}^2,(2)$

which is calculated using the squared modulus of the DFT. Then the sum over the values of the power spectrum $\sum P(f_k),k=1,\dots ,N/2-1$ equals the variance of the signal ${\sigma }^2$.

If the considered time series can be modeled as a superposition of harmonic functions from these frequencies $f_k$, then the frequency $f_{\text{max}}$ that corresponds to the maximum of the power spectrum is the dominant frequency of the arrhythmia.

In this paper, the fraction of explained variance, denoted by ${\rho }^2$ will be used as a measure for the goodness of the dominant frequency estimation and is defined as:

${\rho }^2=P\left(f_{\text{max}}\right)/{\sigma }_x^2=P\left(f_{\text{max}}\right)/\sum P\left(f_k\right)(3)$

Where $P(f_{\text{max}})$ is the power at the dominant frequency. Values of ${\rho }^2$ close to 1 indicate a good description of the signal, while small values suggest problems with the concept of dominant frequency for the considered signal.

The periodogram approach is known to have limitations, such as spectral leakage and false estimates, particularly in the presence of trends, non-sinusoidal components, or frequencies that are not multiples of the fundamental frequency $f_1$. As illustrated in Figure 4A, the amplitude spectrum of an ECG signal obtained with the FFT may exhibit multiple frequency peaks, some of which are attributable to the leakage effect. Despite these constraints, the periodogram remains a widely used method for dominant frequency estimation, and it serves as a useful point of reference for comparison in this study.

Figure 4

Figure 4. Illustration of seven techniques for estimating the arrhythmia frequency: (A) Spectral Techniques: maximum of the periodogram (blue curve) at $f_{\text{MaxP}}$ (blue dashed vertical line) and maximum of the Lorentz fit to the Fourier amplitudes at $f_{\text{Lorentz}}$ (red dashed vertical line). (B) Adaptive Variable Period (AVP) technique: shown are the ECG (blue curve) and three periods (indicated by vertical lines of different colors). The period that corresponds to the grey dashed vertical lines matches best the oscillation. (C) Time-domain techniques: shown is the ECG (blue curve) and (i) the best fitting sinusoidal model (light blue thick line), (ii) the maxima identified by the sensing algorithm implemented in ICDs (beige squares), (iii) the maxima identified by the M2M algorithm (green filled circles) and (iv) upward zero-crossings (ocher filled triangles) of the filtered ECG. Figure created using the ggplot2 library of R (Wickham et al., 2019).

3.1.2 Adaptive variable period (AVP)

The leakage of the FFT routine is specifically addressed by this estimation method. In instances where the signal duration does not correspond to an exact integer multiple of the fundamental period of the periodic signal investigated, the FFT spectrum will be sampled at frequencies that do not align with the fundamental frequency to be detected.

In order to mitigate this issue, the measured signal can be shortened by a certain number of samples, and the FFT from the shorter sequence can be calculated. In the event that the duration of the shorter sequence becomes a integer multiple of the fundamental period, the peak in the FFT spectrum will stand out from the other values in the spectrum and the peaks in all other spectra with wrong assumptions of the duration, making it easier to identify one dominant frequency.

The AVP method essentially investigates how good an assumed fundamental period of the signal matches the observed data. Unlike the approach described above, the AVP does not require discarding any samples by shortening the sequence.

To this end, first the correlation function $R_k,k=0\dots p-1$ is derived using the entire sequence with $N$ values, where a $p$ samples long period, shorter than $N$, can be assumed in the calculation. In essence, for each $k$ between 0 and $p-1$, all pairs of values $x_i$ and $x_j$ that have the appropriate distance $k$ are considered, under the assumption, the signal would be periodic with a period of $p$. Formally, this can be expressed in the following way, where the indices $i$ and $j$ comply with $(j-i-k)\equiv 0\mathrm{mod}p$, if they have the correct distance.

$R_k=\frac{1}{N_k}{\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j=1}^N}\delta \left(d\right)x_i{x}_j(4)$

with

$N_k={\displaystyle \sum _{i=1}^N}{\displaystyle \sum _{j=1}^N}\delta \left(d\right)(5)$

and

$d=\left(j-i-k\right)\mathrm{mod}p(6)$

and with the discrete Dirac function $\delta (d)$ being one if the argument $d$ is zero and being zero otherwise.

For each assumed period $p=1\dots N$, a corresponding correlation function $R_k$ is derived with $k=0\dots p-1$ as described above. Then, a Fourier transform is applied resulting in a power spectrum that corresponds to the specific $p$ chosen.

The spectral peaks are identified and the maximum peak is obtained among all $p$ values investigated. In order to perform a quantitative comparison of the peak heights between the spectra for the various values of $p$, an additional scaling factor is applied to the spectra. Specifically, all spectra are multiplied by their period length $p$, prior to identifying the final maximum $f_{\text{AVP}}$ among the various spectra.

Figure 4B illustrates for three different values $p$ (solid yellow, solid green and dashed gray), how good a periodic continuation of a $p$ steps long sequence fits the measured data. Yellow and green colored boundaries correspond to a bad choice of the sequence length, while dashed gray color corresponds to a better choice.

The concept of explained variance (see Equation 3) is applied to quantify the goodness of fit of the dominant frequency, as previously described for the MaxP algorithm.

3.1.3 Lorentz fit to the amplitude spectrum (Lorentz)

The two aforementioned methods aim to identify the frequency with the greatest power within the estimated power spectra. However, in the context of complex arrhythmia, there may be multiple frequencies that correspond to the movement of different rotors around the heart, making it challenging to identify a single dominant frequency. To address this issue, the authors propose a methodology that examines the trend of the data rather than attempting to describe the arrhythmia with a single frequency. This approach involves fitting a model spectrum to the calculated amplitude spectrum. In order to achieve this, a model that accurately reflects the overall trend of the data is necessary, aligning closely to the peaks obtained from the FFT. The Gaussian, Lorentz, and Voigt functions are useful for modeling spectra, with the Voigt function being the result of convolving a Gaussian and a Lorentz function. In this study, the Lorentz function is applied because of its ability to generate a pronounced peak at the center and elongated wings, allowing for the inclusion of frequencies with diminished amplitude in the fitting process. The Lorentz model is defined as:

$f\left(|X_k|,A,{\mu }_{\text{L}},{\sigma }_{\text{L}}\right)=\frac{A}{\pi }\left(\frac{{\sigma }_{\text{L}}}{{\left(|X_k|-{\mu }_{\text{L}}\right)}^2+{\sigma }_{\text{L}}^2}\right)(7)$

where $|X_k|$ is the amplitude spectrum from the ECG and represents the data to be fitted, $A$ the amplitude of the fitted peak ${\mu }_{\text{L}}$, the center position of the peak and ${\sigma }_{\text{L}}$, the half-width at half-maximum (HWHM). Having defined the model function, the next step is to utilize the non-linear least-squares method for curve fitting to $|X_k|$. In order to promote a more rapid convergence of the algorithm, initial estimates for the amplitude and center can be provided by utilizing the peak power $P(f_{\text{max}})$ of the power spectrum and its corresponding frequency.

The Lorentz model is fitted using the gslnls algorithm (Chau, 2025), and the estimated ECG dominant frequency is provided by $f_{\text{Lorentz}}={\mu }_{\text{L}}$, as shown by the red dot in Figure 4A Additionally, the width of the resulting Lorentz curve $({\sigma }_{\text{L}})$ is related to the frequency content of the ECG. In the case of VT, the curve is typically narrow (small values of ${\sigma }_{\text{L}}$), whereas in the case of VF, it is usually broad (large values of ${\sigma }_{\text{L}}$), as demonstrated in Diaz-Maue et al. (2022).

Similar to Equation 3, the explained variance of the Lorentz fit, as a measure of the goodness of fit, is given by:

${\rho }^2={\displaystyle \sum _k}{f}^2(|X_k|,A,{\mu }_{\mathrm{L}},{\sigma }_{\mathrm{L}})/|X_k{|}^2(8)$

This quantity ranges between 0 and 1, with values of ${\rho }_{\text{L}}^2$ close to 1 indicating that the Lorentz function is a suitable model for the spectrum of the ECG record under consideration.

3.2 Time-domain methods for the estimation of arrhythmia frequency

In this section, the methods for estimating the frequency of arrhythmia as the inverse of the obtained CL, which are included in the R-Library, will be presented. These methodologies are based on the morphological analysis of the recorded ECG waveform.

3.2.1 Sinus fit (SinFit)

As shown in Section 3.1.1, a given signal can be decomposed into a sum of sine and cosine signals by converting it to the frequency-domain using the FFT. Here, the investigated ECG time series is modeled as a sinusoidal function:

$f\left(t,a_0,a_1,a_2,a_3\right)=a_0+a_1\cdot \sin \left(2\pi a_{2}t+a_3\right)(9)$

where $a_0$ represents an offset, $a_1$ the amplitude of the ECG, $a_2$ the modeled frequency and $a_3$ the phase shift. For each test frequency $a_2$, the parameters $a_0,a_1$ and $a_3$ are optimized by linear regression, i.e., by minimizing the squared distance between the ECG signal and the sinusoidal model (see Equation 9). The frequency $a_2$ that corresponds to the maximum explained variance ${\rho }^2={\sigma }^2(f(t,a_0,a_1,a_2,a_3))/{\sigma }_x^2$ (i.e., the time-domain representation of Equation 3) is considered the best fitting frequency $f_{\text{SinFit}}$ with respect to the chosen set of test frequencies.

For the analysis of the murine ECGs, the value of the frequency, $a_2$, can be constrained to 15 Hz $<a_2<$ 40 Hz, as this arrhythmia frequency interval has been previously observed in the mouse heart (Kass et al., 1998; Betsuyaku et al., 2004). In contrast, the frequency range for frequency analysis of the human arrhythmia is 2 Hz $<a_2<$ 7 Hz. In the case of quasi-sinusoidal-shaped ECGs, the explained variance is ${\rho }^2=1$. Moreover, if the test frequencies $a_2$ are equal to the frequencies obtained from the MaxP algortihm, the best fitting frequency $f_{\text{SinFit}}$ coincides with the frequency that corresponds to the maximum of the periodogram.

3.2.2 Average Distance between upward zero-crossings (UZC)

Another intuitive frequency estimation method is based on the analysis of upward (or downward respectively) zero-crossings. Since the signal pre-processing includes mean value subtraction there are multiple zero-crossings. We identify the time points of upward zero-crossings $t_{n_1},t_{n_2},\dots t_{n_k}$ (illustrated as ocher-colored filled triangles in Figure 4C). These correspond to return intervals $T_i=t_{n_{i+1}}-t_{n_i}$. In the case of a periodic signal with a frequency $f$, these return intervals are equal to $1/f$. In the more general case of irregular oscillations the frequency estimator can be expressed as the inverse of the average return interval:

$f_{\text{UZC}}={\left(\frac{1}{k-1}{\displaystyle \sum _{i=1}^{k-1}}{T}_i\right)}^{-1}(10)$

In contrast to the SinFit method, this technique is applicable to a broader range of non-sinusoidal and periodic signals. However, challenges may arise when dealing with polymorphic tachyarrhythmia that are modeled as superpositions of periodic signals. To quantify the uncertainty associated with this method, the coefficient of variation $c_{\text{v}}$ for the intervals between upward zero-crossings is suggested

$c_{\text{v}}=\sigma \left(T\right)/\mu \left(T\right)(11)$

with $\sigma (T)$ being the standard deviation and $\mu (T)$ being the mean of the return intervals. Small values of $c_{\text{v}}$ indicate that the time series has a dominant frequency.

3.2.3 Average Maximum-to-Maximum Distance (M2M)

The arrhythmia time series can be considered as a smooth and oscillating signal, and therefore, the local maxima (or minima, respectively) $t_{n_1},t_{n_2},\dots t_{n_k}$ and the resulting return intervals $T_i=t_{n_{i+1}}-t_{n_i},i=1,\dots ,k-1$ can be identified with peak finding algorithms (illustrated as green filled circles in Figure 4C). In this study, the peaks routine of the R-library splus2R (Constantine and Hesterberg, 2024), was applied with the parameters (span=11, strict=TRUE, endbehavior=0). In the context of a sinusoidal signal with frequency $f$, these maxima are equidistant, with a distance of $1/f$.

Consequently the inverse value of the average temporal distance between succeeding local maxima, $f_{\text{M}2\text{M}}$, which equation is similar to Equation 10, is defined as an additional estimator of arrhythmia frequency. The coefficient of variation $c_{\text{v}}$ is also calculated as a quantitative measure for the M2M uncertainty. In the presence of irregular signals such as pVT or VF, $c_{\text{v}}$ will yield large values.

3.2.4 ICD algorithm (ICD)

In order to facilitate a comparison with a clinically utilized algorithm, the authors implemented the sensing algorithm used by Biotronik in ICDs (Brüggemann et al., 2016), and included it in the provided R-Library. A notable advantage of this method is its versatility, as it can be applied to both sinus rhythm and to arrhythmia cases. The ICD algorithm similar to the M2M method, is based on the determination of the maximum values of the ECG. However, it constrains the maxima to those that have a minimum distance to the preceding maximum and a maximum height of at least a certain percentage of the height of this preceding maximum.

As a brief description of the method, the algorithm involves the use of distinct threshold levels, which are calibrated dynamically in relation to the peak amplitude of the signal over predefined time periods (refractory period). The local maximum of the ECG is measured and employed to establish the upper and lower thresholds. The upper threshold is initially set to 50% of the detected peak and maintained for a specified refractory period. Subsequent to this interval, the lower threshold is set to 25% of the measured peak for an additional defined period. Subsequently, the algorithm proceeds to decrease the sensing threshold by 12.5% for the stipulated period. This process continues until either the minimum threshold is reached or a new maximum value is encountered. More details can be found in the original publication by Brüggemann et al. (2016).

Similar to the M2M and UZC techniques the time intervals between the identified maxima (illustrated as beige filled squares in Figure 4C) will be determined and used to compute the frequency $f_{\text{ICD}}$ and its uncertainty $c_{\text{v}}$. The strength of this technique lies in the incorporation of a refractory period, during which no maxima will be detected. Consequently, a reduced number of identified maxima in complex arrhythmia are observed in comparison to the M2M and UZC techniques, resulting in a lower frequency $f_{\text{ICD}}$.

4 Benchmarks

A comparative analysis of seven distinct algorithms for determining arrhythmia frequency was conducted in this work. Three of the algorithms operate in the frequency-domain using Fourier techniques, while the remaining four exploit the morphological properties of the ECG in the time-domain.

In this section the statistical properties of the seven frequency estimators will be studied by applying them to model data with well-known properties. To this end, two sets of model time series have been generated, one exhibiting regular and the other both regular and irregular oscillations. To ensure the comparability of the model data to the considered murine ECGs, the model data will contain 1,000 data points each.

The first class of model time series contains sinusoidal functions with frequencies $f$ in the range from $15.0\text{\hspace{0.17em}Hz}$ to $30.0\text{\hspace{0.17em}Hz}$. The evaluation of these data demonstrated that all seven techniques achieved the expected result within the limits of accuracy of the specific estimator, e.g., the MaxP technique will result in the integer value closest to the frequency of the model time series. The explained variances of the spectral techniques are equal to or close to 1. Furthermore, the coefficients of variation for the time-domain techniques are close to 0.

The second class of model time series encompasses regular and irregular oscillations, which have been modeled as stochastic oscillations (Gardiner, 2009; Pikovsky et al., 2003). Stochastic oscillators are described by stochastic differential equations. In this instance, the following system has been integrated by applying an Euler–Maruyama (Maruyama, 1955; Kloeden and Platen, 1992) scheme:

$\ddot{x}\left(t\right)=-\gamma \dot{x}\left(t\right)-{\omega }_0^{2}x\left(t\right)+\sigma \zeta \left(t\right)(12)$

with ${\omega }_0$ denoting the average angular frequency, $\sigma$ the amplitude of the stochastic force and $\gamma$ the strength of the damping, which is controlling the irregularity of the oscillation.

A total of 1,000 time series have been generated with a frequency of $f^{*}=22\text{\hspace{0.17em}Hz}$, a noise amplitude $\sigma =1$, and a uniformly distributed $\gamma$ value within the interval (0,15). The length of the time series length was set to $1\text{\hspace{0.17em}s}$ and the sampling frequency to $f_{\text{s}}=1\text{\hspace{0.17em}kHz}$. Moreover, to ensure optimal compatibility with the ECG data, the model data underwent filtering using identical bandpass filter settings. Depending on the value of $\gamma$, the oscillations may exhibit slight or moderate irregularity, as shown in Supplementary Figure S1 (Supplementary Material), which presents three exemplary time series for reference.

The frequency estimation results for the seven algorithms are graphically presented in Figure 5, where panel A shows boxplots of the estimated frequencies along with the results for the statistical tests. Here, the median frequency (i.e., horizontal line inside the box) of all estimators, with the exception of M2M, is found to be approximately 22 Hz. The height of the boxes represents the random error of the estimation techniques in terms of the inter-quartile range, which is particularly large for the M2M and the ICD estimator and small for the Lorentz fit. It is noteworthy that the MaxP algorithm exhibits a coarser resolution, a consequence of the spectral resolution attained with $N=1000$ points.

Figure 5

Figure 5. Results of the seven frequency estimators applied to model time series with irregular oscillations. All time series have a frequency of $f^{*}=22\text{\hspace{0.17em}Hz}$. (A) Boxplots vizualise results of the individual techniques: The frequency estimates are represented by gray dots. Each of the boxes gives the mean value (diamond symbol), the median (horizontal line in middle of the box) as well as 25th and 75th percentiles of the estimated frequencies (upper and lower edges of the box). Significance brackets present the results (adjusted p-values) of pairwise t-tests: * stands for $p<0.05$, ** for $p<0.01$, *** for $p<0.001$ and **** for $p<0.0001$. (B) Pairwise comparison of the seven estimators of arrhythmia frequency: The histograms of the results of the seven techniques are shown along the diagonal. Bland-Altman plots of pairwise arrhythmia frequency estimates are presented in the upper right triangular part of the matrix and the Pearson correlation coefficients with their levels of statistical significance and their graphical representation are given in the lower left triangular part. Figure created using the ggplot2 library of R (Wickham et al. (2019)

A series of statistical tests (Fisher, 1925; Montgomery, 2017), were applied to determine the significant differences between the outcomes derived from the seven techniques. Initially, a repeated measures analysis of variance (ANOVA) with a Greenhouse-Geisser correction was conducted. This analysis indicated statistically significant differences between the frequency estimators $(p<10^{-4})$. Given the findings of the Shapiro–Wilk test, which indicated violations of the normality assumption for all estimators, a nonparametric Friedman test was performed as a complement. This additional analysis yielded a value of $p<10^{-4}$, thereby confirming the results of the ANOVA test.

Subsequently, pairwise t-tests with Bonferroni corrections were conducted, to enable a comparison between the individual frequency estimators. The results (see Figure 5A) demonstrated that the spectral estimators, in conjunction with FitSin, did not differ significantly from one another. However, a significant difference was observed between the group of frequency-domain estimators and FitSin, and the UZC, M2M, and ICD estimators. Additionally, the latter three estimators were also found to be mutually distinct. These findings were then corroborated by pairwise Wilcoxon signed-rank tests.

The pairwise agreement between the estimators was visualized using Bland-Altman plots (upper triangular matrix of Figure 5B). Here, each scatter point represents the difference between the paired measurements plotted against the average value derived from the identical time series. The horizontal lines in each Bland-Altman plot represent the mean and the 95% confidence intervals of the differences (calculated as the 2.5th and the 97.5th percentile). In the event that all points were situated along the x-axis, the two considered techniques would exhibit a high degree of agreement. With the exception of a few outliers, this is the case when a comparison is made between AVP and FitSin. An offset in the Bland-Altman plot is indicative of a systematic difference or fixed bias between the considered techniques, as evidenced in all comparisons to M2M.

Finally, the Pearson correlation coefficients for pairwise techniques were calculated and are presented in the lower triangular matrix of Figure 5, panel B. All correlation coefficients, with the exception of the one between M2M and ICD, are determined to be significantly positive ($p$-values based on bootstrapping).

Large correlation coefficients between the spectral estimators indicate that they result in higher (or lower, respectively) than average frequencies when applied to the same model time series. In contrast, small correlation coefficients are indicative of more independent frequency estimation results, as evidenced for all pairs of frequency-domain techniques and UZC, M2M, and ICD.

In summary, the performance of the frequency estimators was benchmarked, and it was demonstrated that in case of regular oscillations all techniques yielded the expected outcome within the limits of the accuracy of the specific estimator. In the context of time series with irregular oscillations, the application of the frequency-domain estimators in conjunction with the FitSin method ensures the absence of bias, while exhibiting a high degree of inter-correlation. Notably, the Lorentz-fit approach exhibits the smallest random error among all considered methods. Additionally, it was also demonstrated that the UZC, M2M and ICD techniques are weakly correlated with the spectral techniques and FitSin.

5 Results

The performance of the seven algorithms for measured data was evaluated using a set of 112 ECGs from perfused arrhythmic murine hearts and was further applied to a short set of human arrhythmia.

5.1 Mouse data

The results for the frequency estimation of the mouse data are summarized in Figure 6 in the same manner as the results for the benchmarks were visualized in Figure 5. Figure 6A presents boxplots for the results of the seven frequency estimation techniques and the results of the corresponding statistical tests. The estimated frequencies range from 7.5 Hz to 36.5 Hz with the results of the M2M technique being consistently larger than those obtained by the other six techniques: The main frequency ${\bar{f}}_{\text{M}2\text{M}}=26.0\text{\hspace{0.17em}Hz}$ is at least $3\text{\hspace{0.17em}Hz}$ (equivalent to 13%) higher than the mean frequencies observed in the other estimators. The Lorentz fit yielded the smallest mean frequency of ${\bar{f}}_{\text{Lorentz}}=21.2\text{\hspace{0.17em}Hz}$. The ANOVA test yielded p-values of $p<10^{-4}$, indicating statistically significant differences between the seven methods. Subsequently, pairwise t-tests were performed, whose results are displayed as test brackets above the corresponding boxplots. The t-test results indicate similarities in the mean values of the frequency-domain estimators MaxP, AVP, Lorentz, and the FitSin (${\bar{f}}_{\text{MaxP}}=22.2\text{\hspace{0.17em}Hz}$, ${\bar{f}}_{\text{AVP}}=22.5\text{\hspace{0.17em}Hz}$, ${\bar{f}}_{\text{Lorentz}}=21.2\text{\hspace{0.17em}Hz}$ and ${\bar{f}}_{\text{FitSin}}=22.5\text{\hspace{0.17em}Hz}$), whereas the mean values of the remaining estimators are significantly different (${\bar{f}}_{\text{UZC}}=23.04\text{\hspace{0.17em}Hz}$, ${\bar{f}}_{\text{M}2\text{M}}=26.0\text{\hspace{0.17em}Hz}$ and ${\bar{f}}_{\text{ICD}}=23.1\text{\hspace{0.17em}Hz}$). In addition, the interquartile range, as illustrated by the height of the boxes in Figure 6A, demonstrates that the ICD technique yields the least dispersion among the seven examined techniques, with a value of $5.81\text{\hspace{0.17em}Hz}$. The two spectral techniques MaxP, AVP and the time-domain method FitSin yield almost identical results for the majority of ECGs.

Figure 6

Figure 6. Results of the seven frequency estimators applied to the dataset with 112 ECGs of arrhythmic murine hearts. (A) Boxplots of the seven techniques: The frequency estimates are represented by gray dots. Each of the boxes gives the mean value (diamond symbol), the median (horizontal line in middle of the box) as well as 25th and 75th percentiles of the estimated frequencies (box upper and lower edges of the box). Significance brackets present the results (adjusted p-values) of pairwise t-tests: * stands for $p<0.05$, ** for $p<0.01$, *** for $p<0.001$ and **** for $p<0.0001$. (B) Pairwise comparison of the seven estimators of arrhythmia frequency: The histograms of the results of the seven techniques are shown along the diagonal. Bland-Altman plots of pairwise arrhyhtmia frequency estimates are presented in the upper right triangular part of the matrix and the Pearson correlation coefficients with their levels of statistical significance and their graphical representation are shown in the lower left triangular part. Figure created using the ggplot2 library of R (Wickham et al. (2019).

Figure 6B shows the distribution of the estimated frequencies for each method along the diagonal, Bland-Altman plots, which are used to assess the agreement between methods, are displayed above the diagonal, and Pearson’s correlation coefficients and their visualizations are shown below the diagonal. As previously discussed, the three frequency-domain methods and the FitSin exhibit strong agreement as the points in the Bland-Altman are grouped horizontally around the line $y=0$ with only a few outliers. In instances where a method from the frequency-domain or SinFit is compared against time-domain techniques, a proportional bias is observed, manifesting as a downward or upward slope of the scatter plot points.

In agreement with the results shown in Figure 6A, the analysis in Figure 6B further demonstrates that the M2M algorithm estimates larger frequencies $f_{\text{M}2\text{M}}$ than the six other techniques. This is evidenced by a fixed systematic bias as displayed by the horizontal solid cyan line, which is shifted toward the negative direction in subpanel 6g, otherwise shifted in the positive direction for all subpanels in column f. In addition, the M2M method exhibited good alignment with UZC and ICD, demonstrating low variability, narrow limits of agreement, and elevated correlation coefficients (0.79 and 0.81). On the other hand, the limits of agreement are more extensive in comparison with those of MaxP, AVP, Lorentz, and FitSin. This discrepancy is further confirmed by the low correlation coefficients observed in the respective cases.

Figure 6B also shows a strong correlation $(0.83)$ between the MaxP, AVP, and FitSin methods, suggesting they use related information in the recorded signals, specifically a dominant harmonic oscillation. The AVP and FitSin methods have Bland-Altman plots where all means align at 0 and the correlation coefficient is 1, indicating that these two methods can be used interchangeably.

In contrast, when comparing these techniques with the Lorentz fit, a stronger variability is observed. Furthermore the Lorentz fit tends to yield larger frequency estimates for $f_{\text{FitSin}}<15.0\text{\hspace{0.17em}Hz}$ and smaller frequencies $f_{\text{Lorentz}}$ for $f_{\text{FitSin}}>20.0\text{\hspace{0.17em}Hz}$. The Bland-Altman plots 1c, 2c, and 3d reveal a substantial number of dots along zero, indicating that the results are generally similar, except in complex cases where the signal contains multiple frequencies. In such instances, the assumption of a single dominant frequency is inappropriate. Consequently, the MaxP, AVP, and FitSin algorithms fail to handle the superimposed combination of multiple frequencies, while the Lorentz algorithm can effectively address this challenge by using a wider distribution of the function fitted to the amplitude spectrum. Notably, the correlation coefficients between the Lorentz fit and the other three methods are lower, with values of 0.79 and 0.78, respectively.

A direct performance comparison of frequency- and time-domain techniques reveals a moderate degree of correlation between the frequency-domain techniques, the FitSin and the ICD technique, with values of 0.76 and 0.81. Notably, the ICD method and the Lorentz fit showed the strongest correlation among the aforementioned methods, with a correlation coefficient of 0.84. Furthermore, the UZC technique, which is relatively simple, also shows high correlation coefficients with the spectral techniques, FitSin, and the ICD estimator. However the Bland-Altman plots show a proportional bias against all other methods. Conversely, the M2M time-domain technique exhibits lower correlations with spectral techniques, with all corresponding coefficients falling below 0.63.

The performance of the seven estimators was further analyzed for the murine ECGs with less complex arrhythmia, such as ventricular tachycardia (VT). Here, the identification of the arrhythmia frequency can be achieved with a high degree of precision through the application of both time- and frequency-domain methods. The 36 ECG time series with high values of the explained variance with respect to the MaxP estimator $({\rho }^2>0.5)$ were identified, and the correlation coefficients exceeded 0.95 for all pairs of techniques, with the exception of the M2M method. A similar set of results was obtained in 31 ECG time series, in which the time intervals between the maxima identified by the ICD algorithm exhibited a low coefficient of variation $(c_{\text{v}}<0.2)$. Consequently, the CLs exhibited a high degree of similarity, with all correlation coefficients exceeding 0.90. This finding confirms the results achieved for sinusoidal functions (Section 4).

5.2 Human data

The analysis of human arrhythmia data revealed frequencies of $f_{\text{VT}}=2.61\text{pm\hspace{0.17em}}0.04\text{\hspace{0.17em}Hz}$, $f_{\text{VFlu}}=3.53\text{pm\hspace{0.17em}}0.03\text{\hspace{0.17em}Hz}$, and $f_{\text{VF}}=4.82\text{pm\hspace{0.17em}}0.82\text{\hspace{0.17em}Hz}$ for events of VT, VFlu and VF respectively. It is noteworthy that all methods produced nearly identical results for the instances of VT and VFlu (as shown in Table 1), while the inter-method variability was markedly higher for VF. In this case, M2M yielded the largest frequency estimate.

Table 1

Table 1. Frequency estimates and the corresponding dispersions (explained variance and coefficient of variation) for the three examples of human arrhythmic ECGs shown in Figure 3.

In the case of VT, the low coefficient of variations for the time-domain methods $(c_{\text{v}}<0.04)$, and the high explained variances for the frequency-domain methods and the SinFit ${\rho }^2>0.8$ suggest the presence of a single dominant frequency. The analysis of VFlu demonstrates high degree of agreement between methods while the findings in Table 1 substantiate the single dominant frequency hypothesis as well. Accordingly, as stated by Viskin et al. (2003), VFlu is characterized by sinusoidal morphology, which is driven by a rotating spiral with a defined frequency (Davtdenko, 1993; Comtois et al., 2005; Lim et al., 2006). Conversely, the frequency estimates for VF range from 4.2 Hz to 6.35 Hz, while the explained variances fall below 0.27 and the coefficient of variation increases to more than 0.26. These findings suggest that the arrhythmia may not be adequately defined by a single dominant frequency.

6 Discussion

The analysis of ECG signals from arrhythmic hearts often requires the identification of heart rates, dominant frequencies or cycle lengths for diagnostic or therapeutic purposes. Arrhythmic ECGs exhibit dynamic, repeating, and non-repeating patterns that can be analyzed using the algorithms presented in this study. These algorithms are designed to quantify the periodicity of ECG patterns, thereby yielding a numerical value that represents the dominant frequency or dominant period. A thorough examination of the mouse dataset revealed that the algorithms presented herein yielded distinct outcomes, particularly when applied to complex signals. This divergence does not necessarily imply an issue with the performance of the algorithms. Rather, it highlights a more fundamental question: Does the information contained within the ECG data genuinely correspond to a single dominant frequency or period?

To illustrate this point, consider a signal that initially comprises two harmonic frequencies ($f_1$ and $f_2=2f_1$) with the same amplitude. In this case, the waveform will exhibit a clear repeating pattern, and the identified frequency will unambiguously correspond to a fundamental period $(T_1=1/f_1)$. This phenomenon can be attributed to the repetitive nature of the signal pattern, which resumes after a specific interval, $T_1$. This characteristic is commonly observed in VT, particularly in cases of mVT. The frequency-domain techniques previously introduced will identify the nature of the signal, one of the two frequencies will be identified. In contrast, some of the time-domain techniques previously examined will lead to $f_2$ as dominant frequency. In particular, we find $f_{\text{M}2\text{M}}=f_2$ and $f_{\text{UZC}}=f_2$. Now, introducing a third frequency, $f_3$, to the original signal, which is not in resonance with $f_1$ but has a comparable amplitude, would significantly alter the overall waveform. This phenomenon manifests, for instance, when a stationary spiral begins to meander within the cardiac substrate, as it is observed in pVT. Consequently, the new dominant period determined with the presented algorithms would shift to a rational multiple of the previous $T_1$. The presence of intermittent dynamic patterns, such as those observed during ventricular fibrillation, complicates further the determination of a dominant frequency from an ECG signal.

In this study, we also proposed two quantitative measures for identifying the presence of a dominant frequency. The first involves the computation of the explained variance, denoted by ${\rho }^2$, for the three spectral techniques and the sinus fit method. Values of ${\rho }^2$ that are close to one represent ECG signals which can be accurately modeled by sinusoidal functions and can be therefore represented by a single repeating pattern. In contrast, non-sinusoidal ECG signals, despite their periodicity, will result in lower values of ${\rho }^2$. Therefore, it is also crucial to consider the second quantitative measure, the coefficient of variation $c_{\text{v}}$ of the time interval lengths between the maxima of the ECG signal, as it is calculated for the M2M, ICD and UZC methods. A narrow range of these interval lengths demonstrates a regular oscillation and thus the presence of a dominant frequency even in the case of non-sinusoidal oscillations. In this context, the assumption of a single dominant frequency is challenged by complex rhythms, such as atrial fibrillation (AF), which are characterized by multifrequency spectra resulting from multiple interacting rotors (Sanders et al., 2005; Gadenz et al., 2017). The usage of the two independent quantitative measures presented in this work enables the rigorous evaluation of the validity of the single dominant frequency assumption in any given ECG segment. By considering both ${\rho }^2$ and $c_{\text{v}}$, a robust framework is provided to determine when a single frequency is a suitable descriptor for the arrhythmia segment, thereby strengthening the interpretation of the results beyond a simple frequency estimation. This approach is essential, given the morphological variability observed in cardiac arrhythmia.

It was also observed that time-domain-based algorithms show a substantial dispersion of frequency estimates, resulting in a more pronounced mutual deviation when compared to the frequency-domain methods. This discrepancy is attributable to their ability to identify localized patterns within the ECG, while neglecting the overall course of the signal. As discussed in Section 5, the M2M algorithm results in the largest frequency estimate for the time-domain methods. This can be explained as follows: Each threshold crossing is associated with a local maximum; however, because multiple local maxima may occur above the threshold and multiple local minima may occur below it, the total number of local maxima can exceed the number of threshold crossings. Therefore, the frequency estimates derived from M2M are larger than or equal to the UZC-based estimates. The frequency estimate derived from the ICD techniques excludes maxima, which occur during the refractory period, resulting in a value that is smaller than or equal to the M2M estimate.

In the frequency-domain methods, Fourier analysis has been incorporated into the MaxP algorithm, as this method has a long history of application in ECG frequency determination (Nolle et al., 1980; Herbschleb et al., 1980; Caldwell et al., 2007). In this case it is important to consider that when a segment of a continuous signal undergoes Fourier analysis, there is a high probability that the selected segment will contain incomplete periods. It is inevitable that incomplete periods will result in leakage effects, which, in turn, will lead to a widened Fourier spectrum. Consequently, pinpointing a single sharp peak corresponding to a dominant frequency becomes challenging. In previous studies, alternative methods such as the Organization Index (Everett et al., 2001), Organization Analysis (Barquero-Perez et al., 2010), Spectral Flatness (Eftestøl et al., 2000), and Ensemble Average (Ciaccio et al., 2011) have been employed to address this issue. While the application of a window function is another recognized method applied for mitigating the leakage effect, we opted against it to avoid the introduction of an additional amplitude modulation effect. We also present a novel method for addressing the aforementioned problem, namely, the AVP algorithm, which can handle an assumed fundamental period of the signal under investigation independent of the duration of the signal recorded, calculating corresponding FFTs for various signal periods assumed. This process identifies the maximum peak of the spectrum that corresponds to an integer multiple of the fundamental period, and consequently to the dominant frequency.

The SinFit and Lorentz algorithms were developed based on fitting techniques. While the dispersion of the frequency values found with the SinFit algorithm appears to perform equally to the MaxP and AVP methods, it was observed that the frequencies found with the Lorentz technique are more similar to the time-domain techniques in particular to the ICD and UZC methods. This is because the fitting can better manage the superposition of two or more frequencies in the ECG signal.

Within the time-domain methods, the developed R-library includes an algorithm for identifying heart rate, based on the ICD algorithm used by Brüggemann et al. (2016). Here it is important to note that the intracardiac electrogram (EGM) obtained via an intravenous catheter exhibits a distinct frequency spectrum (Requena-Carrión et al., 2013) when compared to the frequency observed in surface ECGs. Consequently, it is necessary to apply proper signal preconditioning to EGM data, including bandpass filtering (different settings for unipolar and bipolar EGMs), and under certain circumstances rectification (Ng and Goldberger, 2007), prior to applying the remaining six algorithms presented in this study. ECG data from arrhythmic hearts can be analyzed directly using the ICD method and the other algorithms presented. Although the pre-processing of surface ECG varies slightly between species, these variations do not impede the direct application of the aforementioned methods.

The remaining two time-domain methods M2M and UZC are particularly well-suited for scenarios where low computational power is a prerequisite, such as in lower-architecture microcontroller-aided applications. As demonstrated in our results, the dominant period obtained is of a diminished duration, leading to a dominant frequency that exceeds those yielded by the alternative algorithms. Nonetheless, these two methods have been shown to provide rapid performance and a reasonably accurate estimate of the dominant frequency.

In instances where the estimation of a dominant frequency is necessary prior to the application of a particular arrhythmia therapy, such as overdrive pacing (for VT (Ellenbogen et al., 2008), for VF (Luther et al., 2011)) or multistage therapies (Li et al., 2009), our findings suggest that the utilization of both time- and frequency-domain algorithms ensures a more precise frequency estimation for complex arrhythmia. However, a critical consideration is that the ECG signal of an arrhythmic heart undergoes constant changes. Therefore, it would be prudent to consider the implementation of the presented algorithms within sliding windows, thereby facilitating the monitoring of dynamic changes.

The present study is constrained to the statistical comparison of the seven methods for frequency estimation that are currently or that could be employed in real-time, low-power processing environments, such as those found in implantable cardioverter-defibrillators (ICDs). These life-critical embedded systems require algorithms that are mathematically simple, highly reliable, and computationally efficient. Nonetheless, the implementation of alternative signal analysis techniques, such as wavelet analysis (Grossmann and Morlet, 1984; Daubechies et al., 2011; Zhuravlev et al., 2025), which offer powerful capabilities for characterizing non-stationary events and the dynamic frequency components associated with complex arrhythmia transitions, remains a promising avenue for future research. In addition, recent studies have demonstrated the efficacy of advanced signal processing methodologies, such as machine learning (ML) algorithms, empirical mode decomposition (EMD), and sophisticated autoregressive models, in the domain of biosignal processing (Mohanty et al., 2021; Ahmadi and Ekhlasi, 2019; Baselli et al., 1985). Subsequent research endeavors will involve the integration and comparison of other analysis approaches within this analytical framework.

7 Conclusion

The objective of this study was to provide interested readers and the scientific community with a comparative analysis of algorithms that determine the dominant frequency of ECGs of arrhythmic hearts. We have investigated seven time- and frequency-domain techniques, which are made available in the R-library FibFreq. The R-library has been demonstrated to be suitable for the analysis of both an experimentally obtained mouse ECG dataset and a human dataset (MIT-BIH-malignant-ventricular-ectopy-database-1.0.0). Despite the awareness that heart rates vary significantly among different species (Schüttler et al., 2020), we have demonstrated that by employing the appropriate pre-processing techniques and utilizing our R-Library, it is possible to ascertain a reliable estimate for the dominant frequency for less complex arrhythmia. Moreover, the R-Library can be also applied to other oscillatory phenomena, as demonstrated in Supplementary Material for the frequency determination of sunspots.

In addition to providing a frequency estimate, the algorithms under consideration also yield secondary quantities that can be used to evaluate whether the given signal can be adequately modeled using a single periodic pattern. These are the explained variance for all frequency-domain techniques and the SinFit, and the coefficient of variation of the identified time intervals for the three time-domain techniques. The data analysis presented herein indicates that a low coefficient of variation or an explained variance that approaches to one are indicative of a signal that is essentially periodic and can be described using a single dominant frequency. Therefore, in order to ascertain the presence of a single dominant frequency in an ECG time series, it is recommended that both time-domain and frequency-domain techniques be applied.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. The repository for the R-Library: https://github.com/origami-dss/FibFreq.

Ethics statement

Ethical approval was not required for the study involving humans in accordance with the local legislation and institutional requirements. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements. The animal study was approved by Federation of European Laboratory Animal Science Associations. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

LD-M: Conceptualization, Data curation, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing. AW: Data curation, Formal Analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review and editing. LE: Investigation, Writing – review and editing. HN: Formal Analysis, Investigation, Methodology, Validation, Writing – original draft, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This work was supported by the Max Planck Society and by the DZHK (German Centre for Cardiovascular Research), funding code: 81Z0300114.

Acknowledgements

The authors want to thank Marion Kunze and Tina Althaus for their excellent technical assistance during experimental setup and experiment conduction. Many thanks to Claudia Richter and Stefan Luther for fruitful discussions and thought-provoking impulses. We also thank the reviewers, who provided in-depth reviews and comments which substantially improved this manuscript.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. DeepL AI Labs (write function for editing).

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsip.2025.1707422/full#supplementary-material

References

Ahmadi, H., and Ekhlasi, A. (2019). “Types of emd algorithms,” in 2019 5th Iranian conference on signal processing and intelligent systems (ICSPIS), 1–5. doi:10.1109/ICSPIS48872.2019.9066155

CrossRef Full Text | Google Scholar

Ambrosi, C. M., Ripplinger, C. M., Efimov, I. R., and Fedorov, V. V. (2011). Termination of sustained atrial flutter and fibrillation using low-voltage multiple-shock therapy. Heart rhythm. 8, 101–108. doi:10.1016/j.hrthm.2010.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Annoni, E. M., Arunachalam, S. P., Kapa, S., Mulpuru, S. K., Friedman, P. A., and Tolkacheva, E. G. (2017). Novel quantitative analytical approaches for rotor identification and associated implications for mapping. IEEE Trans. Biomed. Eng. 65, 273–281. doi:10.1109/TBME.2017.2763460

PubMed Abstract | CrossRef Full Text | Google Scholar

Barquero-Perez, O., Rojo-Alvarez, J. L., Caamano, A. J., Goya-Esteban, R., Everss, E., Alonso-Atienza, F., et al. (2010). Fundamental frequency and regularity of cardiac electrograms with fourier organization analysis. IEEE Trans. Biomed. Eng. 57, 2168–2177. doi:10.1109/TBME.2010.2049574

PubMed Abstract | CrossRef Full Text | Google Scholar

Baselli, G., Bolis, D., Cerutti, S., and Freschi, C. (1985). Autoregressive modeling and power spectral estimate of rr interval time series in arrhythmic patients. Comput. Biomed. Res. 18, 510–530. doi:10.1016/0010-4809(85)90027-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Betsuyaku, T., Kanno, S., Lerner, D. L., Schuessler, R. B., Saffitz, J. E., and Yamada, K. A. (2004). Spontaneous and inducible ventricular arrhythmias after myocardial infarction in mice. Cardiovasc. Pathol. 13, 156–164. doi:10.1016/S1054-8807(03)00152-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Brüggemann, T., Dahlke, D., Chebbo, A., and Neumann, I. (2016). Tachycardia detection in modern implantable cardioverter–defibrillators. Herzschrittmachertherapie + Elektrophysiologie 27, 171–185. doi:10.1007/s00399-016-0449-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Cakulev, I., Wit, A. L., and Waldo, A. L. (2021). The history of mapping. Cham: Springer International Publishing, 27–40. doi:10.1007/978-3-030-63355-4_3

CrossRef Full Text | Google Scholar

Caldwell, J., Burton, F. L., Smith, G. L., and Cobbe, S. M. (2007). Heterogeneity of ventricular fibrillation dominant frequency during global ischemia in isolated rabbit hearts. Journal of Cardiovascular Electrophysiology 18, 854–861. doi:10.1111/j.1540-8167.2007.00867.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chau, J. (2025). Gslnls: GSL multi-start nonlinear least-squares fitting.

Google Scholar

Ciaccio, E. J., Coromilas, J., Wit, A. L., and Garan, H. (2011). Onset dynamics of ventricular tachyarrhythmias as measured by dominant frequency. Heart Rhythm 8, 615–623. doi:10.1016/j.hrthm.2010.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Clayton, R. H., Murray, A., and Campbell, R. W. (1993). Comparison of four techniques for recognition of ventricular fibrillation from the surface ecg. Medical and Biological Engineering and Computing 31, 111–117. doi:10.1007/BF02446668

PubMed Abstract | CrossRef Full Text | Google Scholar

Comtois, P., Kneller, J., and Nattel, S. (2005). Of circles and spirals: bridging the gap between the leading circle and spiral wave concepts of cardiac reentry. EP Europace 7, S10–S20. doi:10.1016/j.eupc.2005.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Constantine, W., and Hesterberg, T. (2024). splus2R: Supplemental S-PLUS functionality in R. R Package Version 1, 3–5. Available online at: https://github.com/spkaluzny/splus2r.

Google Scholar

Cooley, J. W., and Tukey, J. W. (1965). An algorithm for the machine calculation of complex fourier series. Mathematics of Computation 19, 297–301. doi:10.2307/2003354

CrossRef Full Text | Google Scholar

Daubechies, I., Lu, J., and Wu, H.-T. (2011). Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool. Applied and Computational Harmonic Analysis 30, 243–261. doi:10.1016/j.acha.2010.08.002

CrossRef Full Text | Google Scholar

Davtdenko, J. M. (1993). Spiral wave activity: a possible common mechanism for polymorphic and monomorphic ventricular tachycardias. Journal of Cardiovascular Electrophysiology 4, 730–746. doi:10.1111/j.1540-8167.1993.tb01258.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Diaz-Maue, L., and Nobach, H. (2022). “Correlation function and power spectrum from almost periodic signals with incomplete periods. In,” in 29. Fachtagung experimentelle strömungsmechanik (german association for laser anemometry GALA e.V.,), 33, 1–33.9

Google Scholar

Diaz-Maue, L., Steinebach, J., and Richter, C. (2022). Patterned illumination techniques in optogenetics: an insight into decelerating murine hearts. Frontiers in Physiology 12, 750535. doi:10.3389/fphys.2021.750535

PubMed Abstract | CrossRef Full Text | Google Scholar

Eftestøl, T., Sunde, K., Ole Aase, S., Husøy, J. H., and Steen, P. A. (2000). Predicting outcome of defibrillation by spectral characterization and nonparametric classification of ventricular fibrillation in patients with out-of-hospital cardiac arrest. Circulation 102, 1523–1529. doi:10.1161/01.cir.102.13.1523

PubMed Abstract | CrossRef Full Text | Google Scholar

Ellenbogen, K. A., Wood, M. A., and Gold, M. R. (2008). “The implantable cardioverter defibrillator,” in Cardiac pacing and ICDs Editors K. A. Ellenbogen, and M. A. Wood 5 edn. (Malden, MA: Blackwell Publishing), 382–414.

Google Scholar

Everett, T. H., IV, Moorman, J. R., Kok, L.-C., Akar, J. G., and Haines, D. E. (2001). Assessment of global atrial fibrillation organization to optimize timing of atrial defibrillation. Circulation 103, 2857–2861. doi:10.1161/01.cir.103.23.2857

PubMed Abstract | CrossRef Full Text | Google Scholar

Fenton, F. H., Luther, S., Cherry, E. M., Otani, N. F., Krinsky, V., Pumir, A., et al. (2009). Termination of atrial fibrillation using pulsed low-energy far-field stimulation. Circulation 120, 467–476. doi:10.1161/CIRCULATIONAHA.108.825091

PubMed Abstract | CrossRef Full Text | Google Scholar

Fisher, R. A. (1925). Statistical methods for research workers. Edinburgh: Oliver & Boyd.

Google Scholar

Gadenz, L., Hashemi, J., Shariat, M. H., Gula, L., and Redfearn, D. P. (2017). Clinical role of dominant frequency measurements in atrial fibrillation Ablation–a systematic review. Journal of Atrial Fibrillation 9, 1548. doi:10.4022/jafib.1548

PubMed Abstract | CrossRef Full Text | Google Scholar

Gardiner, C. W. (2009). Handbook of stochastic methods: for physics, chemistry and the natural sciences. (Berlin: Springer), 4 edn.

Google Scholar

George, S. A., and Efimov, I. R. (2019). Optocardiography: a review of its past, present, and future. Current Opinion in Biomedical Engineering 9, 74–80. doi:10.1016/j.cobme.2019.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, R. A., Jalife, J., Panfilov, A., Baxter, W. T., Cabo, C., Davidenko, J. M., et al. (1995a). Nonstationary vortexlike reentrant activity as a mechanism of polymorphic ventricular tachycardia in the isolated rabbit heart. Circulation 91, 2454–2469. doi:10.1161/01.cir.91.9.2454

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, R. A., Jalife, J., Panfilov, A. V., Baxter, W. T., Cabo, C., Davidenko, J. M., et al. (1995b). Mechanisms of cardiac fibrillation. Science 270, 1222–1223. doi:10.1126/science.270.5239.1222

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenwald, S. D. (1986). The development and analysis of a ventricular fibrillation detector. Ph.D. thesis, Cambridge, United States: Massachusetts Institute of Technology.

Google Scholar

Grossmann, A., and Morlet, J. (1984). Decomposition of hardy functions into square integrable wavelets of constant shape. SIAM Journal on Mathematical Analysis 15, 723–736. doi:10.1137/0515056

CrossRef Full Text | Google Scholar

Herbschleb, J., Heethaar, R., Tweel, I., and Meijler, F. (1980). Frequency analysis of the ecg before and during ventricular fibrillation. Computers in Cardiology, 365–368. Available online at: https://dspace.library.uu.nl/handle/1874/209848.

Google Scholar

Hornung, D., Biktashev, V., Otani, N., Shajahan, T., Baig, T., Berg, S., et al. (2017). Mechanisms of vortices termination in the cardiac muscle. Royal Society Open Science 4, 170024. doi:10.1098/rsos.170024

PubMed Abstract | CrossRef Full Text | Google Scholar

Jalife, J., and Gray, R. (1996). Drifting vortices of electrical waves underlie ventricular fibrillation in the rabbit heart. Acta Physiologica Scandinavica 157, 123–132. doi:10.1046/j.1365-201X.1996.505249000.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, Y. C., Uzelac, I., Otani, N., Luther, S., Gilmour Jr, R. F., Cherry, E. M., et al. (2017). Synchronization as a mechanism for low-energy anti-fibrillation pacing. Heart Rhythm 14, 1254–1262. doi:10.1016/j.hrthm.2017.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Kass, D. A., Hare, J. M., and Georgakopoulos, D. (1998). Murine cardiac function: a cautionary tail. Circulation Research 82, 519–522. doi:10.1161/01.res.82.4.519

PubMed Abstract | CrossRef Full Text | Google Scholar

Kastor, J. A., Spear, J. F., and Moore, E. N. (1972). Localization of ventricular irritability by epicardial mapping. Circulation 45, 952–964. doi:10.1161/01.CIR.45.5.952

PubMed Abstract | CrossRef Full Text | Google Scholar

Kinoshita, O., Kamakura, S., Ohe, T., Yutani, C., Matsuhisa, M., Aihara, N., et al. (1992). Spectral analysis of signal-averaged electrocardiograms in patients with idiopathic ventricular tachycardia of left ventricular origin. Circulation 85, 2054–2059. doi:10.1161/01.cir.85.6.2054

PubMed Abstract | CrossRef Full Text | Google Scholar

Kloeden, P. E., and Platen, E. (1992). in Stochastic modelling and applied probability (Springer). doi:10.1007/978-3-662-12616-5

CrossRef Full Text | Google Scholar

Li, W., Ripplinger, C. M., Lou, Q., and Efimov, I. R. (2009). Multiple monophasic shocks improve electrotherapy of ventricular tachycardia in a rabbit model of chronic infarction. Heart Rhythm 6, 1020–1027. doi:10.1016/j.hrthm.2009.03.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Lim, Z. Y., Maskara, B., Aguel, F., Emokpae Jr, R., and Tung, L. (2006). Spiral wave attachment to millimeter-sized obstacles. Circulation 114, 2113–2121. doi:10.1161/CIRCULATIONAHA.105.598631

PubMed Abstract | CrossRef Full Text | Google Scholar

Luther, S., Fenton, F. H., Kornreich, B. G., Squires, A., Bittihn, P., Hornung, D., et al. (2011). Low-energy control of electrical turbulence in the heart. Nature 475, 235–239. doi:10.1038/nature10216

PubMed Abstract | CrossRef Full Text | Google Scholar

Maria, E. D., Giacopelli, D., Borghi, A., Modonesi, L., and Cappelli, S. (2017). Antitachycardia pacing programming in implantable cardioverter defibrillator: a systematic review. World Journal of Cardiology 9, 429–436. doi:10.4330/wjc.v9.i5.429

PubMed Abstract | CrossRef Full Text | Google Scholar

Maruyama, G. (1955). Continuous markov processes and stochastic equations. Rendiconti del Circolo Matematico di Palermo 4, 48–90. doi:10.1007/BF02846028

CrossRef Full Text | Google Scholar

Mohanty, M., Dash, M., Biswal, P., and Sabut, S. (2021). Classification of ventricular arrhythmias using empirical mode decomposition and machine learning algorithms. Progress in Artificial Intelligence 10, 489–504. doi:10.1007/s13748-021-00250-6

CrossRef Full Text | Google Scholar

Montgomery, D. C. (2017). Design and analysis of experiments (Hoboken, NJ: John Wiley and Sons), 9 edn.

Google Scholar

Moody, G. B., and Mark, R. G. (2001). The impact of the mit-bih arrhythmia database. IEEE Engineering in Medicine and Biology Magazine 20, 45–50. doi:10.1109/51.932724

PubMed Abstract | CrossRef Full Text | Google Scholar

Moss, A. J., Hall, W. J., Cannom, D. S., Daubert, J. P., Higgins, S. L., Klein, H., et al. (1996). Improved survival with an implanted defibrillator in patients with coronary disease at high risk for ventricular arrhythmia. N. Engl. Journal of Medicine 335, 1933–1940. doi:10.1056/NEJM199612263352601

PubMed Abstract | CrossRef Full Text | Google Scholar

Moss, A. J., Zareba, W., Hall, W. J., Klein, H., Wilber, D. J., Cannom, D. S., et al. (2002). Prophylactic implantation of a defibrillator in patients with myocardial infarction and reduced ejection fraction. N. Engl. Journal of Medicine 346, 877–883. doi:10.1056/NEJMoa013474

PubMed Abstract | CrossRef Full Text | Google Scholar

Ng, J., and Goldberger, J. J. (2007). Understanding and interpreting dominant frequency analysis of af electrograms. Journal of Cardiovascular Electrophysiology 18, 680–685. doi:10.1111/j.1540-8167.2007.00832.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Nolle, F., Ryschon, K., and Zeencka, A. (1980). Power spectrum analysis of ventricular fibrillation and imitative artifacts. Comput. Cardiol 209, 212.

Google Scholar

Pandit, S. V., and Jalife, J. (2013). Rotors and the dynamics of cardiac fibrillation. Circulation Research 112, 849–862. doi:10.1161/CIRCRESAHA.111.300158

PubMed Abstract | CrossRef Full Text | Google Scholar

Pikovsky, A., Rosenblum, M., and Kurths, J. (2003). Synchronization: a universal concept in nonlinear sciences. in Cambridge Nonlinear Science series (Cambridge: Cambridge University Press).

Google Scholar

Requena-Carrión, J., Alonso-Atienza, F., Everss, E., Sánchez-Muñoz, J. J., Ortiz, M., García-Alberola, A., et al. (2013). Analysis of the robustness of spectral indices during ventricular fibrillation. Biomedical Signal Processing and Control 8, 733–739. doi:10.1016/j.bspc.2013.06.013

CrossRef Full Text | Google Scholar

Ricci, R., Pignalberil, C., Disertori, M., Capucci, A., Padeletti, L., Botto, G., et al. (2001). Antitachycardia pacing therapy to treat spontaneous atrial tachyarrhythmias: the 7250 dual defibrillator italian registry. European Heart Journal Supplements 3, P25–P32. doi:10.1016/S1520-765X(01)90069-0

CrossRef Full Text | Google Scholar

Sanders, P., Berenfeld, O., Hocini, M., Jaïs, P., Vaidyanathan, R., Hsu, L.-F., et al. (2005). Spectral analysis identifies sites of high-frequency activity maintaining atrial fibrillation in humans. Circulation 112, 789–797. doi:10.1161/CIRCULATIONAHA.104.517011

PubMed Abstract | CrossRef Full Text | Google Scholar

Schüttler, D., Bapat, A., Kääb, S., Lee, K., Tomsits, P., Clauss, S., et al. (2020). Animal models of atrial fibrillation. Circulation Research 127, 91–110. doi:10.1161/CIRCRESAHA.120.316366

PubMed Abstract | CrossRef Full Text | Google Scholar

van Boxtel, G., Short, T., Kienzle, P., Abbott, B., Aguado, J., Annamalai, M., et al. (2021). Gsignal: signal processing. Available online at: https://cran.r-project.org/web/packages/gsignal/gsignal.pdf.

Google Scholar

Viskin, S., Ish-Shalom, M., Koifman, E., Rozovski, U., Zeltser, D., Glick, A., et al. (2003). Ventricular flutter induced during electrophysiologic studies in patients with old myocardial infarction: clinical and electrophysiologic predictors, and prognostic significance. Journal of Cardiovascular Electrophysiology 14, 913–919. doi:10.1046/j.1540-8167.2003.03082.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Weisfeldt, M. L., Sitlani, C. M., Ornato, J. P., Rea, T., Aufderheide, T. P., Davis, D., et al. (2010). Survival after application of automatic external defibrillators before arrival of the emergency medical system: evaluation in the resuscitation outcomes consortium population of 21 million. Journal of the American College of Cardiology 55, 1713–1720. doi:10.1016/j.jacc.2009.11.077

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H., and Bryan, J. (2023). R Packages: Organize, Test, Document, and Share Your Code. 2nd edn. Sebastopol, CA: O’Reilly Media.

Google Scholar

Wickham, H., Averick, M., Bryan, J., Chang, W., McGowan, L. D., François, R., et al. (2019). Welcome to the tidyverse. Journal of Open Source Software 4, 1686. doi:10.21105/joss.01686

CrossRef Full Text | Google Scholar

Zaglia, T., Pianca, N., Borile, G., Da Broi, F., Richter, C., Campione, M., et al. (2015). “Optogenetic determination of the myocardial requirements for extrasystoles by cell type-specific targeting of channelrhodopsin-2,” 112. Proceedings of the National Academy of Sciences, E4495–E4504. doi:10.1073/pnas.1509380112

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuravlev, M. O., Runnova, A. E., Mironov, S. A., Zhuravleva, J. A., and Kiselev, A. R. (2025). Continuous estimation of heart rate variability from electrocardiogram and photoplethysmogram signals with oscillatory wavelet pattern method. Sensors 25, 5455. doi:10.3390/s25175455

PubMed Abstract | CrossRef Full Text | Google Scholar