Assess Sleep Stage by Modern Signal Processing Techniques

JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 1 Assess Sleep Stage by Modern Signal Processing T echniques Hau-tieng W u, Ronen T almon, Y u-Lun Lo* Abstract —In this paper , two modern adaptive signal processing techniques, Empirical Intrinsic Geometry and Synchrosqueezing transform, are applied to quantify different dynamical features of the respiratory and electroencephalographic signals. W e show that the proposed features ar e theoretically rigorously supported, as well as captur e the sleep inf ormation hidden inside the signals. The features ar e used as input to multiclass support vector machines with the radial basis function to automatically classify sleep stages. The effectiveness of the classification based on the proposed features is shown to be comparable to human expert classification – the proposed classification of awake, REM, N1, N2 and N3 sleeping stages based on the respiratory signal (resp. respiratory and EEG signals) has the overall accuracy 81 . 7% (resp. 89 . 3% ) in the relativ ely normal subject group. In addition, by examining the combination of the respiratory signal with the electroencephalographic signal, we conclude that the respiratory signal consists of ample sleep information, which supplements to the inf ormation stored in the electroencephalographic signal. Index T erms —Sleep Stage; Empirical Intrinsic Geometry; Synchrosqueezing transform; br eathing pattern variability I . I N T R O D U C T I O N I N human beings, sleep is a uni versal recurring dynam- ical and physiological activity , and the quality of sleep influences our daily liv es in di verse ways. Ho wever , it was not until recently that sleep became a brach of medicine and found its role in several seemingly unrelated clinical problems. Physiologically , it is divided into two broad stages: rapid eye movement (REM), and non-rapid eye movement (NREM) [1]. Normally , sleep proceeds in cycles in between REM and NREM. The NREM stage is further divided into shallow sleep (stage N1 and N2) and deep sleep (stage N3). In all procedures identifying sleep stages, we need a sleep scoring process with the help of polysomnography (PSG), which includes electroencephalography (EEG), electromyogram (EMG), and electrooculogram (EOG), etc. Among these physiological signals, EEG signals are the most concentrated ones since the clinically acceptable stage of the sleep is majorly determined by reading the recorded EEG based on the R&K criteria, which were standardized in 1968 by Allan Rechtschaffen and Anthony Kales [2] and further developed by the American Academy of Sleep Medicine on 2007 (AASM 2007) [3]. Howe ver , due to the H.-T . Wu is with the Department of Mathematics, Uni versity of T oronto, T oronto Ontario Canada. (email: hauwu@math.toronto.edu) R. T almon is with the Department of Electrical Engineering, T ech- nion - Israel Institute of T echnology , Haifa, 32000, Israel. (email: ro- nen@ee.technion.ac.il) *Correspondence: Y .-L. Lo is with the Department of thoracic medicine, Chang Gung Memorial Hospital, Chang Gung Univ ersity , School of Medicine, T aipei, T aiwan. (email: loyulun@hotmail.com) subjectiv e judgement and different training background, the agreement of manual sleep scoring among trained clinicians and professionals has been known to be limited [4], thereby motiv ating the development of an objecti ve and automatic scoring system. Based on these clinical findings, various features of the EEG signals hav e been proposed to study the sleep dynamics, for example, time domain summary statistics, spectral analysis, coherence, time-frequency analysis, entropy , to name but a few [5], [6], [7], [8]. Recently , a theoretically solid approach suitable to model the underlying dynamics of the brain activity and estimate the ev olutionary dynamics from recorded EEG signal was proposed in [9], [10], and had been successfully applied to predict the pre-seizure state from the intra-cranial EEG signals [11], [12]. Howe ver , it is well kno wn that sleep is a global and system- atic behavior not localized solely in the brain. For example, the muscular atonia and low amplitude EMG are related to the significant changes in the breathing pattern during normal sleep: during NREM sleep, especially stage N3, breathing is remarkably regular , while during REM sleep, breathing is irregular with sudden changes. The abov e physiological facts hint that the respiratory pattern of the recorded breathing signal during sleep might well reflect the sleep stage. There have been some reported studies of the sleep stage from analyzing the respiratory signal [13], [14], [15], [16], [17]. In [13] (resp. [14]), an av eraged respiratory rate over a fixed window is used to estimate the REM and NREM (resp. aw ake and sleep). In [15], a notch filter based instantaneous frequency estimator is applied to extract features to differentiate a wake, REM and NREM. In [16], [17], the adaptive harmonic model and a modern time-frequency analysis technique hav e been applied to further quantify the notion of respiratory dynamic. In particular , the instantaneous respiratory rate has been related to awak e, REM, shallow and deep sleep stages, with a rigorous mathematical foundation. The above-mentioned physiological patterns inside the EEG and the respiratory signals are actually outcomes of the in- tricate deformation of the underlying sleep dynamics, which we call intrinsic dynamical features of the sleep, that are not directly accessible to us. Although it is not an easy task to fully model or estimate the dynamical system underlying sleep, we might expect to benefit if we are able to quantify and integrate these hidden intrinsic dynamical features. In this paper , we propose to combine two modern adaptiv e signal processing techniques, Empirical Intrinsic Geometry (EIG) and Sync hrosqueezing transform (SST) , to estimate these intrinsic dynamical features of sleep guiding the observ ed EEG JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 2 and respiratory signals – we define an index, referred to as Sleep Index , to quantify these features. Then, by applying the suitable classifier algorithm, we sho w that the extracted features are highly correlated to the sleep stage determined by reading the EEG by the AASM 2007 criteria. Indeed, the proposed classification based on the respiratory signal (resp. respiratory and EEG signals) has the o verall accuracy 81 . 7% (resp. 89 . 3% ) in the relativ ely normal subject group, which is comparable to human expert classification. The article is organized in the following way . In Section II, we summarize the theoretical background of EIG and SST and the associated models. Then the Sleep Index is introduced in Section III. In Section IV, the proposed Sleep Index is applied to study the whole night sleep signals. W e conclude with discussion in Section V. I I . T W O A L G O R I T H M S – S Y N C H RO S Q U E E Z I N G T R A N S F O R M A N D E M P I R I C A L I N T R I N S I C G E O M E T RY The work presented in this paper is an application of the modern signal processing techniques to study the sleep dynamics. In particular , we will extract different features from the respiratory and EEG signals by the well studied EIG [9], [18] and SST [19], [17]. As such, the theoretical material will be presented in a compact, informal manner emphasizing on the intuitions. W e provide a formal and rigorous summary without proof of the details. Those interested in the proofs are encouraged to read the associated references. A. Adaptive Harmonic Model and Synchr osqueezing T rans- form The major characteristic pattern of the respiratory signal is that it is almost periodic. W e call the movement of air from the environment into the lungs inspiration and the movement of air in the opposite direction expir ation . An inspiration and an expiration constitute a respir atory cycle . Br eathing process is a physiological process consists of a sequential respiratory cycles. In this paper , we focus on the breathing process and call the time-varying volume occupying the lung space the physiological respir atory signal . This general observation leads us to the following phenomenological model for the respiratory signal R ( t ) (without noise): R ( t ) = A ( t ) s ( φ ( t )) , (1) where we shall call s ( · ) the wav e shape function; it is a 1 -periodic real function that satisfies some mild technical conditions. See [16] for the details. The respiratory signals recorded from dif ferent devices, like the airflow measuring device or the chest wall mov ement, shall be understood as observations of the respiratory system. Different observations lead to different shape functions. W e call the deriv ati ve φ 0 ( t ) of the function φ ( t ) the instantaneous frequency (IF) of the respiratory signal R ( t ) . W e require IF to be positiv e, but it does not required to be constant as long as the variations are slight from one period to the next, i.e. | φ 00 ( t ) | ≤ φ 0 ( t ) for all time t , where  is some small, pre-assigned positi ve number . Like wise, W e call A ( t ) the amplitude modulation (AM) of R ( t ) , which should be positive, but is allowed to vary slightly as well, i.e. | A 0 ( t ) | ≤ φ 0 ( t ) for all time t . W e refer the interested reader to [17] for the technical details and a further discussion of the well-definedness of the definition of AM and IF . Note that our treatment of the respiratory signal is purely phenomenological; that is, the parameters and indices we will deriv e from the signal will be based solely on these signals themselves, and not on explicit, quantitative models of the underlying mechanisms. Physiologically , the quantities φ 0 ( t ) , A ( t ) and s in the model (1) quantify the dynamics of the breathing process, which we refer to as phenomenological dynamical featur es . For example, one way to quantify the widely used notion breathing rate variability (BRV) [20], [21], [22] is considering the IF and AM [23]. Indeed, if φ 0 ( t 0 ) > φ 0 ( t 1 ) where t 1 6 = t 0 , we know that the subject breaths faster at time t 0 than at time t 1 . W e mention that this “fast-slo w” momentary behavior in the respiratory signal has been shown to be clinically informative and can be helpful in the ventilator weaning prediction [22], [23] and sleep stage estimation [16], [17]. Due to the inevitable measurement error and other outliers appearing inside the system, we model the recor ded r espira- tory signal as Y ( t ) = R ( t ) + σ ( t ) ξ ( t ) , (2) where ξ ( t ) is a stationary generalized random process and σ is a smooth function which varies slowly . Here ξ ( t ) models the noise and other outliers and σ models the possible non- stationary nature of the noise. A particular example for ξ frequently encountered in practice is the Gaussian white noise. W e refer the read ha ving interest to [17] for further information about this noise model and its mathematical details. T o estimate the phenomenological dynamical features of R ( t ) , φ 0 ( t ) and A ( t ) , from Y ( t ) , we apply the Synchr osqueez- ing transform (SST), which is a special reallocation technique [19], [17]. In a nutshell, we ev aluate any linear time-frequency analysis on the observation Y ( t ) , for example the short time Fourier transform or the continuous wav elet transform, and we take the phase information hidden inside the chosen linear time-frequency analysis into account to obtain a sharpened time-frequency representation, which is denoted as S Y ( t, ξ ) . In addition to capturing the phenomenological dynamical features of R ( t ) , the SST provides a sharper time-frequency repre- sentation compared with the other traditional time-frequency analyses. See Figure 1 for an example of the respiratory signal and its instantaneous frequency . W e refer the reader to [19], [17] for more details. B. Empirical Intrinsic Geometry and its underlying Mathe- matical Model In many real-world applications, the seeming complicated time series collected from the system is controlled by a relativ ely simple underlying process. In some situations, when the underlying ev olutionary process lies on a low-dimensional Riemannian manifold, it can be parameterized through a manifold learning framew ork, which was first introduced and JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 3 Fig. 1. The time-frequency representation of the respiratory signal determined by the Synchrosqueezing transform (SST) superimposed by the respiratory signal (the blue curv e). The dominant black curve sho wn in the time-frequency representation indicated by the red arro ws is the instantaneous frequency (IF) of the respiratory signal. It is clear that the subject breathes faster during the time indicated by the blue arrow than that indicated by the green arrow . This observation is captured by the IF indicated by the blue and green arrows. studied in [24]. The main idea in [24] 1 and its extensions to time-series [9], [10] is to bridge between data mining, and in particular manifold learning, and dynamical systems. The authors’ observation that the accessible measurements at hand do not necessarily conv ey the true essence of the system led to the dev elopment of a more generalized model, which separates between measurements and intrinsic hidden variables. One particular example for such a dynamical system is the respiratory signal recorded during sleep – we consider the model that the ev olutionary process gov erning the respiratory signals is restricted to a low-dimensional Riemannian mani- fold, which is fundamentally different from the phenomeno- logical model (2). This dependency is encoded using the state- space formalism and the model of the recorded respiratory signal (2) is extended as follows:  Y ( t ) = R θ ( t ) + σ ( t ) ξ ( t ) , [measurement equation] d θ i ( t ) = a i ( θ i ( t )) d t + d w i ( t ) , [state equation] (3) where θ ( t ) := ( θ 1 ( t ) , . . . , θ d ( t )) forms the inaccessible in- trinsic state at time t that governs the respiratory signal R θ ( t ) and e volv es in time with unkno wn drifts a i and independent standard Brownian motions w i , i = 1 , . . . , d . The idea that lies behind the model (3) is twofold. First, the measured signal Y ( t ) has typical (unknown) dynamics, modeled here by the state equation, which need to be taken into account and encoded in the desired features. Second, the accessible signal is viewed as a measurement of the neural system controlling the sleep cycle. While it can be ef fected by numerous factors relating to the measurement modality (e.g., measurements of airflow or chest movements), the used equipment (e.g., the type of sensors and their exact positions), and noise, the true intrinsic variable we hav e interest in is the intrinsic states controlling the respiratory signal (represented 1 In the original paper [24], the method was referred to as Nonlinear Independent Component Analysis. Ho wev er, for the sake of avoiding pos- sible confusion with Independent Component Analysis, the name Empirical Intrinsic Geometry (EIG) was adpated in [9], [10]. here by θ ( t ) ). Indeed, the notation R θ ( t ) implies that the res- piratory signal depends upon the “real” physiological variable θ in an unkno wn way , possibly through its amplitude A θ ( t ) , wa ve shape s θ ( · ) , or instantaneous frequency φ 0 θ ( t ) . In the above model (3), howe ver , due to noise and other nui- sance factors, the measurement Y ( t ) might be too redundant to faithfully describe the dependency of the the respiratory signal on the underlying state and its temporal ev olutionary . Thus, to improv e the underlying state observability , we introduce a high dimensional (possibly nonlinear) observer Φ to the measured signal [12], i.e., Z ( t ) = Φ( Y ( t )) ∈ R m , [observation equation] (4) where Φ is a map from the suitable scalar valued functional space to the R m -valued functional space and m ≥ 1 is an integer specified by the observer . W ith the sampled observation set Z := { Z ( t i ) } N i =1 , a natural question is ho w to estimate θ ( t ) , namely , the system intrinsic state and dynamics of interest. Such an analysis may complement the phenomenological dynamical features pro- vided by the SST . While the SST mainly carries instantaneous information, recovering the intrinsic state of the dynamical system θ ( t ) provides a characterization of coarser , slower dynamical changes of the shape and structure of the signal, especially when the observer Φ is implemented as a transform that relies on short time frames analysis. It was shown in [24] that if the observations Z ( t ) can be written as a regular deterministic function f : R d → R m of the samples of the underlying state, i.e., Z ( t ) = f ( θ ( t )) , then, by It ˆ o’ s formula, we hav e d Z j ( t ) = d X i =1  1 2 f j ii + a i f j i  d t + d X i =1 f j i d w i ( t ) (5) where f j i = ∂ f j /∂ θ i and f j ii = ∂ 2 f j /∂ θ 2 i . By a direct calculation, the cov ariance matrix C ( t ) ∈ R m × m of the observation at time t define by C j,k ( t ) := Cov ( d y j ( t ) , d y k ( t )) , (6) satisfies C ( t ) = J ( t ) J T ( t ) , where J = ∇ f ∈ R m × d is the Jacobian of f . This key result, along with the assumption that θ ( t ) is locally stationary ev olving much more slo wly than the observation scale so that the it stays closely on a low- dimensional manifold M embedded in R d , which is referred to as the intrinsic state manifold , as well as the assumption that f is stably inv ertible on its range, allow the authors in [24] to estimate the inaccessible state through the solution of an eigenv ector problem, which will be described later in this section. The main step leading to the solution theory is the fol- lowing estimation. Suppose θ ( t ) , θ ( τ ) ∈ M , Z ( t ) = f ( θ ( t )) and Z ( τ ) = f ( θ ( τ )) . By the T aylor expansion of f , up to the error term O ( k Z ( t ) − Z ( τ ) k 4 ) [24], we hav e: k θ ( t ) − θ ( τ ) k 2 R d = 1 2 ( Z ( t ) − Z ( τ )) T × [ C − 1 ( t ) + C − 1 ( τ )]( Z ( t ) − Z ( τ )) . (7) Note that in our example, the function f leading to the observation depends on the observ er Φ . As a result, with the estimated covariance matrix from the accessible collected JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 4 data Z , we can build a graph Laplacian associated with the intrinsic state manifold from the finite observations Z using the estimated Euclidean distance between the corresponding underlying samples θ ( t i ) (7). This graph Laplacian gives rise to re-parametrization of the intrinsic manifold through diffusion maps (DM) [25]. This re-parametrization procedure aiming to extract the intrinsic dynamics of the observation is referred to as Empirical Intrinsic Geometry (EIG). The remaining key question is the choice or design of a “proper” observer Φ in (4) to the system. In particular , in order to accommodate the ine vitable noise in real-world signals, estimates of the conditional probability density p ( Z | θ ) (e.g., histograms) were proposed as observers in [9]. The analysis relies on the following facts: (a) any measurement noise is translated to a linear transformation in the conditional densities domain, and (b) the distance (7), which is the Mahalanobis distance , is inv ariant under linear transformations. Indeed, these two facts allo w for the estimation of the distance between two nearby samples on the intrinsic manifold in adverse noisy conditions. Howe ver , estimating the conditional probability densities requires a large amount of data and often is not feasible. Unfortunately , standard representations based on the Fourier transform are also inadequate for respiratory signals. By linear approximation of the function φ ( t ) around a nearby sample at t 0 , the respiratory signal in (1) can be approximated by R ( t ) ≈ A ( t 0 ) s ( φ ( t 0 ) − t 0 + φ 0 ( t 0 ) t ) , (8) As a result, the modulus of the Fourier transform of R ( t ) around t 0 is approximated by | ˆ R ( t 0 , ω ) | ≈ | A ( t 0 ) ˆ s ( ω /φ 0 ( t 0 )) | , (9) where ω is the frequency , ˆ R ( t 0 , ω ) is the Fourier transforms of R ( t ) around t 0 and ˆ s ( ω ) is the Fourier transform of s ( t ) , respectively , assuming A ( t ) changes slowly with time. The approximation in (9) implies that ev en an almost linear function φ ( t ) (i.e., when the IF is φ 0 ( t ) ≈ 1 ) is translated to large deformations in the Fourier domain in high frequencies [26]. Consequently , if we take the short time Fourier transform as the observer , the observations exhibit instabilities, thereby leading to irregular f and a poor estimation of the Euclidean distance on the underlying intrinsic state manifold (7). T o overcome the instability of the Fourier representation, following [12], we use the scattering transform as an ob- server . The scattering transform has a low variance because it is based on first order moments of contractiv e operators, it linearizes deformations, and it can represent effecti vely intermittent behavior [26], [27]. The scattering transform is computed based on a cascade of w av elet transforms and nonlinear modulus operators [26]. Here, we briefly revie w the construction procedure of its first and second order levels, since they were empirically shown to provide a sufficient representation of the signals considered in this paper . Let ψ ( t ) be a complex wa velet, whose real and imag- inary parts are orthogonal and have the same L 2 norm. Let ψ j ( t ) denote the dilated wav elet, defined as ψ j ( t ) := 2 − j ψ (2 − j t ) , ∀ j ∈ Z . Let Φ s ( R θ ( t )) denote the observations computed by applying the first and second le vel scattering transform to the signal samples R θ ( t ) , which are giv en by Φ s ( R θ ( t )) = ( || R θ ( t ) ∗ ψ j 1 ( t ) | ∗ ψ j 2 ( t ) | ∗ w ( t ) ∀ ( j 1 , j 2 ) ∈ Z n , n ∈ { 1 , 2 } ) j 1 ,j 2 where w ( t ) is a smoothing window , i.e., a scaling function associated with the mother wa velet. The scattering transform has been shown to be an observer that is especially suitable for deformations and intermittencies [12]. In particular, it was shown that it is regular with respect to time deformations. Therefore, the application of the scattering transform to the respiratory signal is particularly suitable. Building on the generality of the described analysis, in this study , we use it to represent the EEG signals as well. As the respiratory signal, the EEG signal measures a physiological phenomenon (“the brain activity”), but, it is subject to noise, interferences, and nuisance factors. Likewise, it can be repre- sented using a state-space model, similar to (3), giv en by  X ( t ) = E ζ ( t ) + V ( t ) [measurement equation] d ζ i ( t ) = α i ( ζ i ( t )) d t + d u i ( t ) , [state equation] where E ζ ( t ) ∈ R m 0 and X ( t ) ∈ R m 0 are the clean and noisy EEG signals, V ( t ) is a measurement noise, and ζ ( t ) := ( ζ 1 ( t ) , . . . , ζ d 0 ( t )) denotes the inaccessible intrinsic state representing the brain acti vity that gov erns the EEG signal E ζ ( t ) and evolv es in time with unknown drifts α i and independent standard Brownian motions u i , i = 1 , . . . , d 0 . By applying EIG to the recorded EEG signals, we may reconstruct the intrinsic states ζ . W e remark that this approach w as applied to identify the pre-seizure state from intracranial EEG data [11], [12]. W e refer the interested reader to [9], [10], [12] for more technical details and references. Before closing this section, we summarize the construction of the graph Laplacian parametrization. In a nutshell, the main ingredient is integrating local similarities at different scales, which leads to a global description of the data set. Unlike linear methods such as principal component analysis (PCA), a graph Laplacian parametrization embodies nonlinear relation- ships among the variables. In addition to the mathematical analysis results [25], [28], [29], it has been shown to be robust to noise perturbation [30], [31] and it is computationally efficient. W e outline the algorithm here and refer the readers to these literatures for further theoretical details. T ake N multiv ariate measurement samples Z = { Z ( t i ) } N i =1 and build a complete graph with vertices Z . W e first build an affinity matrix (or adjacency matrix) W of size N × N . The affinity between a pair of samples is defined by a metric d in the following way: W ij = e − d 2 ( Z ( t i ) , Z ( t j ))  , for i, j = 1 , . . . , N , i 6 = j . (10) Note that according to the noise analysis in [31], when the signal to noise ratio is small, it is beneficial to set the diagonal terms of the affinity matrix to 0 . In the present work, following the analysis in [24], the metric we choose is the Mahalanobis distance (7). It is clear that the matrix W is symmetric. Note that theoretically (and practically) we can choose a more general kernel function, but we focus on the Gaussian kernel to simplify the exposition. Then we define the diagonal JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 5 degree/density matrix D of size n × n , consisting of the sum of rows of W : D ii = N X j =1 W ij , for i = 1 , . . . , N . Based on W and D , the graph Laplacian is defined by L := I − D − 1 W . Note that under the manifold assumption, D − 1 exists. Also note that D − 1 W can be viewed as a transition matrix of a Markov chain on the samples. Since L is similar to the sym- metric matrix I − D − 1 / 2 WD − 1 / 2 , it has a complete set of right eigen vectors ϕ 1 , ϕ 2 , . . . , ϕ N with corresponding eigenv alues 0 = λ 1 < λ 2 ≤ · · · ≤ λ N ≤ 1 , where ϕ 1 = (1 , 1 . . . , 1) T [25]. By the abov e construction, the eigen vectors ϕ 1 , . . . , ϕ N are vectors in R N . Through the eigen vectors, the measurement samples are mapped into R ˆ d via Z ( t i ) 7→ ( ϕ 2 ( t i ) , . . . , ϕ ˆ d ( t i )) , for i = 1 , . . . , N . (11) where ˆ d is an estimate of the dimension of the intrinsic state of the system and is usually ˆ d  N . Estimating the intrinsic dimension of the system d extends the scope of the paper and is empirically set according to the spectral gap in the decay of the eigen values, as will be described in Section III. In (11), we obtain a ˆ d -dimensional parameterization of the measurements. In particular, we view the j th coordinate of the parameterization of Z ( t i ) , i.e., ϕ j +1 ( t i ) , as the j th coordinate of the recov ered hidden intrinsic state θ j ( t i ) , which we view as the features associated with the sleep stage in this analysis. An illustration of the DM reparametrization process with the first 3 non-trivial eigen vectors is shown in Figure 2. I I I . M AT E R IA L A N D M E T H O D A. Data Collection Standard polysomnography was performed with at least 6 hours of sleep to confirm the presence or absence of OSA from the clinical subjects suspicious of sleep apnea in the sleep center at Chang Gung Memorial Hospital (CGMH), Linkou, T aoyuan, T aiwan. The institutional revie w board of the CGMH approved the study protocol (No. 101-4968A3) and the enrolled subjects provided written informed consent. Four channel EEG signals (C3A2, C4A1, O1A2 and O2A1), two channel EOG signals and chin EMG were recorded at the sampling rate 200 Hz for sleep staging. Chest and abdominal motions are recorded by the piezo-electric bands and airflow was measured using thermistors and nasal pressure, both at the sampling rate 100 Hz. All signals were acquired on the Alice 5 data acquisition system (Philips Respironics, Murrysville, P A). Apneas and hypopneas were defined using AASM 2007 guidelines [3], and the apnea-hyponea index (AHI) provided is the value determined during sleep. T ake the recorded EEG signals, denoted as E k , k = 1 , . . . , 4 , and the respiratory signal, denoted as R . Suppose the recording time period is T = [0 , T ] . W e divide T into contiguous subintervals T i of τ seconds long, i = 1 , . . . , N ; that is, T = ∪ N i =1 T i and T i ∩ T j = ∅ for all i 6 = j . W e call T i the i -th epoch . W e will extract p > 0 features out of the recorded respiratory and EEG signals for each epoch. B. F eatur es from the r espiratory signal Giv en a recorded respiratory signal R ( t ) , we extract its phenomenological dynamical featur es , including the instan- taneous frequency φ 0 ( t ) and the amplitude modulation A ( t ) by applying the SST . Denote the estimated instantaneous frequency by ˜ φ 0 ( t ) and the amplitude modulation by ˜ A ( t ) . The mean of ˜ A ( t ) restricted to the i -th epoch, denoted as the A i , and the mean of ˜ φ 0 ( t ) restricted to the i -th epoch, denoted as φ 0 i , form the first two features for the respiratory signal for the i -th subinterval. The third feature, denoted as v i , is obtained by e valuating the standard deviation of ˜ φ 0 ( t ) on the interval of length 30 seconds centered on the middle of the i -th epoch. W e apply the analysis described in Section II-B to R ( t ) in order to complement the phenomenological dynamical features and to obtain a characterization of the structural, slower underlying variables of the data. Here as well we obtain the graph Laplacian L ( R ) ∈ R N × N . Then, the eigen vectors and eigen values of L ( R ) are gi ven by L ( R ) ϕ ( R ) j = λ ( R ) j ϕ ( R ) j . The first ˆ d ( R ) ≥ 1 nontrivial eigenv ectors are chosen based on the following “spectral gap” thresholding criteria λ ( R ) ˆ d ( R ) +1 < δ and λ ( R ) ˆ d ( R ) +2 ≥ δ, (12) where 0 < δ < 1 is the threshold chosen by the user . Thus, using (11), we obtain ˆ d ( R ) intrinsic dynamical featur es of the respiratory system. C. F eatur es from the EEG signal Giv en the EEG signal E k ( t ) recorded from the k -th channel, we run the analysis described in Section II-B and obtain the graph Laplacian L ( E ,k ) ∈ R N × N . Then, the eigen vec- tors and eigen v alues of L ( E ,k ) are given by L ( E ,k ) ϕ ( E ,k ) j = λ ( E ,k ) j ϕ ( E ,k ) j with 0 = λ ( E ,k ) 1 ≤ λ ( E ,k ) 2 ≤ . . . . The first ˆ d ( E ,k ) ≥ 1 nontrivial eigen vectors are chosen based on the thresholding criteria (12) with the same δ . Using the eigen- vectors, each subinterv al of the EEG signal E k ( t ) is mapped into a sub-vector of ˆ d ( E ,k ) dimensions according to (11). By collecting the low dimensional vectors of all the channel, we obtain a vector consisting of P 4 k =1 ˆ d ( E ,k ) intrinsic dynamical featur es of the cortical activity for each subinterval. D. Sleep Inde x W e consider the following two feature vectors. The first one is extracted only from the respiratory signal and is referred as the Respiratory Index : r i :=  A i , φ 0 i , v i , ϕ ( R ) 2 ( i ) , . . . , ϕ ( R ) ˆ d ( R ) +1 ( i )  . The second one is extracted only from the EEG signals and is referred as the EEG Index : e i :=  ϕ ( E , 1) 2 ( i ) , . . . , ϕ ( E , 1) ˆ d ( E , 1) +1 ( i ) , . . . , ϕ ( E , 4) 2 ( i ) , ϕ ( E , 4) ˆ d ( E , 4) +1 ( i )  . An analysis result of the O1A2 EEG signal, denoted as  ϕ ( E , 1) 2 ( i ) , . . . , ϕ ( E , 1) 4 ( i )  , is shown in Figure 2. Clearly dif- ferent sleep stages represented in dif ferent colors ha ve dif ferent features and are well clustered. In addition, these different JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 6 Fig. 2. The intrinsic dynamical features of the cortical activity extracted from the O1A2 EEG signal by the scattering Empirical Intrinsic Geometry (EIG). On the left, the scattering EIG is illustrated – the graph Laplacian is built up from the Mahalanobis distance from the EEG signal via the scattering operator . On the right, the top three nontri vial eigen values of the graph Laplacian are used to sho w the underlying evolutionary dynamics. The blue circles (resp. c yanid circles, yellow circles and red circles) represent the awake (resp. REM, N1 and N2 and N3) sleep stage. It is clear that the extracted dynamical features well parametrize the sleep stages in the sense that different sleep stages are located in different places. sleep stages are organized in a continuous but nonlinear way – from the right hand side of the figure to the left hand side we hav e awake, REM, N1 and N2 and deep sleep stages. Next, the 3 phenomenological respiratory features, the ˆ d ( R ) intrinsic respiratory features and the intrinsic dynamical fea- tures of the cortical acti vity at the i -th epoch are combined together to comprise the Sleep Index with p = P 4 k =1 ˆ d ( E ,k ) + 3 + ˆ d ( R ) : s i :=  r i , e i  . E. Sleep Stag e Classifier Support vector machine (SVM) is a commonly used tech- nique for the purpose of classification in statistical learning theory [32]. In a nutshell, SVM determines a hyperplane in the space separating the data set into two disjoint subsets, such that each subset is lying in a different side of the hyperplane. W ith the help of the reproducing kernel Hilbert space theory , SVM is generalized to the kernel SVM , which allo ws for classi- fication with nonlinear relationship; that is, a nonlinear surface separating the data set into two disjoints subsets may be used. W e refer the interested reader to [32] for technical details. For the sake of identifying the (possible) nonlinear relationship between dif ferent sleep stages, in this work we choose the radial based function (RBF), K ( x , x 0 ) = exp( − k x − x 0 k 2 2 2 σ 2 ) , where σ > 0 , as the kernel function. Note that our dataset is multi-class – the response has more than 2 categories – therefore, we need to further generalize the kernel SVM to the multi-class SVM to complete our mission. T o this end, we apply the one versus all (O V A) classification scheme [33]. Despite its simplicity , the O V A classification scheme is highly effecti ve and useful, as was extensi vely sho wn and discussed in [33]. Group data will be reported as mean ± standard de viation unless otherwise specified. I V . R E S U L T T en subjects without sleep apnea (AHI less than 5 ) were chosen for this study . The demographic characteristics of the individuals whose data was used are as follows: 6 males and 4 females, age: 45 . 9 ± 12 . 3 years, range 28 − 61 years; BMI: 23 . 6 ± 1 . 9 kg/m 2 , range 21 . 5 − 28 kg/m 2 ; AHI: 1 . 9 ± 1 . 1 , range 0 . 4 − 3 . 4 . The total recorded time are of length 384 ± 27 . 8 minutes with range 363 − 443 minutes and we hav e a sleep period time of 367 ± 27 . 5 minutes with range 338 − 428 minutes for the sleep stage estimation. W e divide the whole night sleep into contiguous epochs of length 2 . 56 seconds. W e take δ = 0 . 01 and the dimension of the Sleep Index p is 11 . 2 ± 1 . 69 . W e consider the sleep stages in this study: R = { A wake , REM , N1 , N2 , N3 } =: { 1 , 2 , . . . , 5 } . Here to simplify the notation, we reindex the set of sleep stages and use the teletype-font to avoid confusion; that is, 1 is the aw ake stage, etc. Then we generate the different indices, { s i } N i =1 , { r i } N i =1 and { e i } N i =1 , from the recorded EEG and respiratory signals. The sleep stages in R are determined by the sleep expert as the ground truth. The O V A kernel SVM with the RBF kernel with σ = 1 is applied to classify the different sleep stages. Suppose there are n ` subintervals with sleep stage ` , where ` = 1 , . . . , 5 , in the validation dataset. Denote n i,j to be the number of subintervals with the sleep stage i as the gold standard, but classified as the sleep stage j . W e call the 5 × 5 matrix N with the ( i, j ) -th entry n i,j the confusion matrix . W e also define the confusion per centage matrix P as a 5 × 5 matrix with its ( i, j ) entry n i,j P 5 j =1 n i,j . W e will call P i,i the sensitivity (SE) of the sleep stage i prediction, which is denoted as SE( i ). W e will also report the ov erall accuracy (A C) denoted as A C := JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 7 P 5 i = 1 n i , i P 5 i,j = 1 n i , j and the specificity (SP) of the sleep stage i denoted as SP ( i ) := P j 6 = i n j , j P k P j 6 = i n j , k . Note that these definitions are direct generalizations of the A C, SE and SP of the binary categorical response data. T o prev ent ov er-fitting and confirm the classification result, we run the repeated random sub-sampling validation 25 times and ev aluate the av erage. T o be more precise, we randomly partition the data into the training dataset and the validation dataset – the training dataset comprises 80% of the features and the rest are used to form the validation dataset. The trained classifier based on the training dataset is applied to predict the sleep stages of the validation dataset. W ith the above preparation, first, we show that the proposed features capturing the sleep information hidden inside the res- piratory signal are not only theoretically rigorously supported, but also useful in practice. The overall A C is 81 . 7% . The error bar of SE and SP of correlating the Respiratory Index and the sleep stages R over 25 repeated random sub-sampling validation for the 10 subjects is shown in the light gray curve in Figure 3. The a verage SE’ s (resp. SP’ s) over 10 subjects for the awak e, REM, N1, N2 and N3 stages are 82% , 89% , 72% , 82% and 62% (resp. 81% , 81% , 83% , 82% and 82% ). Second, we show that the EEG Index also correlates with the sleep stages. The overall A C is 71 . 6% . The error bar of the SE and SP over 25 repeated random sub-sampling v alidation for the 10 subjects is shown in the dark gray curve in Figure 3. The av erage SE’ s (resp. SP’ s) ov er 10 subjects for the awake, REM, N1, N2 and N3 stages are 70% , 67% , 50% , 74% and 54% (resp. 70% , 71% , 73% , 65% and 71% ). Next, we combine all the features extracted from the respi- ratory signal and the EEG signals and show the result is better than simply using the Respiratory Indices or EEG Indices. The ov erall A C is 89 . 3% . The error bar of the SE and SP over 25 repeated random sub-sampling validation for the 10 subjects is shown in the black curve in Figure 3. The av erage SE’ s (resp. SP’ s) ov er 10 subjects for the aw ake, REM, N1, N2 and N3 stages are 85% , 94% , 79% , 90% and 68% (resp. 89% 89% , 90% , 87% and 90% ). W e then apply the Mann-Whitney U test to test if the Sleep Index better predicts sleep stage than the Respiratory Index under our setup. The p v alue less than 0 . 01 is considered significant. For the 25 realizations of sub-sampling validation from 10 subjects, we obtained 250 SE’ s and 250 SP’ s for different indices respectiv ely . The Mann-Whitney U test is applied to see if the SE’ s and SP’ s are significantly different. The performance of the Sleep Index compared with the Respiratory Index on the awake, REM, N1, N2 and N3 stages in the sense of SE (resp. SP) are all significant with p-values < 0 . 001 ( < 0 . 001 ). Lastly , to better present the classification result, the av eraged confusion percentage matrices over all subjects and sub- sampling realizations based on the Respiratory Index, EEG Index and the Sleep Index are shown in Figure 4. Note that the diagonal entries are the SE’ s of sleep stage prediction. V . D I S C U S SI O N The results in Section IV show that an accurate estimation of all sleep stages by solely analyzing the respiratory signal is possible by combining EIG and SST . Indeed, in addition to the ov erall AC 81 . 7% , the average SE is greater than 72% except N3, and the av erage SP is greater than 81% . On the other hand, we mention that while the features of the respiratory signal extracted by EIG and SST are complementary , only EIG can be applied to the EEG signal, since the EEG signal can not be modeled by the adaptiv e harmonic model. The ov erall A C based on EIG applied to the EEG signals is 71 . 8% , the SE is greater than 67% , except N1 and N3, and the av erage SP is greater than 65% . Namely , the performance of the EEG Index is not better than the Respiratory Index. Nev ertheless, we see an impro vement in the classification performance based on the Sleep Index, which contains information from both the respiratory and EEG signals. The overall A C is increased to 89 . 3% , the average SE is now greater than 79% except N3, and the average SP is greater than 87% . In addition, it has been shown that the SE and SP of the Sleep Index are significantly better than those of the Respiratory Index. Moreov er, the confusion percentage matrices also indicate that except N3, the mis-classification does not land in any specific sleep stage. The above findings lead to the following two tentati ve conclusions: 1. in addition to the EEG signals, the respiratory signal contains ample information about the sleep stage; 2. combining the rele v ant but dif ferent information hidden inside the respiratory and the EEG signals leads to a better result. The main innov ation in our sleep depth analysis is the com- bination of the clinical observation and modern adaptive signal processing techniques. From the clinical standpoint, we take the well known physiological fact that in addition to the brain activity , sleep is a global dynamical process in volving dif ferent sub-system dynamics, in particular the significant changes in the respiratory pattern among different sleep stages. From the signal processing standpoint, we emphasize the importance of the nonlinearity controlling the sleep cycle and focus on finding suitable mathematical tools not only adaptiv e to the signal but also with sufficient rigorousness to quantify the clinical observation. Indeed, since the unaccessible intrinsic sleep dynamics is reflected in the nonlinear behavior of the respiration, and the two modern signal processing techniques, EIG and SST , hav e being theoretically studied to well quantify these nonlinearity , we obtain effecti ve features by analyzing the recorded respiratory signal, which surrogate the intrinsic sleep dynamics. The meaning of accuracy deserves some discussion. It is well known that the sleep stage determination agreement between dif ferent sleep experts is limited to 85% e ven when the subjects under e xamination are normal, and it is e ven worse on the abnormal subjects [4] 2 . Although our cases are not diagnosed as sleep apnea, they cannot be considered as in the normal population either , thereby attaining accuracy rates 2 It is reported in [4] that the mean agreement in the normal subset is higher (mean 76%, range 65-85%) than in the subset of sleep disordered breathing (mean 71%, range 65-78%). JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 8 Fig. 3. The error bar of the performance of each features for predicting the sleep stage. The upper (resp. lo wer) subfigure is the sensitivity (resp. specificity) of predicting different sleep stages by dif ferent indices over 25 repeated random sub-sampling v alidation. The black (resp. light gray and dark gray) curve is for the Sleep Index (resp. Respiratory Index and EEG Index). The x-axis is the subject index ranging from 1 to 10 . Fig. 4. The averaged confusion percentage matrices over all subjects and sub-sampling realizations based on the Respiratory Index (resp. EEG Index and Sleep Index) is shown in the left (resp. middle and right) subfigure. The percentage is represented by the color . The darker the entry is, the higher the value is. The precise v alue is shown in the color bar on the top of each matrix. Here, 1 (resp. 2 , 3 , 4 and 5 ) in the x- and y-axis tick label stands for awake (resp. REM, N1, N2 and N3). It is clear to see the inclination of mis-classifying N3 into N2. higher than 80% in our subjects may not be meaningful. On the other hand, we found that the classification of N3 stage is consistently worse and its mis-classification tends to land in N2, as is sho wn in Figure 4. Notice that the subjects in our study are on average 48 years old, and the distribution of N3 sleep stages in the normal population of this age is 4 − 20% . Howe ver , the N3 sleep stages in our study cases is 3 . 1% ± 3 . 26% with 25% and 75% quantiles 0 . 8% and 5% respectiv ely , which is much fewer than those in the normal population. Since the number of N3 in the training set is relativ ely small, ev en by applying the weighted SVM to handle the unbalanced data, we do not expect to attain a compatible classification rate of N3. This unbalanced training set issue, combined with the stable breathing pattern during N2 and N3, might explain the inclination of mis-classifying N3 into N2. Furthermore, while the accuracy of our classification is compatible with/better than the state-or-art reported results, we are able better classify between different sleep stages. Indeed, in [14], the ov erall accuracy of classifying awak e and sleep is 83 . 6% based on the respiratory signal; in [13] an av eraged respiratory rate classifies REM and NREM with accuracy over 85% ; in [15], a notch filter based IF estimator is applied to extract respiratory features, which classifies awak e, REM and NREM with mean accuracy approximately 70% ; in [17], the IF estimated by SST is shown to be able to distinguish awak e, REM, shallow and N3 with statistical significance. W ith the above discussions, we conclude that our features and the selected classifier are accurate. The sleep depth estimation by the EEG Index is inferior with respect to the traditional EEG analysis. T o understand this result, we briefly revisit how a sleep expert determines the sleep stage. Based on the protocol criteria, in addition to an EEG signal of duration that exceeds 30 seconds, the expert also takes into account past and future EEG signals to determine the sleep stage. Howe ver , in our study , the EEG Index is based on the signal in epochs of length 2 . 56 seconds. The choice of 2 . 56 -second interval is for the sake of balancing between the dimension and number of data points in EIG. Although the local covariance structure of the EEG signal is taken into account in the EIG analysis, this relationship is different from the protocol criteria. As a result, we do not expect to obtain a compatible stratification power . Howe ver , we see that even if we only focus on these short-term EEG signals, we still can predict the sleep stage up to some accuracy and it does help to attain a better classification rate when combined with the Respiratory Index. This hints the possibility that some useful information is hidden inside a finer scale EEG signals. This interesting potential will be reported in the future study . The discussion would not be complete without mentioning the shortcomings of our study . First, we focus on a small database containing only 10 relativ ely normal subjects in this study . T o confirm the usefulness of the proposed features, JOURNAL OF L A T E X CLASS FILES, VOL. 11, NO. 4, DECEMBER 2012 9 we need to study a larger database with different types of subjects. Second, the chosen features, in particular the features selected by EIG, are subject-dependent. Indeed, different sub- jects might hav e different dynamical systems and the number of dominant components determined by EIG might vary . A theoretical and practical study of integrating the proposed features among different subjects is undergoing. In conclusion, by applying modern signal processing tech- niques to EEG and respiratory signals, we find a set of suitable features, which allow us to predict the sleep stages accurately . In addition to gaining insight into the dynamics controlling the sleep dynamics, the automatic annotation system based on the analysis might lead to an objectiv e classification as well as reduce the required human expert analysis in volv ed in sleep ev aluation. A C K N O W L E D G E M E N T S Hau-tieng W u and Ronen T almon thank the helpful dis- cussions with Professor Ronald Coifman. Ronen T almon ac- knowledges the support by the European Union’ s - Sev enth Framew ork Programme (FP7) under Marie Curie Grant No. 630657. Y u-Lun Lo ackno wledges the support by T aiwan National Science Council grant 101-2220-E-182A-001 and 102-2220-E-182A-001. R E F E R E N C E S [1] T . Lee-Chiong, Sleep Medicine: Essentials and Re view . Oxford, 2008. [2] A. Rechtschaf fen and A. Kales, A Manual of Standardized T erminology , T echniques and Scoring System for Sleep Stages of Human Subjects . W ashington: Public Health Service, US Government Printing Office, 1968. [3] C. Iber, S. Ancoli-Isreal, A. Chesson Jr ., and S. Quan, The AASM Manual for Scoring of Sleep and Associated Events-Rules: T erminology and T echnical Specification . American Academy of Sleep Medicine, 2007. [4] R. Norman, I. Pal, C. Ste wart, J. W alsleben, and D. Rapoport, “Breathing pattern in humans: diversity and indi viduality , ” Inter observer agreement among sleep scorers from differ ent centers in a large dataset , vol. 23, no. 7, pp. 901–908, 2000. [5] V . Bajaj and R. B. Pachori, “ Automatic classification of sleep stages based on the time-frequency image of { EEG } signals, ” Computer Methods and Pr ograms in Biomedicine , vol. 112, no. 3, pp. 320 – 328, 2013. [6] N. Kannathal, M. Choo, U. Acharya, and P . Sadasivan, “Entropies for detection of epilepsy in eeg, ” Computer Methods and Pro grams in Biomedicine , vol. 80, pp. 187–194, 2005. [7] S. Blanco, R. Quiroga, O. Rosso, and S. Kochen, “Time-frequency analysis of electroencephalogram series, ” Physical Review E , vol. 51, no. 3, pp. 2624–2631, 1995. [8] S. Geng, W . Zhou, Q. Y uan, D. Cai, and Y . Zeng, “Eeg non-linear feature extraction using correlation dimension and hurst e xponent, ” Neur ological Resear ch , vol. 33, no. 9, pp. 908–912, 2011. [9] R. T almon and R. Coifman, “Empirical intrinsic geometry for nonlinear modeling and time series filtering, ” Pr oc. Nat. Acad. Sci. , vol. 110, no. 31, pp. 12 535–12 540, 2013. [10] R. T almon and R. R. Coifman, “Intrinsic modeling of stochastic dynam- ical systems using empirical geometry , ” Appl. Comput. Harmon. Anal. , 2014, T ech. Report Y ALEU/DCS/TR1467. [11] D. Duncan, R. T almon, H. P . Zaveri, and R. R. Coifman, “Identifying preseizure state in intracranial eeg data using diffusion kernels, ” Math- ematical Biosciences and Engineering , vol. 10, pp. 579–590, 2013. [12] R. T almon, S. Mallat, H. Zaveri, and R. R. Coifman, “Manifold learning for latent variable inference in dynamical systems, ” submitted , 2014, T ech. Report Y ALEU/DCS/TR1491. [13] G. S. Chung, B. H. Choi, K. K. Kim, Y . G. Lim, J. W . Choi, D.- U. Jeong, and K.-S. Park, “Rem sleep classification with respiration rates, ” in Information T echnology Applications in Biomedicine, 2007. IT AB 2007. 6th International Special T opic Conference on , 2007, pp. 194–197. [14] G. Guerrero-Mora, P . Elvia, A. Bianchi, J. Kortelainen, M. T enhunen, S. Himanen, M. Mendez, E. Arce-Santana, and O. Gutierrez-Navarro, “Sleep-wake detection based on respiratory signal acquired through a pressure bed sensor, ” in Engineering in Medicine and Biology Society (EMBC), 2012 Annual International Conference of the IEEE , 2012, pp. 3452–3455. [15] J. Sloboda and M. Das, “ A simple sleep stage identification technique for incorporation in inexpensive electronic sleep screening devices, ” in Aer ospace and Electr onics Confer ence (NAECON), Proceedings of the 2011 IEEE National , 2011, pp. 21–24. [16] H.-T . Wu, “Instantaneous frequenc y and wav e shape functions (I), ” Appl. Comput. Harmon. Anal. , vol. 35, pp. 181–199, 2013. [17] Y .-C. Chen, M.-Y . Cheng, and H.-T . W u, “Nonparametric and adaptiv e modeling of dynamic seasonality and trend with heteroscedastic and dependent errors, ” J. Roy . Stat. Soc. B , vol. 76, pp. 651–682, 2014. [18] R. T almon, I. Cohen, S. Gannot, and R. R. Coifman, “Diffusion Maps for Signal Processing: A Deeper Look at Manifold-Learning T echniques Based on Kernels and Graphs, ” IEEE T rans. Signal Process. , vol. 30, no. 4, pp. 75–86, 2013. [19] I. Daubechies, J. Lu, and H.-T . Wu, “Synchrosqueezed W avelet Trans- forms: an empirical mode decomposition-like tool, ” Appl. Comput. Harmon. Anal. , pp. 243–261, 2011. [20] M. Engoren, “ Approximate entropy of respiratory rate and tidal volume during weaning from mechanical ventilation, ” Critical Care Medicine , vol. 26, pp. 1817–1823, 1998. [21] G. Benchetrit, “Breathing pattern in humans: di versity and individuality , ” Respir . Physiol. , vol. 122, no. 2-3, pp. 123 – 129, 2000. [22] M. W ysocki, C. Cracco, A. T eixeira, A. Mercat, J. Diehl, Y . Lefort, J. Derenne, and T . Similowski, “Reduced breathing v ariability as a pre- dictor of unsuccessful patient separation from mechanical ventilation, ” Critical Care Medicine , vol. 34, pp. 2076–2083, 2006. [23] H.-T . Wu, S.-S. Hseu, M.-Y . Bien, Y . R. Kou, and I. Daubechies, “Evaluating physiological dynamics via synchrosqueezing: Prediction of ventilator weaning, ” IEEE T ransactions on Biomedical Engineering , vol. 61, pp. 736–744, 2013. [24] A. Singer and R. R. Coifman, “Non-linear independent component analysis with diffusion maps, ” Appl. Comput. Harmon. Anal. , vol. 25, no. 2, pp. 226 – 239, 2008. [25] R. R. Coifman and S. Lafon, “Diffusion maps, ” Appl. Comput. Harmon. Anal. , vol. 21, no. 1, pp. 5–30, 2006. [26] S. Mallat, “Group in variant scattering, ” Pure and Applied Mathematics , vol. 10, no. 65, pp. 1331–1398, 2012. [27] J. Bruna, E. Mallat, S. Bacry , and J. F . Muzy , “Multiscale intermittent process analysis by scattering, ” submitted, , 2013. [28] M. Belkin and P . Niyogi, “Con ver gence of laplacian eigenmaps, ” in Advances in Neural Information Pr ocessing Systems 19: Proceedings of the 2006 Conference , vol. 19. The MIT Press, 2007, p. 129. [29] A. Singer and H.-T . W u, “Spectral con ver gence of the connection laplacian from random samples, ” submitted , 2013. [30] N. El Karoui, “On information plus noise kernel random matrices, ” The Annals of Statistics , vol. 38, no. 5, pp. 3191–3216, 2010. [31] N. El Karoui and H.-T . W u, “Connection graph Laplacian methods can be made robust to noise, ” submitted , 2014. [32] B. Scholkopf and A. Smola, Learning with K ernels . MIT Press, 2002. [33] R. Rifkin and A. Klautau, “In Defense of One-Vs-All Classification, ” Journal of Machine Learning Researc h , vol. 5, pp. 101–141, 2004.