Sleep Apnea Classification Using the Mean Euler–Poincaré Characteristic and AI Techniques

Abstract

Sleep apnea is a sleep disorder that disrupts breathing during sleep. This study aims to classify sleep apnea using a machine learning approach and a Euler–Poincaré characteristic (EPC) model derived from electrocardiogram (ECG) signals. An ensemble K-nearest neighbors classifier and a feedforward neural network were implemented using the EPC model as inputs. ECG signals were preprocessed with a polynomial-based scheme to reduce noise, and the processed signals were transformed into a non-Gaussian physiological random field (NGPRF) for EPC model extraction from excursion sets. The classifiers were then applied to the EPC model inputs. Using the Apnea-ECG dataset, the proposed method achieved an accuracy of 98.5%, sensitivity of 94.5%, and specificity of 100%. Combining machine learning methods and geometrical features can effectively diagnose sleep apnea from single-lead ECG signals. The EPC model enhances clinical decision-making for evaluating this disease.

1. Introduction

Sleep apnea is a sleep disorder. During sleep, decreased breathing occurs in periods of a few seconds to minutes; this condition is known as sleep apnea. Generally, there is a high prevalence of undiagnosed cases among the population suffering from obstructive sleep apnea (OSA). This prevalence is higher in men than women; nearly one billion people have OSA worldwide [1]. OSA affects people with cardiovascular disease (arrhythmia, heart failure, coronary syndrome) and increases sympathetic activity that compromises the heart [2,3,4]. Also, OSA affects people with comorbidities with hypertension, as mentioned in [5].

OSA can be diagnosed using polysomnography; this method examines sleep and respiration parameters using electroencephalograms (EEGs), electrocardiograms (ECGs), and other measures, e.g., pulse oximetry ($S_P{O}_2$), airflow measurement, and so on  [6,7]. This process is complex and time-consuming because it requires a complete test in a controlled environment to monitor the patient’s sleep and other signals. Also, this diagnosis is unfeasible and extremely expensive for a large population. Some alternative methods were developed using single ECG signals and machine learning techniques to diagnose OSA problems.

Studies, including [8,9,10], have demonstrated the capability of ECG signals to identify respiratory events. Furthermore, research results such as in [11,12], indicate that heart rate variability (HRV) signals can be used to measure and detect apnea events. Atrial fibrillation (AF) is widely related to OSA problems, as mentioned in [13,14], where the authors conducted questionnaires and home sleep studies to check this relation.

Researchers have explored two approaches to detect sleep apnea events in the literature: feature engineering-based techniques and deep learning-based methods. In the first category, the ECG signal is analyzed, and different features are extracted from the P-wave (atrial depolarization), QRS complex (ventricular depolarization), and T-wave (ventricular repolarization), which reflect different cardiac cycle phases. For example, Babaeizadeh et al. presented a cost-effective, noninvasive algorithm for detecting sleep apnea using single-lead ECG. By analyzing heart rate variability and QRS complexes, the algorithm accurately classifies apnea episodes, achieving 84.7% accuracy and demonstrating high potential for sleep medicine applications [15]. Yilmaz explored the utility of using a single-lead ECG to determine sleep stages and detect OSA across 30-s epochs throughout the night, potentially replacing more complex polysomnography systems. Data from 18 subjects, including 10 diagnosed with OSA, revealed promising results: support vector machines (SVMs) and quadratic discriminant analysis (QDA) achieved over 82% accuracy in classifying sleep stages and around 88% in detecting OSA. Using the RR-interval (time interval between two successive R-peaks in an ECG), three features were extracted, resulting in high classification performance. This suggests that single-lead ECG, combined with SVM or QDA, could effectively monitor sleep and OSA. However, integrating additional physiological features could further enhance accuracy [16]. Zarei et al. described an algorithm that detects sleep apnea using ECG features and a combined convolutional neural network (CNN) and long short-term memory (LSTM) recurrent network [17]. Using the same combination of networks, we looked at the work by Bahrami and Forouzanfar [18], which reported accuracy, sensitivity, and specificity of 88.13%, 84.26%, and 92.27%, respectively.

Feature extraction using machine learning (ML) methods was employed to detect OSA events. In the study by Pant et al. [19], wavelet transform was used to extract features from the sub-bands and classify them using various ML methods. In the work of Yang et al., the HRV analysis and ECG-derived respiration (EDR) were used to extract complementary information using a residual group network [20]; their results had an accuracy of 90.3%, a sensitivity of 87.6%, and a specificity of 91.9%. Liu et al. demonstrated a CNN that learned features from the ECG, reaching a sensitivity of 88.2% in classification [21].

Instead of relying on an overall diagnostic analysis of the entire biomedical signal, certain studies focus on determining if a patient is experiencing apnea by analyzing specific indicators. One approach measures the variation in oxygen saturation over a short interval (12 s) [22,23]; other indicators include non-linear methods mentioned in [24]. Maier et al. used Holter-ECG to screen for sleep-disordered breathing (SDB) by comparing it with nocturnal polygraphy (PG) in 50 cardiology patients. The ECG estimates demonstrated strong agreement with PG, exhibiting high specificity (0.96) and good sensitivity (0.77) for detecting significant SDB, promising an effective, cost-efficient method for SDB screening [25].

Some studies have adopted an approach of using fragments of the signal, such as analyzing it minute by minute [26] or using windows of varying sizes [27], to identify occurrences of apnea instead of making an overall diagnosis of OSA, which is generally regarded as the responsibility of physicians.

Various biomedical signals are utilized for the classification and diagnosis of OSA. For instance, the photoplethysmography signal ([28], the $S_p{O}_2$ signal [24,29,30], the oronasal air flow signal [31], the ECG signal [23,26], and so on. These works use techniques that involve RR variability [26,32,33], frequency domain features using wavelet transform [29,34], and machine learning techniques as neural networks [17,29,34,35,36], and so on.

Artificial intelligence models currently play a huge role in medicine, especially artificial neural networks, which have been used for the classification and detection of diseases such as diabetes [37], heart diseases [38], chronic kidney disease [39], and in oncological diagnostics [40]. Due to the efficiency and utility of these tools, we use neural networks in this work.

Our study makes significant contributions to the field of ECG analysis by introducing an innovative algorithm that integrates geometric models as features and inputs for AI models, specifically designed to classify and detect sleep apnea and heart diseases. This algorithm leverages geometric features that closely match the behavior of the ECG, enhancing the accuracy of classifications. We ensure relevance and continuity in ECG research by addressing the challenges posed by limited and imbalanced datasets and aligning our work with historical benchmarks established in the 2000 challenge [41] using the Apnea-ECG database [42]. Overall, this comprehensive approach not only deepens our understanding of ECG signal behavior but also lays the groundwork for future studies involving larger and more balanced datasets, ultimately improving diagnostic capabilities and patient outcomes.

Existing OSA detection methods often rely on handcrafted features and struggle to achieve high accuracy with single-lead ECG signals. To address these limitations, this article introduces a novel algorithm that uses random algebraic polynomials (RAP) for signal smoothing and creates a non-Gaussian physiological random field (NGPRF) to capture intricate, non-Gaussian ECG features. This innovative framework models critical ECG peaks (P-waves, QRS-complex, T-waves) and extracts their geometrical characteristics as Euler–Poincaré characteristics (EPCs). By capturing unique geometrical forms in the ECG signal, the EPCs provide robust inputs to our classification system. Unlike conventional approaches, we integrate these advanced geometric descriptors with machine learning classifiers, including ensemble K-nearest neighbors (EKNN) and feedforward multi-layer neural networks (FNNs), to significantly improve OSA detection accuracy. This novel combination of RAP-based signal processing and advanced classification methods bridges existing gaps, making single-lead ECG-based OSA detection more accurate and potentially adaptable for real-time and diverse clinical applications.

This article is organized as follows: In Section 2, we describe the mathematics involved. In Section 3, we describe the data used. In Section 4, we discuss the methodology applied to the cases. Then, in Section 5, we present the results and discussion. Finally, a conclusion and future works are provided in Section 6.

2. Mathematical Framework

This section can help one understand the main steps of the methodology mentioned in [43].

2.1. Polynomization

The first step of the method involves pre-processing to smooth the signal. In this step, we model the ECG signals using RAP to reduce noise and construct the NGPRF. RAP is denoted by $\phi (\mathit{X},\xi )$, for an ECG cycle j from 1 to N. The elemental waves in the ECG are P-waves, QRS complex, and T-waves. In this case, each cycle is divided by a threshold $R\left(j\right)$, resulting in two intervals: $I_1=\left[1R\left(j\right)\right]$ and $I_2=\left[R\mathit{E}\left(j\right)\right]$, with $\mathit{E}\left(\mathit{j}\right)$ representing the final index in each cycle. Additionally, we define the RAP of order D representing the P-waves, as follows:

$${\phi }_B(\mathcal{S},\overline{\mathsf{\Omega }},X\left(j\right))={\sum }_{k=0}^D{a}_k(j,\overline{\mathsf{\Omega }})X{\left(j\right)}^k,$$

(1)

where $\mathcal{S}=\left(\mathsf{\Omega },\mathcal{F},\mathcal{P}\right)$ represents a complete probability space, $\mathsf{\Omega }$ denotes the sample space, $\mathcal{F}$ denotes the $\sigma$-algebra of events, and $\mathcal{P}$ denotes the probability measure on $\mathcal{F}$. Moreover, $a_k(j,\overline{\omega })$ denotes a sequence of independent and identically distributed random variables for each cycle. The selected interval includes both the R-peak and the P-wave, represented as $X\left(j\right)=[1,R(j\left)\right]$, where $R\left(j\right)$ indicates the index of the R-peak.

2.2. Building Non-Gaussian Physiological Random Field

This framework defines a two-dimensional NGPRF, $\mathsf{\Phi }=(\mathit{X},\xi )$, within a complete probability space $\mathcal{S}$. It models physiological signals across cardiac cycles, mapping $\theta :\xi \to {\mathcal{R}}^{\mathcal{D}}$. Here, $\mathit{X}\in \mathcal{D}\subset {\mathcal{R}}^2$, and $\theta$ denotes a two-dimensional NGPRF defined over N cardiac cycles as follows:

$$\mathsf{\Phi }(X_1,X_2,\overline{\xi })=\sum _{c=1}^N{\phi }_B(\mathcal{S},\overline{\xi },X\left(j\right))\delta (j-X_2)$$

(2)

with ${\phi }_B(\mathcal{S},\overline{\xi },X\left(j\right))$ as the random polynomial that represents the ECG behavior, and $\delta$ is the Dirac’s distribution. According to Equation (2), polynomials ${\phi }_B(\mathcal{S},\overline{\xi },X\left(j\right))$ are used to decompose each ECG signal at a fixed point into an NGPRF in the coordinate $X_2$.

2.3. Excursion Set

In the random field theory, the excursion set is the main concept. Consider $\theta (\mathcal{X},\xi )$ as NGPRF, where $\mathcal{X}\in \mathcal{D}\subset {\mathcal{R}}^2$ is defined within a set $\mathcal{D}$. The excursion set is a geometric entity defined as the collection of points where the field exceeds a certain threshold value. Specifically, for a given threshold $\lambda$, the excursion set is given by the following:

$${\mathcal{A}}_{\lambda }(\theta ,\mathcal{D})=\{\mathcal{X}\in \mathcal{D},\theta (\mathcal{X},\xi )\ge \lambda \}$$

(3)

The set ${\mathcal{A}}_{\lambda }$ of $\theta (\mathcal{X},\xi )$ above a threshold $\lambda$ consists of points in $\mathcal{D}\subset {\mathcal{R}}^2$, where $\theta (\mathcal{X},\xi )$ exceeds $\lambda$. Stratified manifolds $\mathcal{D}\in {\mathcal{R}}^2$ can be partitioned into disjoint unions of $k-$dimensional manifolds, expressed as $\mathcal{D}={\cup }_{k=0}^{dim\mathcal{D}}{\partial }_k\mathcal{D}$, where each stratum ${\partial }_k\mathcal{D}$ represents several $k-$dimensional manifolds. The non-Gaussian properties of $\theta$ necessitate the consideration of generalized random fields, expressed as follows:

$$\theta (\mathcal{X},\xi )=F\left({\theta }_{\ast }(\mathcal{X},\xi )\right)=F({\theta }_{\ast }^1(\mathcal{X},\xi ),\dots ,{\theta }_{\ast }^m(\mathcal{X},\xi ))$$

(4)

where $\theta (\mathcal{X},\xi )$ represents a collection of independent and identically distributed (i.i.d.) Gaussian random fields (GRFs), each defined over a topological space $\mathcal{D}$. The function $F:{\mathcal{R}}^m\to \mathcal{R}$ is a smooth, piecewise twice continuously differentiable function with appropriate boundary conditions. The excursion set of a real-valued non-Gaussian field $\theta =F\circ {\theta }_{\ast }$ above a threshold $\lambda$ is equivalent to the excursion set for a vector-valued Gaussian field ${\theta }_{\ast }$ in $F^{-1}[\lambda ,\infty )$, and is given by the following:

$${\mathcal{A}}_{\lambda }(\theta ,\mathcal{D})={\mathcal{A}}_{\lambda }(F\circ {\theta }_{\ast },\mathcal{D})=\{\mathcal{X}\in \mathcal{D},(F\circ {\theta }_{\ast })(\mathcal{X},\xi \ge \lambda \}$$

(5)

In our context, let ${\mathbf{S}}_S=(s_1,\dots ,s_{S^M})$ denote the set of signatures, where each $s_k$ represents a signed unit weight associated with the sub-Gaussian random field ${\theta }_{s_k}(\mathcal{X},\xi )$. Furthermore, we define $S_M={\mathsf{\Sigma }}_{k=1}^M{(-1)}^k{s}_k$. We shall characterize the NGPRF, defined by Equation (2), as follows:

$$\theta (\mathcal{X},\xi ,S_M)=(s_1{\theta }_{S_1}(\mathcal{X},\xi ),\dots ,s_{S_M}{\theta }_{S_{S_M}}(\mathcal{X},\xi ),$$

(6)

where $s_{S_M}{\theta }_{s_{S_M}}(\mathcal{X},\xi )$ represents a set of i.i.d.-signed GRFs with $\mathcal{X}\in {\mathcal{R}}^2$, $\xi \in \mathsf{\Omega }$, and $s_{2n}=+1$, $s_{2n+1}=-1$.

2.4. Euler–Poincaré Characteristic

Let $\varphi \left(A_{\lambda }(\theta ,\mathcal{D})\right)$ represent the EPC of the set $A_{\lambda }(\theta ,\mathcal{D})$. For multiple trials, if the expected value $\mathcal{E}\left\{\varphi \left(A_{\lambda }(\theta ,\mathcal{D})\right)\right\}$ is computable, the random field theory imposes specific regularity conditions on $\theta (\mathit{X},\xi )$. These conditions ensure that the realizations of $\theta$ are smooth and that the boundary $\partial \mathcal{D}$ is smooth as well. Thus, $\mathcal{D}$ is regarded as a regular $C^2$ domain within a compact subset of ${\mathcal{R}}^2$, bounded by a smooth, one-dimensional $C^2$ manifold $\partial \mathcal{D}$.

Consider $\theta (\mathcal{X},\xi )$, where $\mathcal{X}=[{\chi }_1,{\chi }_2]\subset {\mathbb{R}}^2$, as a stationary, non-isotropic random field. Define ${\dot{\theta }}_i(\mathcal{X},\xi )={\displaystyle \frac{\partial \theta (\mathcal{X},\xi )}{\partial {\chi }_i}}$ and ${\ddot{\theta }}_{ij}(\mathcal{X},\xi )={\displaystyle \frac{{\partial }^2\theta (\mathcal{X},\xi )}{\partial {\chi }_i\partial {\chi }_j}}$, where $i,j=1,2$.

The moduli of continuity of ${\dot{\theta }}_j(\mathcal{X},\xi )$, and ${\ddot{\theta }}_{jk}(\mathcal{X},\xi )$ inside $\mathcal{D}$, are given by the following:

$${\xi }_j\left(h\right)=\underset{\parallel \mathcal{X}-\mathcal{Y}\parallel <h,\xi }{sup}\left|{\dot{\theta }}_j(\mathcal{X},\xi )-{\dot{\theta }}_j(\mathcal{Y},\xi )\right|,$$

and

$${\xi }_{jk}\left(h\right)=\underset{\parallel \mathcal{X}-\mathcal{Y}\parallel <h,\xi }{sup}\left|{\ddot{\theta }}_{jk}(\mathcal{X},\xi )-{\ddot{\theta }}_{jk}(\mathcal{Y},\xi )\right|,$$

respectively.

To ensure that the realizations of $\theta (\mathcal{X},\xi )$ are sufficiently smooth, the following conditions must be met:

• C1: The maximum of ${\xi }_j\left(h\right)$ and ${\xi }_{jk}\left(h\right)$ should be asymptotically negligible, specifically $o\left(h^N\right)$, as $h\to 0$.
• C2: The Hessian $\ddot{\mathit{\theta }}$ of ${\dot{\theta }}_{jk}$ must have a finite variance conditional on $\theta$, where $\dot{\theta }$ is the gradient of ${\dot{\theta }}_j$.
• C3: The densities of $\theta$ and $\dot{\theta }$ must be uniformly bounded for all $\mathcal{X}\in \mathcal{D}$.

At a boundary point $\mathcal{X}$ on $\partial \mathcal{D}$:

• ${\dot{\theta }}_{\perp }$ is the gradient of $\theta$ in the direction perpendicular to $\partial \mathcal{D}$.
• ${\dot{\theta }}_{\mathcal{D}}$ is the gradient in the tangent direction within the plane of $\partial \mathcal{D}$.
• ${\ddot{\theta }}_{\mathcal{D}}$ denotes the $1\times 1$-Hessian matrix in the tangent plane.
• r is the $1\times 1$ matrix representing the internal curvature of $\partial \mathcal{D}$.

Using the sign function, denoted as $sign$, and following Knuth’s notation [44], which assigns a value of 1 to a true expression in parentheses, and 0 if false, the EPC is then given by the following:

$$\begin{array}{cc}\hfill \mathsf{\Psi }\left({\mathcal{A}}_{\lambda }(\theta ,\mathcal{D})\right)& ={\mathsf{\Sigma }}_{\mathcal{X}\in \mathcal{D}}(\theta \ge \lambda )(\dot{\mathit{\theta }}=0)sign[det(-\ddot{\mathit{\theta }})]\hfill \\ & +{\mathsf{\Sigma }}_{\mathcal{X}\in \mathcal{D}}(\theta \ge \lambda )({\dot{\mathit{\theta }}}_{\mathcal{D}}=0)({\dot{\theta }}_{\perp }\le 0)\times \mathrm{sign}[det(-{\ddot{\mathit{\theta }}}_{\mathcal{D}})-\dot{\theta }\perp r)],\hfill \end{array}$$

(7)

with a probability of one. In ${\mathcal{R}}^2$, integral geometry characterizes $\phi \left({\mathcal{A}}_{\lambda }\left(\theta ,\mathcal{D}\right)\right)$ as the difference between the count of connected components and the number of voids in ${\mathcal{A}}_{\lambda }\left(\theta ,\mathcal{D}\right)$. Moreover, the expected value of $\phi \left({\mathcal{A}}_{\lambda }\left(\theta ,\mathcal{D}\right)\right)$ over multiple realizations is as follows:

$$\begin{array}{cc}\hfill \mathbb{E}\left\{\mathsf{\Psi }\right({\mathcal{A}}_{\lambda }& (\theta ,\mathcal{D})\left)\right\}={\int }_{\mathcal{D}}\mathbb{E}\{(\theta \ge \lambda )det(-\ddot{\mathit{\theta }})|\dot{\mathit{\theta }}=0\}\beta \left(0\right)d\mathit{X}\hfill \\ & +{\int }_{\partial \mathcal{D}}\mathbb{E}\{(\theta \ge \lambda )({\dot{\theta }}_{\perp }<0)det(-{\ddot{\mathit{\theta }}}_{\mathcal{D}}-\dot{\mathit{\theta }}\perp r)|{\dot{\theta }}_{\mathcal{D}}=0\}\times {\beta }_{\mathcal{D}}\left(0\right)d\mathcal{X},\hfill \end{array}$$

(8)

where $\beta (.)$ is the density of $\dot{\theta }$ and ${\beta }_{\mathcal{D}}(.)$ is the density of ${\dot{\theta }}_{\mathcal{D}}$. It is shown that, under the condition that the boundary of $\mathcal{D}$ consists of a finite number of piecewise filtered segments, the expected value of the excursion set of a mean-zero, unit-variance random field is given by the following:

$$\mathbb{E}\left\{\mathsf{\Psi }\left({\mathcal{A}}_{\lambda }\left(\theta ,\mathcal{D}\right)\right)\right\}=\sum _{j=0}^N{\rho }_j\left(\lambda \right){\mathcal{L}}_j\left(\mathcal{D}\right)$$

(9)

where ${\rho }_j\left(\lambda \right)={\left(2\pi \right)}^{-(j+1)/2}{H}_{j-1}\left(\lambda \right)exp(-{\lambda }^2/2)$ is the intensity of the EPC per unit volume, $H_j$ denotes the $j-th$ Hermite polynomial, and ${\mathcal{L}}_j\left(\mathcal{D}\right)$ represents the Lipschitz–Killing curvatures.

3. Data Used

The Apnea-ECG database [42] was used in this study. It contains 70 records, divided into two groups: a training set of 35 records, including 20 severe cases, 5 mild cases, and 10 non-apnea cases, and a validation set consisting of the remaining 35 records. All ECG signals in the dataset are annotated with QRS complexes. The data are derived from a population aged between 27 and 60 years, with all records digitized at 100 samples per second. The length of the records ranges from 7 to 10 h.

4. Methodology

Figure 1 represents the proposed methodology used for classifying sleep apnea. The methodology consists of the following five steps:

• Pre-processing: As part of this step, the ECG signals were filtered using random polynomials. This process is referred to as polynomization.
• Random field transformation: The filtered signal was used to build the NGPRF, collapsing c number of cycles into a geometric structure.
• Geometrical property: Each NGPRF was sectioned at several levels, denoted by $\lambda$, known as excursion sets.
• Feature extraction: For each excursion set, a value was calculated that represented the difference between the number of connected components and the number of holes, capturing key aspects of the set’s geometric structure. This feature is referred to as the Euler–Poincaré characteristic.
• Classifier selection: An EKNN model is proposed to classify OSA from the EPC models. The training dataset was used for learning, and the test dataset was used to validate the model. Also, a feedforward multi-layer neural network is proposed for the binary classification.

4.1. Pre-Processing

The pre-processing stage is referred to as ‘polynomization’, and was dealt with as follows. First, the ECG cycle was divided into two interval parts, as our focus was on the part containing the P-waves and R-peaks, only one interval was selected. The function for this process is presented in (1).

We identified the R-peaks using the Pan–Tompkins algorithm [45] to achieve this objective. Afterward, we split the ECG cycle from R-peak to R-peak using the middle points, $m_k$, and a fixed threshold. Each ECG cycle preserved the main waves, i.e., P, Q, R, S, and T. Every cycle $N_j$ was segmented two sections, $P_1\left(C\right)=\left[m_k{R}_{k+1}\right],$ and $P_2\left(C\right)=\left[R_{k+1}{m}_{k+1}\right]$. Upon completion of the process, the result was in the form of polynomials that correspond to the appropriate order D. In these cases, we only used the $P_1\left(C\right)$ polynomials.

We assessed the normalized root mean square error (NRMSE) using Equation (10) to identify the optimal order of the random polynomials that best approximated the data.

$$FIT=100\left(1-{\displaystyle \frac{∥P_1\left(C\right)-{\sum }_{k=0}^{D_1\left(C\right)}{\widehat{a}}_k(k,\overline{\omega })X_j{\left(k\right)}^k∥}{∥P_1\left(C\right)-{\overline{P}}_1\left(C\right)∥}}\right)$$

(10)

where $P_1\left(C\right)$ represents the section of the ECG containing the P-waves and R-peaks, ${\overline{P}}_1\left(C\right)$ is the mean of each interval $P_1$, and the inverse of the Vandermonde matrix is used to calculate the vector of estimated coefficients ${\widehat{a}}_k$ [46]. Our goal is to find the maximum FIT based on NRMSE. To do this, we evaluate a group of samples (100,000) from each of the 35 patients for the training set using polynomial orders ranging from 15 to 25 ($D_1\left(c\right)=15,\dots ,25$). To calculate the order of polynomials for the sections of the ECG from the R-peaks to the T-waves, we used Equation (10), and instead of $P_1\left(C\right)$, we used $P_2\left(C\right)$. Also, we used the integral absolute error (IAE) criterion. The results for $P_1\left(C\right)$ and $P_2\left(C\right)$ are shown in Figure 2a,b with both criteria.

The optimal polynomial order is 15 according to the FIT criterion and 21 based on the IAE criterion. It is later shown that the order of the polynomials significantly influences the final results, particularly during the classification process. Algorithm 1 shows the entire procedure used for selecting the orders and obtaining the polynomials $P_1\left(C\right)$ and $P_2\left(C\right)$. However, we only use the P-wave polynomials $P_1\left(C\right)$ in the following steps.

Algorithm 1 Algorithm for pre-processing the ECG signals

Input:  ECG Output:  Polynomials $P_1\left(C\right)$ and/or $P_2\left(C\right)$      Initialization;  1: Search the locations $M_j$, $j=1,\dots ,N_c$;  2: Locate the middle points $m_j$;  3: Find the best order for $P_1\left(C\right)$ and $P_2\left(C\right)$ using Equation (10);  4: Estimate the polynomials $I_1\left(c\right)$, $c\in [m_j,M_{j+1}]$, and $I_2\left(c\right)$, $c\in [M_{j+1},m_{j+1}]$;  5: Assess the results;

4.2. Transforming into a Random Field and Geometrical Approach

Once the ECG signal is filtered using the polynomization method, we combine the results of each cardiac cycle to form the random field (see Figure 3). The NGPRF can be expressed as Equation (2).

The excursion set, defined in Equation (3), utilizes a threshold level $\lambda$ determined from the range [−1, 1]. This range is used because we normalize the signal in the pre-treatment process.

To ensure homogeneity across the random fields, any empty spaces should be filled with zeros. Figure 4 and Figure 5 show the behaviors of the excursion set at threshold level $\lambda =0$, where we note the R-peaks and the P-waves of each type of patient. We developed Algorithm 2, which measures the vectors with the P-waves and R-peaks to build the non-Gaussian random field. The algorithm first determines the maximum length of the cycles. If the length of any cycle is less than this maximum, the algorithm fills the empty spaces with zeros to ensure homogeneity.

Algorithm 2 Process to construct the non-Gaussian random physiological field

Input:  RAPs $P_1\left(C\right)$, empty vector $A=\left[\right]$ Output:  $\mathsf{\Phi }(\mathit{X},\overline{\mathsf{\Omega }},S^M)$      Initialization;  1: for $i=1to$$M_c$ do  2:    $M\left(i\right)=length\left(P_1{\left(C\right)}_i\right)$  3: end for  4: $Va{l}_M=max\left\{Matrix\right\}|Va{l}_M\in \mathcal{R}$  5: if ($Matri{x}_i<Va{l}_M$) then  6:    $\mathsf{\Phi }([i,1:M],\overline{\mathsf{\Omega }})=[zeros,P_1{\left(C\right)}_k]$  7: else  8:    $\mathsf{\Phi }([i,1:M],\overline{\mathsf{\Omega }})=P_1{\left(C\right)}_k$  9: end if

4.3. Extracting the Feature

In this work, feature extraction relies on the excursion set to calculate the Euler–Poincaré characteristic (EPC), represented as $\phi \left({\mathcal{A}}_{\lambda }(\mathsf{\Phi },\mathcal{D})\right)$. The approach involves transforming the excursion set into a binary image, where each value is set to 1 if it exceeds the threshold level $\lambda$, and 0 otherwise. To compute the EPC values, we applied Gray’s algorithm [47]. The EPC is a scalar quantity that represents the difference between the number of connected components and the number of holes within those components [48].

After that, the EPC is obtained for all apnea and non-apnea patients in the training set. Figure 6 shows the EPC for each apnea patient; the blue lines depict the results obtained with 15th-order polynomials, while the red lines indicate the behaviors associated with 21st-order polynomials. Figure 7 shows the borderline apnea patients (patients 21–25) and non-apnea patients (26–35), where we have the results of both polynomials, i.e., the blue lines denote 15th-order and the red lines denote 21st-order.

4.4. Classification

The last step of the methodology involves searching for the best classifier. To predict the type of apnea, we use a matrix $\mathcal{X}$ with 102 columns; columns 1 to 101 contain the stacked EPC values, and the last column represents the target value, where sleep apnea is denoted as 0 and non-apnea as 1. The matrix described is utilized for training the classifier.

$$\left[\begin{array}{ccccc}{\mathcal{Z}}_1^1& {\mathcal{Z}}_2^1& \cdots & {\mathcal{Z}}_n^1& Y_1\\ {\mathcal{Z}}_1^2& {\mathcal{Z}}_2^2& \cdots & {\mathcal{Z}}_n^2& Y_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ {\mathcal{Z}}_1^m& {\mathcal{Z}}_2^m& \cdots & {\mathcal{Z}}_n^m& Y_m\end{array}\right]$$

where ${\mathcal{Z}}_j^k$ denotes the predictors with $j=[1,101]$, and $Y_k$ denotes the class of each patient.

In our study, we utilized Classification Learner, a Matlab toolbox, to model several classifiers. We chose several models to provide a baseline comparison, varying complexity and interoperability. Our objective was to evaluate simple and sophisticated models to observe the complexity of the OSA classification. The results are shown in

Table 1. Comparing machine learning techniques to classify sleep apnea using EPC models as inputs.

Type of Classifier	Accuracy	Sensitivity	Specificity	F1 Score	MCC
Linear SVM	62.86%	70%	20%	76.36%	−7.75%
Decision Tree	65.71%	78.26%	41.67%	75%	20.94%
KNN	74.29%	78.57%	57.14%	83.02%	31.62%
Logistic Regression	77.14%	86.96%	58.33%	83.33%	47.59%
Ensemble KNN	80%	82.14%	71.43%	86.79%	47.43%

, applying EPC models derived from 15th-order polynomials as inputs for the classifiers. The computations were conducted on an AMD Ryzen 5 processor with 16 GB of RAM, utilizing the parallel pool option. We used k-fold cross-validation with $k=5$ in all cases.

With the results shown previously, we applied an ensemble of K-nearest neighbors (EKNN); this approach can help to improve the prediction with non-informative features [49]. In Figure 6 and Figure 7, we observed all the EPCs from all the patients; analyzing these data presents a challenge, so we employed classifiers to reveal information that aids in their classification.

Figure 8 presents a diagram of the ensemble KNN classifier, where the EPC models from both healthy and unhealthy patients serve as the training data.

We conducted a study where we selected 30 learners, each representing a different subspace. For each learner, we randomly chose 65 variables to predict their behavior, resulting in 35 observations for both the training and test phases. In the next section, we will discuss the results and findings of this model in detail.

Table 2. Parameters used in the ensemble KNN model.

Parameters	Values
Number of Neighbors (k)	10
Weighting Method	‘distance’
Distance Metric	‘euclidean’
Algorithm for Finding Neighbors	‘auto’
Number of Estimators	30
Ponderation of Models	Uniform

shows the parameters used to train the EKNN model using 101 features as inputs. Remember that from the Apnea-ECG database, there are two sets, one called training data and the other called test data, each with 35 records unbalanced in the classes.

4.5. Neural Network Classifier

This work uses a feedforward multi-layer neural network to classify apnea or non-apnea patients. In this architecture, data move in one direction from the input layer through one or more hidden layers to the output layer. FNNs are commonly used for classification tasks, such as determining whether a patient has sleep apnea. The network is fully connected, meaning each neuron in one layer connects to every neuron in the next layer. The model is trained by adjusting the weights between these connections to minimize the difference between the predicted and actual labels using a process called backpropagation. In binary classification, the output typically comes from a single neuron using a sigmoid activation function to produce a probability.

The main equation governing the computation in a neural network is the weighted sum of the inputs at each neuron, followed by an activation function. For neuron j, the pre-activation output $z_j$ is given by the following:

$$z_j=\sum _{i=0}^n{\omega }_{ji}{x}_i+b_j$$

(11)

where ${\omega }_{ji}$ denotes the weights, $x_i$ denotes the input values, and $b_j$ denotes the bias. The weighted sum $z_j$ is then passed through an activation function $f\left(z_j\right)$ to introduce non-linearity. This nonlinearity is essential for the network to model complex patterns in the data. The final output from the network depends on both the network architecture and the activation function choice.

Figure 9 shows the neural network used in this work to classify apnea, where an input layer, two hidden layers, and an output layer were implemented. A Swish activation function was implemented in the input layer, a ReLU activation function was implemented in the first hidden layer, a tanh activation function was used in the second hidden layer, and a softmax was used in the output layer.

Table 3. Summary of feedforward multi-layer neural network parameters.

Parameter	Description	Value
Input Layer	Number of input features	101 (from EPC model)
Hidden Layers	Number of hidden layers	2
Neurons	Number of neurons per layer	100 (Layer 1), 64 (Layer 2), 64 (Layer 3)
Output Layer	Number of neurons in output layer	1
Activation Functions	Activation function per layer	Swish (Layer 1), ReLU (Layer 2), Tanh (Layer 3), Softmax (Output Layer)
Optimizer	Optimization algorithm	Adam
Epochs	Number of epochs	100
Classes	Number of output classes	2 (binary classification)

summarizes the values of the parameters of the FNN builder in Python using the TensorFlow library.

Below are the results obtained using two types of polynomial orders to smooth the signal. They are also compared with other models.

5. Results and Discussion

A machine learning model was developed to detect cases of sleep apnea using an ensemble K-nearest neighbors (KNN) classifier. The first phase of the study involved training the model using a dataset that included 20 cases of severe sleep apnea (1 h of activity), 5 mild cases (<1 h), and 10 cases with non-apnea. The results of this phase are presented in the form of a confusion matrix, as shown in Figure 10. The model’s accuracy was found to be 80%. Also, we found a sensitivity of 82% and a specificity of 71%. To achieve these results, we used the following formulas using the results of the confusion matrix:

$$Accuracy=(TP+TN)/(P+N)=(23+5)/\left(35\right)\times (100\%)=80\%$$

$$Sensitivity=\left(TP\right)/(TP+FN)=\left(23\right)/(23+5)\times (100\%)=82\%$$

$$Specificity=\left(TN\right)/(TN+FP)=\left(5\right)(5+2)\times (100\%)=71\%$$

where P and N represent positive and negative cases, respectively; $TP$ and $TN$ denote true positives and true negatives; and $FP$ and $FN$ represent false positives and false negatives.

Figure 10 shows the total correct predictions of the trained model, where it achieved 92% accuracy in apnea cases and 50% accuracy in non-apnea cases. In this work, we used 15th- and 21st-order polynomials to soften the signal, but the model achieved the same results obtained with the training data; Figure 10 expresses the results for the two polynomial orders in this phase.

In the testing phase, we used the validation data to obtain the results presented in Figure 11 and Figure 12; one set of results uses a 15th-order polynomial and the other uses a 21st-order polynomial. Figure 11 shows that patients with apnea were correctly classified 100% of the time, while those without apnea were correctly classified 88% of the time. In addition, the accuracy, sensitivity, and specificity were calculated with the above-mentioned formulas, yielding results of 97%, 96%, and 100%, respectively.

Figure 12 shows the results using the 21st-order polynomial, where the patients with apnea were correctly classified 100% of the time, and those without apnea were correctly classified 75% of the time. Additionally, the accuracy, sensitivity, and specificity were calculated using the formulas mentioned above, yielding results of 94%, 93%, and 100%, respectively.

A new method for classifying sleep apnea is proposed, which involves performing a geometric transformation and extracting information through the Euler–Poincaré characteristic model from Holter ECG and one lead (lead $V_3$). Using a 15th-order polynomial in the pre-processing phase, the model achieves an average accuracy of 89%, a sensitivity of 89%, and a specificity of 86%.

The results obtained with the feedforward neural network are presented for the training phase in Figure 13 and for the test phase in Figure 14. The summary of these results is shown in

Table 4. Table presenting the accuracy, sensitivity, specificity, and F1 score results of the FNN classifier.

Phase	Accuracy	Sensitivity	Specificity	F1 Score
Training	100%	96%	100%	98%
Test	97%	93%	100%	95%

.

The FNN has a better response than the EKNN. We used a 15th-order polynomial EPC model to obtain the FNN classifier model because it had better results compared to the 21st-order polynomial. The next step is to compare the two approaches with several cases in the literature.

Table 5. Comparison between our approach and the approaches used in the literature.

Literature Work	Efficiency	Analyzed Signal	Proposed Techniques	Classification
Mcnames et al. [26]	Accu = 92.6%	ECG	HRV in RR, T and S ECG pulse Energy(WT)	Threshold (5 min window)
Raymond et al. [34]	Accu = 92%	ECG	EDR signal RR signal	Shared mixture classifier (spectral features)
De Chazal et al. [32]	Accu = 89.4% Sens = 84.1% Spec = 90.0%	ECG	RR variability R-wave amplitude (PSD)	Linear and quadratic discriminants (spectral features)
Lee et al. [29]	Accu = 88% Sens = 98% Spec = 92%	$S_p{O}_2$	WT ADA, DDA, NA	Transform coefficients Threshold
Maier et al. [33]	Accu = 89.8% Sens = 81.3% Spec = 82.8%	ECG	RR series MAV series	Threshold (120 ms. window)
Corthout et al. [50]	Accu = 90% Sens = 84% Spec = 93%	ECG	EDM+HT, EDM+RAS WA	Linear Discriminant classifier (feature set)
Burgos et al. [30]	Accu = 93.03% Sens = 92.35% Spec = 93.52%	$S_p{O}_2$	Desaturation indexes	Bagging with ADTree
Ye et al. [36]	Accu = 99.22% Sens = 99.25% Spec = 99.02%	ECG	R-R intervals	FENet
Rajesh et al. [51]	Accu = 89% Sens = 86% Spec = 92%	ECG	Power Spectrum, Discrete Wavelet Transform, PSO	Random forest classifier
Li et al. [35]	Accu = 97.8% Sens = 98.6% Spec = 93.9%	ECG $Sp{O}_2$ BMI	Feature extraction	Multi-layer FNN
Zarei et al. [17]	Accu = 97.21% Sens = 94.41% Spec = 98.94%	ECG	LSTM-CNN	AHI index
EKNN approach	Accu = 89% Sens = 89% Spec = 86%	ECG	Random Field Euler–Poincaré characteristic	EKNN classifier
FNN approach	Accu = 98.5% Sens = 94.5% Spec = 100%	ECG	Random Field Euler–Poincaré characteristic	Multi-layer FNN

compares our approach with various cases from the 2000 challenge to modern techniques that employ deep learning machines.

The heart rate variability (HRV) analysis has been employed in the literature. For example, [26] used domain time features of the HRV, achieving a performance of 92.6%. In [34], the authors analyzed the frequency domain features from the T-wave, using the wavelet transform, and reached 92.3% accuracy. In [32], the spectral analysis of HRV features achieved an accuracy of 89.36%.

Other works used oxygen saturation ($S_p{O}_2$) data to classify OSA. In [29], three different algorithms using the wavelet transform were compared, and the results yielded an accuracy of 88%.

The decision methods can vary by work. The most used classifiers include random forest, support vector machines, linear or quadratic regression, neural networks, and KNN. In [50], a linear discriminant using frequency domain features extracted from the Hilbert transform achieved 90% accuracy. In [30], a decision tree combined with blood oxygen saturation levels classified the OSA and achieved an accuracy of 93.03%, a sensitivity of 92.35%, and a specificity of 93.52%. In [36], a neural network was used to classify RR interval data and multi-frequency filters, achieving a specificity of 99%.

In [51], a random forest classifier used features from the power spectrum and discrete wavelet transform. The best features were selected as inputs with the aid of the particle swarm optimization method. The reported accuracy was 99.22%.

The authors of [35] and [17] used deep neural networks. A multi-layer feed-forward neural network (FNN) using information from the ECG, $S_p{O}_2$, and BMI yielded an accuracy of 97.8% [35]. In [17], the authors proposed an algorithm that used deep learning layers (convolutional neural networks and long short-term memory) to extract features and classify the OSA. It achieved an accuracy of 97.21%, a sensitivity of 94.41%, and a specificity of 98.94% on the Apnea-ECG database.

The related works and each study’s performance measures, signals used, techniques, and decision methods are summarized in Table 5. Our approach stands out in terms of performance measures (sensitivity, specificity, and accuracy); we also employed a geometrical approach, which is unique among the studies reviewed. We used a different multi-layer FNN and the EPC model as features, obtaining a similar performance to Li et al. [35].

Appendix A outlines our application of the Synthetic Minority Oversampling Technique (SMOTE) to address dataset imbalance in the database used. This enhancement increases the robustness of detecting sleep apnea and heart disease from ECG data.

6. Conclusions and Future Works

This study introduced a novel method for detecting obstructive sleep apnea (OSA) using single-lead ECG signals by leveraging the Euler–Poincaré characteristic (EPC) to represent ECG peaks geometrically. The RAP-based methodology for smoothing the ECG signal and extracting EPC provided a rich feature set that captured the intricate patterns of ECG peaks. We then applied an ensemble K-nearest neighbors (EKNN) classifier to these EPC features, outperforming traditional methods by combining multiple KNN models for improved accuracy. This approach highlights the effectiveness of using geometric representations and advanced ensemble techniques to enhance OSA detection from ECG data.

In conclusion, using the EKNN as a classifier during the training phase, we achieved an accuracy of 80%, with a sensitivity of 82% and a specificity of 71%. In the test phase, we attained an accuracy of 97%, a sensitivity of 96%, and a specificity of 100%.

For the FNN, as a binary classifier, we recorded a perfect accuracy of 100%, with a sensitivity of 96% and a specificity of 100% in the training phase. In the test phase, we achieved an accuracy of 97%, a sensitivity of 93%, and a specificity of 100%.

In future works, the NGRF will be built using smooth signals from the P-waves, QRS-peaks, and T-waves to obtain a model that represents the behavior of the full ECG instead of only one section (P-R section).

Our framework works with the EPC as a direct input to the classifier. In future works, we will include different features in various domains (time, frequency, geometric, and statistical) and larger datasets to classify various diseases.