Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery

Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery | medRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-P4HH5NV'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery View ORCID Profile Hau-Tieng Wu , Thomas M Tolbert , David M Rapoport doi: https://doi.org/10.1101/2025.04.24.25326393 Hau-Tieng Wu 1 Courant Institute of Mathematical Sciences, New York University , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Hau-Tieng Wu For correspondence: hauwu{at}cims.nyu.edu Thomas M Tolbert 2 Division of Pulmonary, Critical Care and Sleep Medicine, Icahn School of Medicine at Mount Sinai , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site David M Rapoport 2 Division of Pulmonary, Critical Care and Sleep Medicine, Icahn School of Medicine at Mount Sinai , New York, NY, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abstract Full Text Info/History Metrics Data/Code Preview PDF Abstract Objective Cardiogenic oscillations (CO) in airflow signals contain valuable physiological information. However, accurately isolating CO from airflow signals, particularly in individuals with sleep apnea, remains a challenging signal processing problem. Method We introduce the Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA), a novel approach for extracting CO from airflow signals while simultaneously recovering a CO-free, noise-free airflow signal, referred to as diaphragm-driven airflow (DDairflow). The algorithm’s performance is quantitatively evaluated using both a semi-real simulated database and real-world data with benchmark comparisons to existing methods, including the bandpass filter (BPF) and Savitzky-Golay smoothing filters (SGF). Result For the semi-real database, OSADA significantly outperforms BPF and SGF across multiple performance indices, including the normalized root mean square error (NRMSE) for CO and DDairflow recovery, as well as spectral energy indices of CO. For real-world data, OSADA also achieves superior performance in the data-driven spectral energy index of CO. Conclusion OSADA is the first algorithm specifically designed for CO recovery from single-channel airflow signals, without relying on additional channels, and is supported by theoretical foundations. Quantitative results suggest robust performance for both CO extraction and DDairflow recovery. 1. Introduction The ventilatory airflow signal is a rich source of valuable physiologic information in many healthcare settings, including polysomnography in sleep centers, ventilation monitoring in operating rooms, and mobile devices in the homecare environment. A common incidental finding in airflow signals is the presence of cardiogenic oscillations (CO) [ 1 - 3 ], an oscillatory pattern at the cardiac frequency. CO is often observed visually during an absence of diaphragmatic activity (e.g. expiration) and the detection of CO has been suggested as a tool for identifying central sleep apnea. The source of CO is debated, with some arguing that the phenomenon arises from the physical transmission of heartbeats to the trachea, mediastinum, or lungs via direct contact [ 1 - 3 ], while other experimental evidence suggests that it originates from cyclic changes in pulmonary artery pressure and flow [ 4 , 5 ]. Despite these debates, analysis of CO has found applications in various clinical contexts. Its physiological implications, such as its role in gas mixing within the lungs [ 6 ] and lung mechanics assessment [ 7 ], have been extensively studied. The detection of CO during apnea has been suggested as an index distinguishing central from obstructive apnea events [ 8 , 9 ]. This property has been applied to scoring apneas as well as having potential in the design of automated CPAP titration systems for sleep apnea patients. Additionally, detection of CO has been leveraged to develop algorithms for flow-triggered mechanical ventilation [ 10 ], as well as having applications in automatic sleep staging [ 11 ], determining anesthetic depth [ 12 ], and estimating heart rate surrogates [ 13 ]. Given the potential applications of CO, there is a strong incentive to extract it from airflow signals for further investigation, particularly in settings where airflow is the only available channel. However, to our knowledge, most existing CO-related algorithms focus either on suppressing the impact of CO as unwanted “noise” or on simply estimating heart rate using the airflow signal. Examples include bandpass filtering (BPF) [ 14 ], time-frequency analysis, [ 12 , 13 ] and adaptive filtering [ 15 , 16 ] (when a cardiac-related channel is available). In contrast, there has been limited effort toward recovering CO when airflow is the only available signal. From a signal processing perspective, two key characteristics of the airflow signal pose significant challenges to CO recovery. First, respiratory signals and CO exhibit time-varying frequency, amplitude, and non-sinusoidal oscillatory patterns due to breathing pattern variability [ 17 ] and heart rate variability (HRV) [ 18 ] respectively. Second, the non-sinusoidal oscillatory patterns lead to spectral interference of the respiratory signal and CO. Specifically, the spectral content of the respiratory signal ranges from about 0.2Hz to 4Hz, and that of CO ranges from about 1Hz to 10Hz. This spectral interference is worsened when the ratio of heart rate and respiratory rate is low or when the CO amplitude is small relative to breath size. These inherent time-varying natures limit the effectiveness of the traditional BPF approach. Modern time-frequency (TF) analysis methods, such as synchrosqueezing transform (SST) [ 12 ] and the de-shape algorithm explored in previous CO studies [ 13 ], are designed to analyze and decompose signals with these time-varying frequency properties. However, the presence of apnea events (i.e., total absence of respiration for ≥10 seconds) introduces additional complexity. These events disrupt regular breathing cycles and are associated with irregular “gasps,” large breaths extraneous to the native respiratory rate that introduce further spectral leakage, significantly reducing the effectiveness of TF analysis tools. Consequently, there is a clear need for an algorithm specifically designed to handle such complex signals. The primary contribution of this paper is the development of a novel airflow decomposition algorithm, termed the Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA), which is designed to decompose CO from the given airflow signal. The main idea underlying OSADA is the addition of a novel low rank matrix denoising algorithm, eOptShrink [ 19 ], to a TF analysis based decomposition algorithm, shape-adaptive mode decomposition (SAMD) [ 20 ], to achieve an accurate CO decomposition. Intuitively, eOptShrink is a data-driven lowpass filter, where the filter’s basis is not the conventional Fourier basis but is instead derived from the airflow signal itself. This airflow-driven basis offers a significant advantage: it reduces the impacts of spectral interference of the respiratory signal and the relatively weak CO, takes care of the time-varying amplitude and frequency, and mitigates the impact of apnea events on the TF analysis, enabling effective recovery of both the respiratory signal and CO. The established theory underlying TF analysis tools, along with recently established theorems [ 19 ] within the framework of random matrix theory [ 21 , 22 ] supporting eOptShrink, provide OSADA with solid and interpretable theoretical foundations. Here, we present a comprehensive series of numerical evaluations conducted on both real and simulated datasets to demonstrate OSADA’s potential for real-world airflow signal analyses. 2. Mathematical model Respiration and its representation by airflow signals is quasi-periodic: Amplitude (the size of individual breaths), frequency (how long individual breaths last), and oscillatory pattern (breath shape) change over time [ 17 ]. The time-varying frequency, or instantaneous frequency, amplitude modulation and the time-varying oscillatory pattern comes from the breathing pattern variability [ 17 ]. The time-varying oscillatory pattern reflects at least the dynamics of the laryngopharyngeal structure [ 23 - 25 ]. For example, compared with normal breathing, during a reduction in flow that lasts several breaths due to upper airway obstruction (i.e., an obstructive hypopnea) inspiratory time gets longer with a flattened peak [ 24 ]. Further, during the inspiratory phase, lung expansion elongates the trachea resulting in stiffening of the trachea [ 26 ], which has been suggested as a potential source of dampening CO in the airflow signal. To encode the above observations, consider the following phenomenological model to describe the airflow signal. Let t l represent the starting time of expiration of the l th breath cycle (defined here as the start of expiration to the end of inspiration), so that … < t l < t l+ 1 < t l+ 2 < …. We choose to model the airflow signal using the expiration initiation time to simplify the upcoming algorithm design, as it is numerically more stable to detect. Note that reflects how fast the l th breath cycle is, or the notion of instantaneous frequency. It is important to note that during a pause in breathing, such as occurs in sleep apnea, the sequence t l +1 − t l can be very nonuniform. The airflow signal is represented as a random process by the following equation: where ∑ l ∈ℤ A l s l ( t − t l ) models the diaphragm-derived airflow (DDairflow) signal—that is, the component of the airflow signal driven by diaphragmatic activation (inspiration) and relaxation (expiration). A c ( t ) s c / ϕ c ( t )1 models the CO and ξ ( t ) models the noise. See Figure 1 for an illustration of this model. Physiologically, the summation of DDairflow and CO is airflow without noise contamination or clean airflow . When respiration alone is the specific phenomenon of interest, DDairflow is the desired signal, and CO is viewed as an artifact to suppress. Below, we entail conditions for this model. Download figure Open in new tab Figure 1. An illustration of the airflow signal and the proposed phenomenological model for a subject experiencing central (i.e. non-obstructive) sleep apnea. The top two gray curves represent the simultaneously recorded abdominal (ABD) and thoracic (THO) movement, indicating that no respiratory effort is detectable during the apnea, making it central. The black curve represents a realization of the random process Y ( t ) that models the airflow signal, with the detected starting time of expirations superimposed as red circles. The red crosses superimposed on the apnea represent where breaths would be expected based on prior respiratory frequency. These are used in Step 1 of the proposed OSADA algorithm. The airflow signal is modeled as the sum of three components: the DDairflow signal (blue), cardiogenic oscillations (CO, cyan), and noise plus artifacts (red). At the bottom, the simultaneously recorded photoplethysmogram (PPG) is shown in gray. From the PPG vertical lines are drawn at maximal slopes during systolic phases to visualize cardiac timing. Within the DDairflow component, A l > 0 describes the magnitude of the l th breath cycle, and s l is a compactly supported smooth function that describes the oscillatory pattern of the associated breathing. We assume s l has unit L 2 norm, that (according to convention) inhalation is upward and expiration is downward, and that s l (0) is the airflow at the starting time of the expiration. We call s l the wave-shape function (WSF) [ 27 ] associated with the l th breathing cycle. s l itself is not random but follows basic rules guiding the breathing pattern, which means it has a locally low-dimensional structure (mathematically, we can model it by the wave-shape manifold [ 28 ]). This is the key assumption guiding our algorithm design. For the CO component, A c ( t ) ≥ 0 is a smooth function describing the size of CO, s c ) is a 1-periodic function that describes the oscillatory pattern of CO. Analogous to the DDairflow component, we designate s c as the WSF of the CO. ϕ c ( t ) is a smooth and monotonically increasing function describing the phase of CO, so that describes the instantaneous heart rate (IHR). Finally, ξ ( t ) is the inevitable noise, which we assume to be locally stationary [ 29 ] with mean 0 and finite fourth moments after discretization at each time. Note that in this model, we use one WSF to model the CO and the noise is assumed to be locally stationary to simplify the discussion. More general models can be considered if needed. Remark We shall comment on a technical detail. When a subject breathes normally (i.e. during stable sleep in the absence of apneas, hypopneas, or arousals from sleep), the model for DDairflow, ∑ l ∈ℤ A l s l ( t − t l ), can be replaced by a simpler model A r ( t ) s r / ϕ r ( t )1 (e.g., [ 30 ]), where A r ( t ), s r and ϕ r ( t ) have the same meaning as those for CO. If the goal is to estimate the instantaneous breathing rate and instantaneous heart rate (IHR), this model is sufficient, and several TF analysis tools are available (e.g. [ 13 ]). However, if the goal is to accurately recover CO, this simpler model and the associated TF analysis algorithms are insufficient due to the spectral interference. 3. Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) The proposed airflow denoise algorithm, coined optimal shrinkage-aided airflow decomposition algorithm (OSADA), is illustrated in Figure 2 . The key technical ingredients are optimal shrinkage and SAMD [ 20 ], a data-driven method to enhance the decomposition accuracy for accurate CO recovery. This data-driven method is constructed using established physiological knowledge and its role is similar to that of Fourier basis commonly used in signal processing. Download figure Open in new tab Figure 2. An illustration of the proposed OSADA algorithm. The algorithm is composed of three main blocks. First, after standard procedures for pre-processing the airflow signal, SST and SAMD are applied to suppress the impact of CO on eOptShrink. Second, the cleaned airflow signal is segmented into breath cycles using a breath detection algorithm [ 31 , 32 ]. Using K-means clustering, breath cycles are then categorized into K groups, and the extended optimal shrinkage (eOptShrink) algorithm [ 19 ] is applied to each group. eOptShrink, a practical implementation of the OptShrink algorithm [ 21 ] with theoretical support, effectively handles complex noise and reconstructs breath cycles. DDairflow is then reconstructed by stitching the resulting “cleaned” breath cycles together. The raw CO is obtained as the difference between the original airflow signal and reconstructed DDairflow. The final step is applying SST and SAMD to the raw CO to obtain the final CO signal. A detailed step-by-step breakdown of the algorithm follows below. Note that the following steps assume the airflow signal has already been reviewed by an expert with scoring for respiratory events (apneas and hypopneas). The Matlab implementation is available upon request. Step 0: Preprocessing Let the airflow signal sampled at F s = 400 Hz be denoted as . If the original sampling rate differs from 400 Hz, resample accordingly. For consistency, adjust the polarization of the airflow signal as necessary so that inhalation is represented as an upward deflection and exhalation as a downward deflection. Visually review the signal to identify dropout segments, characterized by a flat, featureless waveform for > 3 mins indicating an absence of physiological activity due to temporary subject disconnection. Further analysis focuses on the remaining segments. Following the recommendation in [ 31 ], filter the signal using a 4th-order Butterworth high-pass filter (0.1 Hz cutoff) to remove low-frequency drifts and censor artifacts, followed by a 4th-order Butterworth low-pass filter (10 Hz cutoff) to suppress muscle noise and high-frequency interference. The resulting signal is denoted f ∈ ℝ N . To reduce the impact of CO during further algorithm steps, estimate CO using SAMD, which requires an estimate of the amplitude modulation and phase of the fundamental component of the putative CO. To obtain these estimates, apply SST to f to obtain a time-frequency representation. Use the ridge extraction algorithm to detect the dominant ridge in the spectral range from 1.6Hz to 4Hz in the time-frequency representation. Then apply the reconstruction formula [ 33 ] to obtain the desired estimate. Finally, apply SAMD to f with the estimated amplitude modulation and phase of the fundamental component of the CO to obtain f c , SAMD ∈ ℝ N , which is an estimate of CO. The resulting CO-reduced airflow signal is denoted f L ≔ f − f c , SAMD ∈ ℝ N . See Discussion for more details of this design. Step 1: Breath cycle detection and segmentation Apply a breath detection algorithm [ 31 , 32 ] to f L to identify the start of expiration, denoted t l for the l th breath cycle. Suppose we obtain m intact breathing cycles. For the l th breath cycle, EE l designates the interval between t l and t l− 1 (the start of expiration for the previous ( l − 1)th breath cycle). Define p 1 as the rounded integer of 3/4 of the 90% percentile of , where | EE l | is the number of sample points in the interval EE l . Segment the airflow signal into intervals of length , with each interval centered on the start of expiration. Save the l -th breath cycle as a p -dim vector For some pairs of consecutive segments, t l − p 1 > t l− 1 + p 1 , meaning the two consecutive segments will neither abut nor overlap, but have a “gap” between them. This is particularly likely in the presence of apnea events. Suppose there are K such segments. If K > 0 and the k -th segment happens at time , we divide the segment into overlapping segments of equal length with the cut point , where l = 1, …, m k , and construct m k vectors See Figure 1 for an illustration of these cut points. We construct extra p -dim vectors, each of which may be thought of as an “inferred” breath cycle to fill the “gaps” between non-adjoining consecutive native breath cycles. The total number of breath cycles, both native and inferred, is designated n , such that . The n breath cycles are stored in a data matrix called breath matrix where p is the length of a segment. Note that by construction, the ( p 1 + 1)th row of the first m columns of X is associated with the initiation of expiration. The number p depends on the sampling rate and the respiratory rate. Step 2: Denoise the breath cycle matrix by optimal shrinkage X can be viewed as a noisy data matrix due to contamination by the complicated noise matrix Φ and the residue of CO. The recently developed matrix denoising algorithm eOptShrink [ 19 ] is specifically intended to denoise such a data matrix. Without loss of generality, assume p ≥ n and denote β n ≔ p / n . If p < n , we simply transpose X and apply the same procedure. Denote the singular value decomposition of X as , where the singular values are indexed in decreasing order, and and are the associated left and right singular vectors. Denote . Recall that is the eigenvalue of the positive definite matrix XX T ordered in the decreasing order. eOptShrink is composed of four main steps. First, set and estimate the noise strength, or the rightmost spectral bulk edge under a mild assumption (further discussed below) by where [ n c ] is the closest integer to n c . With the knowledge of the noise strength, we estimate the effective rank , which is the number of greater than . In other words, if is sufficiently large, we view it as a signal. It is called effective rank since, in general, might be much smaller than r given weak signal components could be buried by the noise. Second, take as a base and modify eigenvalues weakly than by applying where . The purpose of this correction is to suppress the impact of signal on these eigenvalues since they are close to the rightmost spectral bulk edge. This modification leads to a more accurate estimate of Φ’s spectral structure. Third, for the top eigenvalues, numerically implement the D-transform associated with Φ (see (15) and (16) in [ 19 ]), and set . Next, numerically implement Stieltjes transforms associated with Φ by Next, compute numerical derivatives of the D transform and Stieltjes transforms associated with Φ by constructing and where Mathematically, ( respectively) approximates the inner products between the true singular vector u i ( v i respectively) and the noise-contaminated singular vector ( respectively). Clearly, when noise gets weaker, we expect this inner product to get close to 1 (see the statement below Proposition 2.6, (28) and Theorem 3.5 in [ 19 ]). The above leads to an estimate of the spectral structure of the high dimensional noise Φ, which depends on the relationship between the clean singular values and singular vectors and the noisy singular values and singular vectors. Last, denoise the matrix X via the formula Here, is called the optimal shrinkage of the singular value with the Frobenius norm as the cost function, which corrects by taking the biased singular vector estimates into account. Note that is the denoised i th breathing cycle. Step 3: Recover the DDairflow signal We recover the DDairflow signal by stitching all denoised breath cycles in with the overlap-and-add method. Denote the recovered DDairflow signal as . The difference between the pre-processed signal and DDairflow, ,is called the raw CO, which includes the putative CO, pure noise, and possible remaining artifacts. Step 4: Recover the cardiogenic oscillation Finally, apply SAMD [ 20 ] to recover CO from the raw CO, . We apply SST to following the same procedure in Step 0 to obtain the phase and amplitude information of the fundamental component of CO. Then, apply SAMD to following the same procedure in Step 0 to obtain the final estimate of the putative CO, denoted as . Before moving on, we would like to make a few points on the philosophy underlying this algorithmic design: In Step 0, we consider the spectral range from 1.6 Hz to 4 Hz when applying the ridge extraction algorithm. While this range is clearly too high for the heart rate, it is intentionally chosen to capture the second harmonic of CO, thereby minimizing spectral interference when the ratio between heart rate and respiratory rate is low. Notably, if PPG is available to provide heart rate information, this ridge extraction can be skipped. Next, consider the essential features of the matrix X produced in Step 1. Based on the proposed phenomenological model, X = D + Φ, where the data matrix D = [ d 1 d 2 … d m 0 … 0] ∈ ℝ p×n consists of DDairflow cycles, and the noise matrix Φ = [ ς 1 ς 2 … ς n= ] ∈ ℝ p×n consists of noise. To be more precise, the l -th detected breathing cycle modeled by A l s l ( t − t l ) is saved as d l , and the bandpassed noise ξ ( t ) over the same interval of d l , from t i − p 1 to t i + p 1 , is saved in ς l . A critical observation is that due to the similarity of breathing cycles, the data matrix D is a low rank matrix; that is, D is a low rank matrix and its singular value decomposition is , where r is much smaller than p and n , and σ i is indexed in decreasing order. Recall that the left singular vector u i is a data-driven orthogonal basis for the breathing cycle, and the right singular vector v i records the coordinates of all breathing cycles; that is, the j -th entry of v i is the coordinate of d j associated with u i . On the other hand, σ i encodes the global strength of the data matrix associated with the i -th data-driven basis. We could therefore view as the i -th signal component and σ i as its strength. The ultimate goal is to recover D from X . Due to the nonstationary nature of ξ ( t ), the noise matrix Φ does not have independent entries. To capture the dependence structure, we assume that Φ has a separable covariance structure; that is, Φ = A 1/2 ΞB 1/2 , where A and B are positively definite and Ξ contains independent entries with mild moment conditions. This structure is chosen to strike a balance between theoretically supported algorithmic design and practical constraints due to noise. 4. Material and statistics Semi-real simulated database To quantify the performance of OSADA and compare OSADA with other algorithms, we constructed a database of semi-real simulated airflow signals. The simulated airflow signals, composed of DDairflow, CO, and noise, were derived from 5 subjects with varying severities of sleep apnea who had simultaneously recorded airflow and PPG signals. The goal was to retain as many real airflow dynamics as possible, including breathing pattern variability and apnea events. The simulated DDairflow, CO and noise are publicly available at https://github.com/hautiengwu . Generate a simulated DDairflow signal Visually review the original 400 Hz (resample as appropriate) airflow signal, denoted , marking segments where airflow is completely absent, out of range, or interrupted by obvious artifacts. Use apnea events to define a baseline of 0 for the airflow signal. For non-apnea segments, including hypopnea, use a breath cycle detection algorithm to determine the start of expiration for each breath cycle, followed by manual verification. (Manual annotation was performed by H.-T. Wu with assistance from D. M. Rapoport and T. M. Tolbert.) The resulting annotated signal is denoted . The marked artifacts are then removed, and the remaining, non-artifact segments are concatenated as follows: Suppose K ≥ 1 artifact segments are removed. For the k th segment, denote s k as the start of expiration for the last intact breath cycle before the start of the k th segment. Denote t k as the start of expiration of the first intact breath cycle after the end of the k th segment. For all K segments, the interval between s k and t k is removed from f 1 . The resulting signal is denoted f 2 ∈ ℝ N , where N ≤N 0 . Next, construct breath cycle segments from f 2 following the same Step 2 of OSADA, but use the manually marked start of expiration timings and let p 1 equal the rounded integer of 3/4 of the 90 th percentile of the total length of consecutive pairs of breath cycles over segments without apnea events. For the i th segment, x i , determine its 50 nearest neighbors in the H 1 norm (a combination of L 2 -norms of the vector together with its derivative), denoted as . Construct y i to be the first column of the top singular pair of the matrix and view y i as the cleaned-up x i . Stitch y i to form the final DDairflow signal f ∈ ℝ N . Next, run a validation to confirm there is no remaining CO: Use the simultaneously recorded photoplethysmogram (PPG) from the PSG to estimate the inverse cardiac phase, , denoted as , by linearly interpolating the data set {( t k , k ), k = 1, … n }, where t k represents the maximal slope during the systolic phase of the k th cardiac cycle from the PPG. Warp f by to obtain f w ∈ ℝ M , sampled also at 400 Hz. Plot the power spectrum of f w . When CO remains, a strong peak is observed at k Hz, where k ∈ ℕ. When there is no obvious peak at k Hz, where k = 1, …,6, take the associated f to be the final simulated DDairflow. Generate a simulated CO signal The CO is simulated by using the simultaneously recorded PPG signal and the CO profile is determined from central apnea events in the following steps. First, apply the signal quality index designed for PPG to determine and concatenate high-quality segments. The maximal slope time point during the systolic phase t k of the k -th cardiac cycle is used as the landmark. Suppose there are N cardiac cycles. Then, estimate the phase of the PPG, denoted as , by linearly interpolating the data set {( t k , k ), k = 1, … N }. Also, estimate the amplitude of PPG by applying the Hilbert transform followed by taking the absolute value and a 5-minute kernel smoothing to avoid potential artifacts. Denote the resulting amplitude as . Warp the airflow signal using the cardiac phase estimated from the simultaneously recorded PPG so that the CO has constant frequency 1Hz. Obtain a segment of CO from a central apnea event with visually clear CO, denoted as h ∈ ℝ m . Estimate the amplitude and phase of harmonics of by the trigonometric regression; that is, fitting 1, sin (2 πlj + γ l ) and cos (2 πlj + γ l ), where l = 1, …, L , and L ∈ ℕ is the number of harmonics, to h by running the linear regression, and denote the estimated amplitude as b l ≥ 0 and global phase as γ l ∈ [0,2 π ). Note that b l and γ l encode the CO harmonic profile; that is, the WSF of CO. Finally, simulate CO, denoted as g ∈ ℝ n , by setting We use 5 harmonics since most CO we have observed, particularly those during the central apnea events, have 5 or fewer harmonics. Generate a noise signal Create simulated noise by applying a BPF with the spectral band 0.1Hz to 10Hz to a white Gaussian noise, denoted as ξ ∈ ℝ n . Add DDairflow,CO, and noise Take ϑ 1 , ϑ 2 ∈ ℝ as the desired signal to noise ratio (SNR) between CO and DDairflow, as well as the desired SNR between CO and noise. Create a semi-real “complete” airflow signal by summing one simulated DDairflow to one simulated CO scaled by C 1 > 0 and one simulated noise signal scaled by C 2 > 0. C 1 and C 2 are chosen so that , where SD indicates the standard deviation, and .In the end, we obtain 25 semi-real simulated signals. Real database We assessed the performance of OSADA on real-world data from 22 subjects from a single-center prospective observational study at the sleep center in Taipei Veterans General Hospital, Taipei (VGHTPE), Taiwan, collected from June to December 2023 with Medical Ethics Committee approval (IRB No: 2023-04-003A). The study adhered to the ethical guidelines of the 1975 Declaration of Helsinki. PSG data were collected using Alice 6 with Sleepware G3 Version 3.9.5 (Philips, Netherlands). Airflow signals were collected using a pressure transducer/nasal cannula with sampling rate 100Hz. Less than 0.5% of all airflow recordings consisted of out-of-range artifacts. We visually inspected the airflow signal of each subject to find segments with visible CO, ultimately obtaining 86 segments with length 20.48 ± 23.33 (range: 2-129) minutes from 3 subjects with severe sleep apnea (apnea-hypopnea index, AHI≥30 events/hr), 1 subject of moderate sleep apnea (30>AHI≥15 events/hr), 7 subjects of mild sleep apnea (15>AHI≥5 events/hr) and 11 normal subjects (AHI<5 events/hr). Comparison While algorithms exist for extracting heart rate information fromairflow, to our knowledge, no algorithms have been specifically designed to decompose CO. We consider two algorithms that aim to reduce CO’s impact on the airflow signal as CO decomposition methods and compare their performance with the proposed OSADA algorithm. BPF: A 4th-order Butterworth filter is applied with a subject-specific spectral band, determined based on the simultaneously recorded PPG signal. The instantaneous heart rate (IHR) is extracted from the PPG, and its range defines the spectral band for each subject. To recover DDairflow, the estimated CO is first subtracted from the airflow signal, followed by applying a high-pass filter (0.1 Hz cutoff) and a low-pass filter (10 Hz cutoff) as suggested in [ 31 ]. Savitzky-Golay smoothing filters (SGF) [ 34 ]: This algorithm, also known as the polynomial smoothing filter, works by minimizing the least-squares error when locally fitting a polynomial to noisy data. SGF can be applied to recover DDairflow [ 35 ] as it is well-known for preserving the high-frequency content of the clean signal while exhibiting some low-pass filtering properties. Here, we use a polynomial order of 2 with a fixed fitting frame size of 20 seconds [ 35 ]. The difference of the airflow signal and the recovered DDairflow can be viewed as a noisy estimate of CO. Finally, this noisy CO estimate is refined by applying a 4th-order Butterworth filter with a subject-specific spectral band, determined based on the simultaneously recorded PPG signal, to produce a final estimated CO. Statistics For the semi-real simulated database, we consider the following metrics to compare the performance of difference algorithms. The normalized root mean square error (NRMSE) of the recovered CO is evaluated over non-overlapping 30-minute windows over the length of the study. For each 30-minute window, the NRMSE is defined as , where g ∈ ℝ n is the true CO and is the estimated CO using different algorithms. Similarly, the NRMSE for the DDairflow recovery is calculated. These metrics quantify the performance of the algorithms in the L 2 -norm sense. An ideal CO or DDairflow recovery should have a small NRMSE. Next, we quantify the adequacy of CO removal from the airflow signal. We apply the same warping procedure described above to each recovered DDairflow, denoted as ,, and denote the warped signal as , also sampled at 400 Hz. Under the phenomenological model, if residual CO remains in the airflow signal, peaks are observed in the power spectrum of concentrated at 1, …,4 Hz. The energy of remaining CO in over each 10-minute segment is defined as the spectral energy of at 1, …,4 Hz, denoted as E a . Then, apply the same warping procedure to the simulated clean DDairflow over each 10-minute segment, and calculate its spectral energy at 1, …,4 Hz, denoted as T a . The final index is denoted ϱ N ≔ E a / T a . An ideal removal of CO should have ϱ a = 1. For the real database, since we do not have a reference standard measure for CO, we propose to quantify how accurately the CO is recovered from the airflow signal by a warping procedure similar to that used in evaluating ϱ a . Warp by to obtain , also sampled at 400 Hz. If CO is well recovered, spikes are observed at 1, …,4 Hz in the power spectrum of . Denote the energy of recovered CO over each segment by evaluating the spectral energy of at 1, …,4 Hz, denoted as E c . To generate a baseline, random permute and warp the resulting signal by to obtain , also sampled at 400 Hz. Denote the spectral energy of at 1, …,4 Hz as E r . The final index is denoted ω c ≔ E c / E r . An ideal removal of CO should have ω c ≫ 1. The performance of various algorithms was evaluated using the Wilcoxon signed-rank test, with statistical significance defined as p < 0.05. To account for multiple comparisons, the Bonferroni correction was applied. 5. Results An example from a subject without sleep apnea is illustrated in Figures 3 and 4 . Figure 3 is a 60-second segment including the original airflow signal, the denoised airflow signal with CO removed, and the reconstructed CO signal. CO is not visibly apparent in the original airflow signal, aside from minor fluctuations resembling noise. After applying OSADA, the raw CO, which contains noise and other artifacts, displays oscillatory behavior that aligns with the oscillatory periods of the simultaneously recorded PPG. Visually, approximately two cycles appear in the raw CO over one cardiac cycle as determined by the PPG. As shown in Figure 4 , this pattern is reflected in the time-frequency representation (TFR) of the raw CO determined by SST: The dominant curve of the TFR coincides with the double of the PPG-derived IHR. Conversely, there is no dominant curve (or at most there is a weak curve) in the TFR coinciding with the IHR. Additionally, prominent curves align with the triple and quintature IHR. These observations suggest that the raw CO has a fundamental frequency corresponding to the IHR measured by the PPG, though its fundamental component is either missing or too weak to be visually detected. See Discussion section for a more detailed analysis of this finding. Download figure Open in new tab Figure 3: A 60-second segment of the simultaneously recorded abdominal movement, thoracic movement and airflow signal recorded from a subject without apnea events are shown as the two gray curves, black curve and blue curve. In the bottom, the denoised airflow signal by OSADA is shown in blue, the raw CO estimate by OSADA is shown in red, the denoised CO with SAMD is shown in magenta, and the simultaneously recorded PPG is shown in blue. The raw CO and final CO are shown with the magnitude amplified by three and nine respectively to enhance the visualization. Gray vertical lines mark the timings of maximal slopes during the systolic phases, overlaid in gray to enhance the visualization of cardiac oscillations. Download figure Open in new tab Figure 4: The TFR determined by SST of the raw CO shown in Figure 3 is shown in both panels. The signal shown in Figure 3 is the first 60 seconds of the signal used to generate the TFR in this plot. The IHR determined by the simultaneously recorded PPG and its double, triple, quadrature and quintature are superimposed in the bottom panel as red, blue, green, magenta and cyan curves. Figures 5 and 6 analogously illustrate an example from a subject with central apneas. In the 60-second segment shown in Figure 5 , CO is visibly present in the airflow signal during the apnea event as minor fluctuations resembling noise. The raw CO obtained by OSADA, which contains noise and other artifacts, displays oscillatory behavior that aligns with the oscillatory periods of the simultaneously recorded PPG, especially in the second half of the signal. Unlike the example in Figure 4 , the TFR of CO shown in Figure 6 exhibits dominant curves that coincide with the IHR and its multiples, particularly the triple. Download figure Open in new tab Figure 5: A 60-second segment of the simultaneously recorded abdominal movement, thoracic movement and airflow signal recorded from a subject with central sleep apnea events are shown as the two gray curves, black curve and blue curve. In the bottom, the denoised airflow signal by OSADA is shown in blue, the raw CO estimate by OSADA is shown in red, the denoised CO with SAMD is shown in magenta, and the simultaneously recorded PPG is shown in blue. The raw CO and denoised CO are shown with the magnitude amplified by three to enhance the visualization. Gray vertical lines mark the timings of maximal slopes during the systolic phases, overlaid in gray to enhance the visualization of cardiac oscillations. Download figure Open in new tab Figure 6: The TFR determined by SST of the raw CO shown in Figure 5 is shown in both panels. The signal shown in Figure 5 is the first 60 seconds of the signal used to generate the TFR in this plot. The IHR determined by the simultaneously recorded PPG and its double, triple and quadrature are superimposed in the bottom panel as red, blue, green and magenta curves. Next, we analyze an example involving obstructive apneas. Since there is no segment containing both obstructive apneas and potential CO in the airflow in the real database, we instead select a representative airflow signal during an obstructive apnea to illustrate the results of OSADA. Figure 7 presents a 60-second segment of the airflow signal during an obstructive apnea, along with the corresponding denoised airflow and raw CO signal. Notably, no visible signs of CO are present in the airflow signal during the apnea. Additionally, the raw CO extracted by OSADA consists only of noise. Note that the quality of the denoised airflow in this case might be considered lower than in Figures 3 and 5 , as some deformations are visually apparent (indicated by magenta arrows). These deformations cause significant fluctuations in the raw CO during breathing cycles before and after obstructive apneas. In the absence of a reference standard measurement for CO, we cannot definitively determine the source of these deformations. However, this observation, along with previous findings, is consistent with an understanding that CO behaves differently in the setting of central and obstructive apneas. Download figure Open in new tab Figure 7: A 60-second segment of the simultaneously recorded abdominal movement, thoracic movement and airflow signal recorded from a subject with obstructive sleep apnea events are shown as the two gray curves, black curve and blue curve. In this example, there is no visually observable CO in the airflow signal. In the bottom, the denoised airflow signal by OSADA is shown in blue, the raw CO estimate by OSADA is shown in red, and the simultaneously recorded PPG is shown in blue. Since there is nothing to denoise, we do not apply SAMD. The raw CO is shown with the magnitude amplified by three to enhance the visualization. Gray vertical lines mark the timings of maximal slopes during the systolic phases, overlaid in gray to enhance the visualization of cardiac oscillations. To quantify the performance of OSADA and compare it with other algorithms, we consider NRMSE, ϱ a , and ω c for the semi-real simulated database and ω c for the real database. Figure 9 presents violin plots of all considered indices applied to the simulated database with different ϑ 1 and ϑ 2 . The indices are listed in Table 1 . In terms of NRMSE, SGF performs consistently better than BPF in the recovery of DDairflow, but worse in the recovery of CO. SGF also better recover CO in terms of ρ c . OSADA outperforms other algorithms in all indices. Viewing p<0.05 as statistical significance, the differences in performance were statistically significant under the Bonferroni correction. The performance of SGF is consistent with a benefit from polynomial fitting in this application, particularly in recovering the high frequency features of DDairflow. Further, in the presence of greater noise or a weaker CO signal, the recovery performance of OSADA is worse, consistent with our above intuition. View this table: View inline View popup Download powerpoint Table 1. A summary of NRMSE of DDairflow, NRMSE of CO, ϱ a and ω c for the semi-real simulated database with various ϑ 1 and ϑ 2 . Download figure Open in new tab Figure 8: The TFR determined by SST of the raw CO shown in Figure 7 is shown in both panels. The signal shown in Figure 7 is the first 60 seconds of the signal used to generate the TFR in this plot. The IHR determined by the simultaneously recorded PPG and its double, triple and quadrature are superimposed in the bottom panel as red, blue, green and magenta curves. Download figure Open in new tab Figure 9: Left to right columns: the violin plot of NRMSE of DDairflow recovery, NRMSE of CO recovery, ρ c and ω c of different algorithms in the simulated database. From top to bottom rows: ϑ 1 = 10, ϑ 2 = 10, ϑ 1 = 10, ϑ 2 = 20, ϑ 1 = 15, ϑ 2 = 10 and ϑ 1 = 15, ϑ 2 = 20. Figure 10 shows the violin plot of ω c for the real database, where ω c = 3.44 ± 0.94, 3.43 ± 0.92, and 8.26 ± 3.01 for BPF, SGF and OSADA, respectively. In all cases, OSADA outperforms other algorithms with statistical significance under Bonferroni correction. Download figure Open in new tab Figure 10: The violine plot of ω c of different algorithms in the real database. 6. Discussion We propose and evaluate the OSADA algorithm, which is designed to decompose CO from airflow in the absence of other signals. Its efficiency is demonstrated in the Results section. OSADA has several potential practical applications. First, it generates a clean DDairflow signal with preservation of high frequency features; that is, peaks and troughs are preserved through matrix denoising. Retention of these features with suppression of non-respiratory related oscillations may enhance the accuracy of breathing dynamics assessment, such as the precise timing of inspiration and expiration and the expiratory-to-inspiratory ratio [ 36 ], as well as techniques for estimating pharyngeal critical closing pressure (Pcrit) or collapsibility based on breath size and shape [ 37 - 40 ]. Further, the combination of the clean DDairflow and CO might help improve automated CPAP titration algorithms (autoCPAP) that analyze the details of the shape of the airflow signal and can be confused by CO. Similarly, the characteristics of CO may be of significant use in assessing if pauses in breathing (apnea) are “central” or “obstructive” in nature, even in the absence of additional monitoring of respiratory effort (typically through motion sensors or plethysmography of the thorax and abdomen). This distinction has utility in the study of sleep-disordered breathing in both research and clinical settings. There is also the possibility that a clean CO signal may contain information useful in deducing sympathetic tone, useful in automated staging of sleep depth. BPF has previously[ 14 ] been used to remove CO from the airflow signal, enabling the extraction of HR information. However, spectral interference caused by breathing pattern variability or spectral leakage caused by the presence of apnea events degrades BPF performance. This limitation is well known in the literature (see, for example, [ 35 ]) and was observed for the semi-real simulated signal analyzed in the current study. Compared with BPF, SGF performs better in smoothing a noisy signal and preserving high frequency information when the clean signal has a broad spectrum. However, SGF also has limited performance in CO decomposition, leading to a suboptimal recovery of a clean DDairflow signal. A potential alternative approach might employ locally weighted polynomial regression [ 41 ] with data-driven bandwidth, but exploration of this possibility is outside the scope of this paper. SAMD is based on time-frequency analysis tools like SST, which can be viewed as a nonlinear variation of the trigonometric regression [ 42 ] that is suitable for extracting CO. Designed to capture time-varying frequencies associated with HRV, it overcomes a key challenge faced by BPF. However, though SAMD outperformed BPF in the current study, spectral interference between DDairflow and CO remains due to the nonsinusoidal oscillation of DDairflow, as well as spectral leakage of DDairflow caused by apnea events. OSADA overcomes these barriers using a data-driven approach that preserves the DDairflow signal as much as possible. This allows SAMD, in Step 4, to achieve a more accurate CO recovery. Among the existing methods for suppressing the impact of CO, the approach most similar to OSADA is [ 35 ], as both fall under the umbrella of template subtraction. In [ 35 ], the authors used R-peaks from the electrocardiogram (ECG) signal (simultaneously recorded with airflow) to segment impedance pneumography and construct CO templates through averaging. These CO templates can be viewed as the WSF of cardiac cycles in CO. The CO impact is then reduced by subtracting the signal using these templates. Compared with that averaging processing, eOptShrink in OSADA is a data-driven way to optimally find the WSF of respiratory cycles from airflow when an ECG or PPG signal is not available. Further, the template subtraction method in [ 35 ] implicitly assumes the WSF of CO is relatively simple and can be recovered through averaging. We hypothesize that by incorporating the low-rank structure of the WSF of CO, eOptShrink may offer a more effective way to handle complex WSF patterns of CO when an ECG or PPG signal is available. Exploring such an approach, along with more general methods like adaptive filtering and sensor fusion, will be the focus of our future work. In the recovered CO shown in Figure 3 , it may appear that the CO signal oscillates twice per cardiac cycle. However, this is misleading, as the two oscillations are not identical. For example, in the latter half of the depicted signal, the paired oscillations differ: the first has a sharper peak than the second. This morphological feature is controlled by higher-order harmonics. Mathematically, a signal is 1-periodic if and only if the greatest common divisor (gcd) of its nonzero harmonic orders is 1. In this case, the CO has visible second, third, and fifth harmonics (see Figure 4 ). Since the gcd of 2, 3, and 5 is 1, CO oscillates at a frequency similar to the IHR. By contrast, the CO example in Figure 5 includes a fundamental component in the TFR that aligns with the IHR, eliminating any visual ambiguity. Several technical details deserve discussion. p 1 in Step 1 of the algorithm was not optimized but was chosen in an ad hoc fashion based on the authors’ previous experience. Optimization may be warranted depending on the application. The main purpose of the upsampling step (i.e., for a signal with a sampling rate below 400Hz) is to enforce the alignment of breathing cycles and hence the low rank structure of the data matrix. To appreciate its significance, assume that the data matrix D = [ d 1 d 2 … d n ] = σ u 1 T ∈ ℝ p×n is rank one, where 1 ∈ ℝ n× 1 is a vector with all entries 1, and u ∈ ℝ p . Denote s k : ℝ p → ℝ p to be a cyclic shift operator by k ∈ ℤ steps. In general, if u is sinusoidal, that is, , a shift of u by k ≠ 0 steps can be linearly spanned by u and w ∈ ℝ p , where ; that is, . Thus, random cyclic shifts of all d j maximally increases the rank from 1 to 2 . However, the situation is significantly different when u is not sinusoidal. A direct expansion with trigonometric expansion shows that if u can be represented by L harmonics using Fourier series—that is, , where α j > 0 and β j ∈ [0,2 π ) —then random cyclic shifts of all d j may maximally increase the rank from 1 to 2 L . Thus, an accurate alignment is necessary to enforce the desired low rank structure for optimal shrinkage. The requirement for a low-rank structure partially accounts for the visual distortions observed in Figure 7 . It is well known that respiratory cycles surrounding obstructive apneas exhibit greater irregularity, though not necessarily chaotic behavior, compared to normal breathing. This irregularity contributes to the observed deformation. However, we emphasize that the exact source of this deformation remains uncertain. Factors such as increased respiratory drive or the body’s compensatory response may also play a role. Further research is needed to clarify these mechanisms. Further, these algorithms may struggle to effectively decompose CO when cardiac cycles are coupled with respiratory cycles in any manner, even in the absence of apneas. Mathematically, this presents an identifiability problem, and without additional information, we are not aware of any practical solution at this time. However, if a PPG or ECG channel is available, the coupling can be quantified, making it possible to modify OSADA to decompose the coupled CO and DDairflow. We leave this as a problem for future research. The construction of the semi-real simulation database warrants discussion. This database was generated using a cleaning procedure based on the proposed phenomenological model to preserve as many clinical respiratory features as possible, particularly in subjects with sleep apnea. The complexity and richness of airflow signals in apnea patients are well recognized in sleep medicine [ 23 - 25 , 37 ]. To our knowledge, while there have been efforts in this direction [ 43 , 44 ], a realistic mathematical model capable of generating such intricate airflow patterns—in addition to CO— has yet to be developed. Rather than focusing on mathematical model development, we chose to simulate the DDairflow and CO signals using a phenomenological approach for the purpose of algorithm evaluation, combined with a validation step and manual visualization. Denoising f 2 by using the top singular pair in the semi-real DDairflow signal generation is a critical step. We employ an approach similar to the nonlocal Euclidean mean (or median) [ 45 ] to cleanup noise and potential CO. However, rather than averaging all nearest neighbors (which often over-smooths the breathing cycle) or taking the median (which often leads to discontinuities) we extract the top left singular vector as the representative breathing cycle. This method strikes a balance between preserving the WSF of breathing cycles, including high-frequency components, and deleting CO and noise. While this technique is not suitable for CO recovery since it still poses a risk of over-smoothing the DDairflow signal, it effectively simulates a reasonably realistic DDairflow from patients with sleep apnea. In this study, segments with potential CO were manually selected through visual review by experts. An automatic system for CO recovery is needed for large-scale analyses, but several technical challenges remain to be addressed. First, an accurate signal quality assessment method must be established, particularly for evaluating the presence of CO and quantifying its quality. Determining whether CO exists based solely on airflow signals is inherently difficult, especially when its quality is obscured by noise. Statistically, developing an algorithm to detect the presence of CO aligns with the change-point detection problem for oscillatory signals—an area that remains underexplored with only a few recent studies addressing a simplified version of the challenge [ 46 ]. Assuming CO is correctly detected, an index is also needed to quantify its quality. Notably, signal quality assessment and change-point detection are closely linked and will likely need to be addressed together. These challenges will be a focus of our future research. This study has several limitations in addition to those mentioned above. First, we emphasize that the proposed phenomenological model-based simulation strategy used for performance evaluation cannot replace the mathematical modeling approach for other research purposes; it is simply an approach suitable for evaluating and comparing CO decomposition algorithms. Developing a realistic mathematical model capturing DDairflow and CO, which could be understood as constructing a digital twin of sleep apnea patients’ respiratory system, particularly the airflow, remains a challenging yet promising area for future research. Data Availability All data produced in the present study are available upon reasonable request to the authors Reference ↵ Arieli , R. : ‘ Cardiogenic oscillations in expired gas: origin and mechanism’ , Respiration Physiology , 1983 , 52 , ( 2 ), pp. 191 – 204 OpenUrl CrossRef PubMed Web of Science Bosman , A.R. , and Lee , G.D. : ‘ The effect of cardiac action upon lung gas volume’ , Clinical science , 1965 , 28 , pp. 311 – 324 OpenUrl PubMed Web of Science ↵ Heckman , J.L. , Stewart , G.H. , Tremblay , G. , and Lynch , P.R. : ‘ Relationship between stroke volume and pneumocardiogram’ , Journal of Applied Physiology , 1982 , 52 , ( 6 ), pp. 1672 – 1677 OpenUrl PubMed Web of Science ↵ Suarez-Sipmann , F. , Santos , A. , Peces-Barba , G. , Bohm , S.H. , Gracia , J.L. , Calderón , P. , and Tusman , G. : ‘ Pulmonary artery pulsatility is the main cause of cardiogenic oscillations’ , Journal of clinical monitoring and computing , 2013 , 27 , pp. 47 – 53 OpenUrl CrossRef PubMed ↵ Tusman , G. , Suarez-Sipmann , F. , Peces-Barba , G. , Climente , C. , Areta , M. , Arenas , P.G. , and Bohm , S.H. : ‘ Pulmonary blood flow generates cardiogenic oscillations’ , Respiratory physiology & neurobiology , 2009 , 167 , ( 3 ), pp. 247 – 254 OpenUrl PubMed ↵ Arieli , R. , Olszowka , A.J. , and Van Liew , H.D. : ‘ Postinspiratory mixing in the lung and cardiogenic oscillations’ , Journal of Applied Physiology , 1981 , 51 , ( 4 ), pp. 922 – 928 OpenUrl PubMed Web of Science ↵ Bijaoui , E. , Baconnier , P.F. , and Bates , J.H.T. : ‘ Mechanical output impedance of the lung determined from cardiogenic oscillations’ , Journal of Applied Physiology , 2001 , 91 , ( 2 ), pp. 859 – 865 OpenUrl PubMed Web of Science ↵ Ayappa , I. , Norman , R.G. , and Rapoport , D.M. : ‘ Cardiogenic oscillations on the airflow signal during continuous positive airway pressure as a marker of central apnea’ , Chest , 1999 , 116 , ( 3 ), pp. 660 – 666 OpenUrl CrossRef PubMed Web of Science ↵ Javaheri , S. , Rapoport , D.M. , and Schwartz , A.R. : ‘ Distinguishing central from obstructive hypopneas on a clinical polysomnogram’ , Journal of Clinical Sleep Medicine , 2023 , 19 , ( 4 ), pp. 823 – 834 OpenUrl PubMed ↵ Imanaka , H. , Nishimura , M. , Takeuchi , M. , Kimball , W.R. , Yahagi , N. , and Kumon , K. : ‘ Autotriggering caused by cardiogenic oscillation during flow-triggered mechanical ventilation’ , Critical care medicine , 2000 , 28 , ( 2 ), pp. 402 – 407 OpenUrl CrossRef PubMed Web of Science ↵ Lin , Y.-Y. , Wu , H.-T. , Hsu , C.-A. , Huang , P.-C. , Huang , Y.-H. , and Lo , Y.-L. : ‘ Sleep apnea detection based on thoracic and abdominal movement signals of wearable piezoelectric bands’ , IEEE journal of biomedical and health informatics , 2016 , 21 , ( 6 ), pp. 1533 – 1545 OpenUrl ↵ Lo , Y.-L. , Wu , H.-T. , Lin , Y.-T. , Kuo , H.-P. , and Lin , T.-Y. : ‘ Hypoventilation paherns during bronchoscopic sedation and their clinical relevance based on capnographic and respiratory impedance analysis’ , Journal of Clinical Monitoring and Computing , 2020 , 34 , pp. 171 – 179 OpenUrl PubMed ↵ Lu , Y. , Wu , H.-t. , and Malik , J. : ‘ Recycling cardiogenic artifacts in impedance pneumography’ , Biomedical Signal Processing and Control , 2019 , 51 , pp. 162 – 170 OpenUrl ↵ Wilson , A.J. , Franks , C.I. , and Freeston , I.L. : ‘ Methods of filtering the heart-beat artefact from the breathing waveform of infants obtained by impedance pneumography’ , Medical and Biological Engineering and Computing , 1982 , 20 , pp. 293 – 298 OpenUrl ↵ Schuessler , T.F. , Gottfried , S.B. , Goldberg , P. , Kearney , R.E. , and Bates , J.H.T. : ‘ An adaptive filter to reduce cardiogenic oscillations on esophageal pressure signals’ , Annals of biomedical engineering , 1998 , 26 , pp. 260 – 267 OpenUrl PubMed ↵ Yasuda , Y. , Umezu , A. , Horihata , S. , Yamamoto , K. , Miki , R. , and Koike , S. : ‘ Modified thoracic impedance plethysmography to monitor sleep apnea syndromes’ , Sleep Medicine , 2005 , 6 , ( 3 ), pp. 215 – 224 OpenUrl PubMed ↵ Benchetrit , G. : ‘ Breathing pahern in humans: diversity and individuality’ , Respiration physiology , 2000 , 122 , ( 2-3 ), pp. 123 – 129 OpenUrl CrossRef PubMed ↵ Electrophysiology , T.F.o.t.E.S.o.C.t.N.A.S.o.P. : ‘ Heart rate variability: standards of measurement, physiological interpretation, and clinical use’ , Circulation , 1996 , 93 , ( 5 ), pp. 1043 – 1065 OpenUrl FREE Full Text ↵ Su , P.-C. , and Wu , H.-T. : ‘ Data-driven optimal shrinkage of singular values under high-dimensional noise with separable covariance structure’ , Applied and computational harmonic analysis , 2024 ↵ Colominas , M.A. , and Wu , H.-T. : ‘ Decomposing non-stationary signals with time-varying wave-shape functions’ , IEEE Transactions on Signal Processing , 2021 , 69 , pp. 5094 – 5104 OpenUrl ↵ Nadakuditi , R.R. : ‘ Optshrink: An algorithm for improved low-rank signal matrix denoising by optimal, data-driven singular value shrinkage’ , IEEE Transactions on Information Theory , 2014 , 60 , ( 5 ), pp. 3002 – 3018 OpenUrl CrossRef ↵ Benaych-Georges , F. , and Nadakuditi , R.R. : ‘ The singular values and vectors of low rank perturbations of large rectangular random matrices’ , Journal of Multivariate Analysis , 2012 , 111 , pp. 120 – 135 OpenUrl ↵ Condos , R. , Norman , R.G. , Krishnasamy , I. , Peduzzi , N. , Goldring , R.M. , and Rapoport , D.M. : ‘ Flow limitation as a noninvasive assessment of residual upper-airway resistance during continuous positive airway pressure therapy of obstructive sleep apnea’ , American journal of respiratory and critical care medicine , 1994 , 150 , ( 2 ), pp. 475 – 480 OpenUrl PubMed Web of Science ↵ Mooney , A.M. , Abounasr , K.K. , Rapoport , D.M. , and Ayappa , I. : ‘ Relative prolongation of inspiratory time predicts high versus low resistance categorization of hypopneas’ , Journal of Clinical Sleep Medicine , 2012 , 8 , ( 2 ), pp. 177 – 185 OpenUrl PubMed ↵ Genta , P.R. , Sands , S.A. , Butler , J.P. , Loring , S.H. , Katz , E.S. , Demko , B.G. , Kezirian , E.J. , White , D.P. , and Wellman , A. : ‘ Airflow shape is associated with the pharyngeal structure causing OSA’ , Chest , 2017 , 152 , ( 3 ), pp. 537 – 546 OpenUrl CrossRef PubMed ↵ Younes , M. : ‘ Role of respiratory control mechanisms in the pathogenesis of obstructive sleep disorders’ , Journal of applied physiology , 2008 , 105 , ( 5 ), pp. 1389 – 1405 OpenUrl CrossRef PubMed Web of Science ↵ Wu , H.-T. : ‘ Instantaneous frequency and wave shape functions (I)’ , Applied and Computational Harmonic Analysis , 2013 , 35 , ( 2 ), pp. 181 – 199 OpenUrl ↵ Lin , Y.-T. , Malik , J. , and Wu , H.-T. : ‘ Wave-shape oscillatory model for nonstationary periodic time series analysis’ , Foundations of Data Science , 2021 , 3 , ( 2 ), pp. 99 – 131 OpenUrl ↵ Silverman , R. : ‘ Locally stationary random processes’ , IRE Transactions on Information Theory , 1957 , 3 , ( 3 ), pp. 182 – 187 OpenUrl ↵ Huang , W.K. , Chung , Y.-M. , Wang , Y.-B. , Mandel , J.E. , and Wu , H.-T. : ‘ Airflow recovery from thoracic and abdominal movements using synchrosqueezing transform and locally stationary Gaussian process regression’ , Computational Statistics & Data Analysis , 2022 , 174 , pp. 107384 OpenUrl ↵ Bates , J.H. , Schmalisch , G. , Filbrun , D. , and Stocks , J. : ‘ Tidal breath analysis for infant pulmonary function testing. ERS/ATS task force on standards for infant respiratory function testing. European respiratory society/American thoracic society’ , European Respiratory Journal , 2000 , 16 , ( 6 ), pp. 1180 – 1192 OpenUrl Abstract / FREE Full Text ↵ Schmidt , M. , Foitzik , B. , Wauer , R.R. , Winkler , F. , and Schmalisch , G. : ‘ Comparative investigations of algorithms for the detection of breaths in newborns with disturbed respiratory signals’ , Computers and biomedical research , 1998 , 31 , ( 6 ), pp. 413 – 425 OpenUrl ↵ Daubechies , I. , Lu , J. , and Wu , H.-T. : ‘ Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool’ , Applied and Computational Harmonic Analysis , 2011 , 30 , ( 2 ), pp. 243 – 261 OpenUrl ↵ Savitzky , A. , and Golay , M.J.E. : ‘ Smoothing and differentiation of data by simplified least squares procedures’ , Analytical chemistry , 1964 , 36 , ( 8 ), pp. 1627 – 1639 OpenUrl CrossRef ↵ Seppä , V.-P. , Hyonen , J. , and Viik , J. : ‘ A method for suppressing cardiogenic oscillations in impedance pneumography’ , Physiological measurement , 2011 , 32 , ( 3 ), pp. 337 OpenUrl CrossRef PubMed ↵ Burh , H.E. : ‘ The Inspiration-Expiration Ratio During Truth and Falsehood’ , Journal of Experimental Psychology , 1921 , 4 , ( 1 ), pp. 1 OpenUrl CrossRef ↵ Azarbarzin , A. , Sands , S.A. , Taranto-Montemurro , L. , Oliveira Marques , M.D. , Genta , P.R. , Edwards , B.A. , Butler , J. , White , D.P. , and Wellman , A. : ‘ Estimation of pharyngeal collapsibility during sleep by peak inspiratory airflow’ , Sleep , 2017 , 40 , ( 1 ), pp. zsw005 OpenUrl PubMed Parekh , A. , Tolbert , T.M. , Mooney , A.M. , Ramos-Cejudo , J. , Osorio , R.S. , Treml , M. , Herkenrath , S.-D. , Randerath , W.J. , Ayappa , I. , and Rapoport , D.M. : ‘ Endotyping sleep apnea one breath at a time: an automated approach for separating obstructive from central sleep-disordered breathing’ , American Journal of Respiratory and Critical Care Medicine , 2021 , 204 , ( 12 ), pp. 1452 – 1462 OpenUrl PubMed Parekh , A. , Kam , K. , Wickramaratne , S. , Tolbert , T.M. , Varga , A. , Osorio , R. , Andersen , M. , de Godoy , L.B.M. , Palombini , L.O. , and Tufik , S. : ‘ Ventilatory burden as a measure of obstructive sleep apnea severity is predictive of cardiovascular and all-cause mortality’ , American Journal of Respiratory and Critical Care Medicine , 2023 , 208 , ( 11 ), pp. 1216 – 1226 OpenUrl PubMed ↵ Staykov , E. , Mann , D.L. , Duce , B. , Kainulainen , S. , Leppänen , T. , Töyräs , J. , Azarbarzin , A. , Georgeson , T. , Sands , S.A. , and Terrill , P.I. : ‘ Increased flow limitation during sleep is associated with increased psychomotor vigilance task lapses in individuals with suspected OSA’ , Chest , 2024 , 165 , ( 4 ), pp. 990 – 1003 OpenUrl PubMed ↵ Ruppert , D. , and Wand , M.P. : ‘ Multivariate locally weighted least squares regression’ , The annals of statistics , 1994 , pp. 1346 – 1370 ↵ Kavalieris , L. , and Hannan , E.J. : ‘ Determining the number of terms in a trigonometric regression’ , Journal of time series analysis , 1994 , 15 , ( 6 ), pp. 613 – 625 OpenUrl ↵ Aittokallio , T. , Virkki , A. , and Polo , O. : ‘ Understanding sleep-disordered breathing through mathematical modelling’ , Sleep Medicine Reviews , 2009 , 13 , ( 5 ), pp. 333 – 343 OpenUrl PubMed ↵ Guerrero , G. , Le Rolle , V. , Loiodice , C. , Amblard , A. , Pepin , J.-L. , and Hernández , A. : ‘ Modeling patient-specific desaturation paherns in sleep apnea’ , IEEE Transactions on Biomedical Engineering , 2021 , 69 , ( 4 ), pp. 1502 – 1511 OpenUrl ↵ Buades , A. , Coll , B. , and Morel , J.M. : ‘ A non-local algorithm for image denoising’ , in Editor (Ed.)^(Eds.): ‘Book A non-local algorithm for image denoising’ ( Ieee, edn .), pp. 60 – 65 ↵ Wu , H.-T. , and Zhou , Z. : ‘ Frequency detection and change point estimation for time series of complex oscillation’ , Journal of the American Statistical Association , 2024 , 119 , ( 547 ), pp. 1945 – 1956 OpenUrl View the discussion thread. Back to top Previous Next Posted April 25, 2025. Download PDF Data/Code Email Thank you for your interest in spreading the word about medRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery Message Subject (Your Name) has forwarded a page to you from medRxiv Message Body (Your Name) thought you would like to see this page from the medRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery Hau-Tieng Wu , Thomas M Tolbert , David M Rapoport medRxiv 2025.04.24.25326393; doi: https://doi.org/10.1101/2025.04.24.25326393 Share This Article: Copy Citation Tools Optimal Shrinkage-aided Airflow Decomposition Algorithm (OSADA) and Cardiac Oscillation Recovery Hau-Tieng Wu , Thomas M Tolbert , David M Rapoport medRxiv 2025.04.24.25326393; doi: https://doi.org/10.1101/2025.04.24.25326393 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Respiratory Medicine Subject Areas All Articles Addiction Medicine (568) Allergy and Immunology (863) Anesthesia (300) Cardiovascular Medicine (4436) Dentistry and Oral Medicine (444) Dermatology (382) Emergency Medicine (608) Endocrinology (including Diabetes Mellitus and Metabolic Disease) (1509) Epidemiology (15229) Forensic Medicine (30) Gastroenterology (1124) Genetic and Genomic Medicine (6600) Geriatric Medicine (668) Health Economics (997) Health Informatics (4538) Health Policy (1368) Health Systems and Quality Improvement (1613) Hematology (542) HIV/AIDS (1264) Infectious Diseases (except HIV/AIDS) (15916) Intensive Care and Critical Care Medicine (1103) Medical Education (623) Medical Ethics (146) Nephrology (667) Neurology (6599) Nursing (346) Nutrition (998) Obstetrics and Gynecology (1144) Occupational and Environmental Health (957) Oncology (3333) Ophthalmology (974) Orthopedics (369) Otolaryngology (420) Pain Medicine (436) Palliative Medicine (130) Pathology (663) Pediatrics (1693) Pharmacology and Therapeutics (691) Primary Care Research (711) Psychiatry and Clinical Psychology (5447) Public and Global Health (9232) Radiology and Imaging (2198) Rehabilitation Medicine and Physical Therapy (1370) Respiratory Medicine (1196) Rheumatology (593) Sexual and Reproductive Health (712) Sports Medicine (530) Surgery (712) Toxicology (99) Transplantation (289) Urology (265) (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'a00f41a878fa73e4',t:'MTc3OTY1NjcwNQ=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();