Estimation of Heart Rate and Respiratory Rate from PPG Signal Using Complementary Ensemble Empirical Mode Decomposition with both Independent Component Analysis and Non-Negative Matrix Factorization

This paper proposes a framework combining the complementary ensemble empirical mode decomposition with both the independent component analysis and the non-negative matrix factorization for estimating both the heart rate and the respiratory rate from the photoplethysmography (PPG) signal. After performing the complementary ensemble empirical mode decomposition on the PPG signal, a finite number of intrinsic mode functions are obtained. Then, these intrinsic mode functions are divided into two groups to perform the further analysis via both the independent component analysis and the non-negative matrix factorization. The surrogate cardiac signal related to the heart activity and another surrogate respiratory signal related to the respiratory activity are reconstructed to estimate the heart rate and the respiratory rate, respectively. Finally, different records of signals acquired from the Medical Information Mart for Intensive Care database downloaded from the Physionet Automated Teller Machine (ATM) data bank are employed for demonstrating the outperformance of our proposed method. The results show that our proposed method outperforms both the digital filtering approach and the conventional empirical mode decomposition based methods in terms of reconstructing both the surrogate cardiac signal and the respiratory signal from the PPG signal as well as both achieving the higher accuracy and the higher reliability for estimating both the heart rate and the respiratory rate.


Introduction
With increasing the life pressure, the cardiorespiratory diseases [1] became the major death reasons of the humans. Thus, it is necessary for the general public to monitor the cardiorespiratory activities [2] so that any abnormal heart situation and any abnormal respiration situation such as the acute physiologic deterioration [3], the cardiovascular diseases [4] and the long term cardiovascular related illnesses [5] can be detected earlier. In addition, both the postoperative treatment [6] and the rehabilitation management [7] can be performed in the early stage. It is worth noting that both the heart rate and the respiratory rate are the important parameters for representing the health conditions. As these signals can be estimated via many consumer electronic devices [8], the general public can monitor the cardiorespiratory activities via these signals with the continuous, noninvasive and comfortable means. However, both the commonly used electrocardiogram (ECG) based heart beat monitoring [9] means and the nasal thermistor based respiratory activity monitoring [10] means are uncomfortable for the patient to use.
The PPG signal is a commonly used signal for measuring the oxygen saturation in the blood. The PPG signal composes of different components. These components are mainly modulated by the heart activities, the respiration activities and other physiological activities [11]. Hence, they are synchronous with both the cardiac rhythm and the respiratory rhythm [12]. By analyzing the PPG signal, the information of the cardiorespiratory activities such as both the heart rate and the respiratory rate could be estimated. Monitoring the cardiorespiratory activities via the PPG signal is a well established noninvasive technique. Comparing with the ECG technique, the hardware implementation cost is much lower. Therefore, the PPG technique is not only used in both the anesthesia and the intensive care units in the hospital, but it is also implemented in a wearable device and used by the general public for monitoring the health condition.
To estimate both the heart rate and the respiratory rate using the PPG signal, the simple digital filtering techniques have been employed [12][13][14]. However, the performance is highly dependent on the cutoff frequency of the filter. Nevertheless, it was empirically selected. Although there are some analytical methods, these methods are very sensitive to the type of the noise. Hence, the result is very poor if the PPG signal is corrupted by a motion artifact. To address this issue, the estimation method based on the time frequency analysis approach [15,16] such as that based on the wavelet transform method [17][18][19][20] was proposed. For example, both the heart rate and the respiratory rate are estimated by the pulse oximeter using the PPG signal. Although this method is less sensitive to both the noise and the motion artifact, these methods require to select more than one parameter such as both the mother wavelet function and the total number of the decomposition level in the filter bank. In practice, these parameters are unknown.
In the recent years, the empirical mode decomposition was proposed [21] as an adaptive tool for processing both the nonstationary and the nonlinear signals. Many signals in the practice can be represented as the sums of their intrinsic mode functions and their residues. These components are localized in the frequency bands with their center frequencies sorted according to their indices. As a result of this nice property, several empirical mode decomposition based methods [22][23][24][25] were developed to decompose the PPG signals for estimating both the heart rates and the respiratory rates. However, due to both the intermittency and the noises corrupted to the PPG signals, the empirical mode decomposition based method may suffer from the occurrence of the mode mixing phenomenon. To address this problem, the ensemble empirical mode decomposition based method is proposed [26]. It is a noise assisted data analysis algorithm which can avoid the occurrence of the mode mixing phenomenon for processing the PPG signals. By applying the ensemble empirical mode decomposition to both the ECG and the PPG signal as well as employing the fusion approach [27], the respiratory rate was estimated using both the second intrinsic mode function and the third intrinsic mode function. However, this method is not robust in terms of applying to other ECG and PPG signals obtained from other subjects or devices. The joint ensemble empirical mode decomposition and the principal component analysis based method was proposed [28] to improve the robustness for estimating both the heart rate and the respiratory rate from the PPG signal by reducing the mode mixing effect. However, it may still yield an inaccurate result if the intrinsic mode functions unrelated to both the heart activity and the respiratory activity are taken into an account.
The complementary ensemble empirical mode decomposition has recently been proposed to avoid the occurrence of the mode mixing phenomenon. Meanwhile, it is found that the reconstruction error can be significantly suppressed [29]. On the other hand, both the independent component analysis [30,31] and the non-negative matrix factorization [32,33] are the powerful techniques for performing the blind source separation. The joint empirical mode decomposition and the independent component analysis based method have been proposed [34] to perform the source separation for the biomedical signal such as the electroencephalogram (EEG) signal. For extracting the fetal ECG signal from a single channel data, the joint empirical mode decomposition and the non-negative matrix factorization based method has also been proposed [35].
The novelty of this paper is to apply the complementary ensemble empirical mode decomposition, the independent component analysis and the non-negative matrix factorization to separate the PPG into two sets of signals as well as the principal component analysis to fuse these two sets of signals to generate two surrogate signals. Here, the surrogate signals refer to the signal components that are used to calculate the corresponding activities. In this paper, these two surrogate signals are the surrogate cardiac signal and the surrogate respiratory signal that are used to calculate the heart rate and the respiratory rate, respectively. It is worth noting that it is not required to select any parameter in the proposed method. Hence, the proposed method is adaptive. The computer numerical simulation results show that our proposed method could achieve the better results in terms of achieving the higher accuracies of both the estimated heart rate and the estimated respiratory rate as well as the more reliable results in terms of achieving the lower variances of the accuracies of both the estimated heart rate and the estimated respiratory rate. These improvements are important because both the accurate and the reliable heart rate as well as both the accurate and the reliable respiratory rate are essential and critical for the medical diagnosis of the cardiorespiratory diseases.
The outline of this paper is as follows. The existing methods for the analysis used in our proposed method are reviewed in Section 2.1. Then, our proposed method for the estimation of both the heart rate and the respiratory rate is presented in Section 2.2. Next, the computer numerical simulation results are shown in Section 3. Finally, a conclusion is drawn in Section 4.

Complementary Ensemble Empirical Mode Decomposition
The complementary ensemble empirical mode decomposition [29] is an improved version of the ensemble empirical mode decomposition. The main advantage of this method is to reduce the reconstruction error. Unlike the conventional ensemble empirical mode decomposition [26], a pair of uniformly distributed white noises with one positive valued and one negative valued is added to the original signal. It has been shown that the reconstruction error can be significantly suppressed. The detail procedures are described as follows.
Step 1 Let N be the total number of the white noises. Denote a set of positive valued white noise sequences as n i (t) and the corresponding negative valued white noise sequences as −n i (t) for i = 1, 2, . . . , N. Let S(t) be the original signal. Let a be the gain multiplied to the white noises. Add an i (t) and −an i (t) to S(t) for i = 1, 2, . . . , N. Denote the sequences corrupted by an i (t) and −an i (t) as r + i (t) and r − i (t) for i = 1, 2, . . . , N, respectively. That is: and Step 2 Decompose both r + i (t) and r − i (t) using the empirical mode decomposition. Let I + ij (t) and I − ij (t) be the j th intrinsic mode function of r + i (t) and r − i (t), respectively. Step 3 Let I j (t) be the j th intrinsic mode function of the reconstructed signal. Here, I j (t) is reconstructed by averaging I + ij (t) and I − ij (t) together. That is: Step 4 Let p be the total number of the intrinsic mode functions used for the reconstruction. Let r(t) be the residue based on the reconstruction of the signal using these p intrinsic mode functions. That is:

Independent Component Analysis
The independent component analysis is a statistical method and widely used in many signal processing applications such as in both the blind source separation application and the feature extraction application [30]. Let X be the observed mixed signal matrix, S be the instantaneously independent source and M be the mixing matrix. Let the dimension of X be m × n. Let both the total number of independent rows and the total number of the independent columns of both M and S be rank. Then, the dimensions of M and S are m × rank and rank × n, respectively. The blind source separation problem is to separate S from X using M. That is: To solve this problem, a common algorithm called the FastICA was used. It employs the maximum negative entropy as the search direction to find S [31].

Non-Negative Matrix Factorization
The non-negative matrix factorization is a kind of linear representations of the non-negative data [32,33]. It is also applied in many signal processing applications such as in the blind source separation application. Denote the observed non-negative signal matrix as V. Let W be the mixed matrix and H be the source signal matrix. Let the dimension of V be m × n. Let both the total number of independent rows and the total number of the independent columns of both W and H be rank. Then, the dimensions of W and H are m × rank and rank × n, respectively. The non-negative matrix factorization problem is to approximate V as the product of W and H. That is: Let W i,j and H i,j be the elements in the ith row and the jth column of W and H, respectively. In this paper, the non-negative matrix factorization with the sparseness constraints [36] is considered. In particular, let ||.|| 0 be the zero norm operator. Here, it refers to the total number of the nonzero elements in the operand. Let ∈ be the specification on the maximum number of the nonzero elements in W. Then, the non-negative matrix factorization problem is formulated as the sparse constrained optimization problem such that the two norm error ||V − WH|| 2 is minimized subject to ||W|| 0 ≤ ∈ as well as both W i,j ≥ 0 and H i,j ≥ 0 for all i,j. It is found that a better performance can be achieved compared to other basic algorithms.

Our Proposed Framework for Extracting the Cardiorespiratory Activity from the PPG Signal
In this paper, we propose a framework for extracting the cardiorespiratory activity from the PPG signal. The framework is shown in Figure 1. It mainly combines the existing methods discussed in Section 2.1. The overall algorithm can be divided into five stages and illustrated as follows.  Figure 1. The proposed framework for extracting the cardiorespiratory activity from the PPG signal.

Decomposition of the PPG Signal Using the Complementary Ensemble Empirical Mode Decomposition
It is worth noting that both the estimated heart rate and the estimated respiratory rate are dependent on the sampling rate of the PPG signal. Thus, both the separation of the cardiorespiratory related signals and the extraction of the cardiorespiratory activities from the PPG signal are also dependent on the sampling rate of the PPG signal. In this paper, first a PPG signal is segmented into a finite number of pieces. Here, the duration of each piece of the PPG signal is 30 s. Then, each piece of the PPG signal is decomposed into a finite number of intrinsic mode functions using the complementary ensemble empirical mode decomposition. It is worth noting the total numbers of the intrinsic mode functions of different pieces of the PPG signal may be different. To address this difficulty, the peak frequency of each intrinsic mode function of each piece of the PPG signal is computed. The peak frequencies of the intrinsic mode functions of two consecutive pieces of the PPG signal are linked together by using the dynamical programming approach. The mean and the variance of the peak frequencies of the intrinsic mode functions of each link are computed. It is found that the intrinsic mode functions of the PPG signal in the fourth link and the seventh link are corresponding to the ECG signal and the respiratory signal, respectively. This is because the means of the peak frequencies of the intrinsic mode functions in these two links fall to the ranges of the heart rate and the respiratory rate, respectively. An example of a piece of a PPG signal, its intrinsic mode functions as well as both the reference ECG signal and the reference respiratory signal are shown in Figure 2. It can be seen from the figure that the fourth intrinsic mode function of this piece of the PPG signal is the intrinsic mode function that its frequency is the closest to that of the reference ECG signal. In addition, the seventh intrinsic mode function of this piece of the PPG signal is the intrinsic mode function that its frequency is the closest to that of the reference respiratory signal. In fact, the frequencies of all the intrinsic mode functions in the fourth link and the seventh link for this subject are the closest to those of the reference ECG signal and the reference respiratory signal, respectively. Additionally, similar results are found for all other subjects.

Decomposition of the PPG Signal Using the Complementary Ensemble Empirical Mode Decomposition
It is worth noting that both the estimated heart rate and the estimated respiratory rate are dependent on the sampling rate of the PPG signal. Thus, both the separation of the cardiorespiratory related signals and the extraction of the cardiorespiratory activities from the PPG signal are also dependent on the sampling rate of the PPG signal. In this paper, first a PPG signal is segmented into a finite number of pieces. Here, the duration of each piece of the PPG signal is 30 s. Then, each piece of the PPG signal is decomposed into a finite number of intrinsic mode functions using the complementary ensemble empirical mode decomposition. It is worth noting the total numbers of the intrinsic mode functions of different pieces of the PPG signal may be different. To address this difficulty, the peak frequency of each intrinsic mode function of each piece of the PPG signal is computed. The peak frequencies of the intrinsic mode functions of two consecutive pieces of the PPG signal are linked together by using the dynamical programming approach. The mean and the variance of the peak frequencies of the intrinsic mode functions of each link are computed. It is found that the intrinsic mode functions of the PPG signal in the fourth link and the seventh link are corresponding to the ECG signal and the respiratory signal, respectively. This is because the means of the peak frequencies of the intrinsic mode functions in these two links fall to the ranges of the heart rate and the respiratory rate, respectively. An example of a piece of a PPG signal, its intrinsic mode functions as well as both the reference ECG signal and the reference respiratory signal are shown in Figure 2. It can be seen from the figure that the fourth intrinsic mode function of this piece of the PPG signal is the intrinsic mode function that its frequency is the closest to that of the reference ECG signal. In addition, the seventh intrinsic mode function of this piece of the PPG signal is the intrinsic mode function that its frequency is the closest to that of the reference respiratory signal. In fact, the frequencies of all the intrinsic mode functions in the fourth link and the seventh link for this subject are the closest to those of the reference ECG signal and the reference respiratory signal, respectively. Additionally, similar results are found for all other subjects.

Filtering on the Intrinsic Mode Functions
The normal ranges of the heart rates and the respiratory rates for the young population (including both children between 2 and 18 years old and young adults) are between 45 and 145 beats per minute as well as between 8 and 45 breaths per minute [23,37], respectively. Therefore, this paper chooses the frequency band between 0.1 Hz and 2.55 Hz as the possible signal band. Then, the intrinsic mode functions with the dominating frequencies lying in the range between 0.75 Hz and 2.55 Hz are selected to form a group of the candidate cardiac intrinsic mode functions. On the other hand, the intrinsic mode functions with the dominating frequency lying in the range between 0.1 Hz and 0.75 Hz are selected as a group of the candidate respiratory modes. Hence, the intrinsic mode functions in these frequency bands are categorized into two groups denoted as the cardiac group and the respiratory group, respectively.

Performing both the Independent Component Analysis and the Non-Negative Matrix Factorization on the Intrinsic Mode Functions
To separate the source signals embedding in each group, the set of intrinsic mode functions in the same group are represented as a matrix X with m rows and n columns. Here, each column of X is an intrinsic mode function. For performing the non-negative matrix factorization, the observed matrix is required to be non-negative. Hence, by letting each column of X subtracted from its minimum as the column of a non-negative matrix V, the non-negative matrix is obtained. Now, the source signal separation problem is formulated as two optimization problems defined in Equations (5) and (6), respectively. Then, by applying the FastICA and the non-negative matrix factorization with the sparseness constraints to X and V, respectively, these two rank × n source signal matrices S and H as well as these two m × rank mixed matrices M and W are found. In this paper, rank is set to 2 to separate four source signals. Then, they are mapped to the range between -1 and 1. Finally, they are denoted as , , ℎ and ℎ . After the source signals , , ℎ and ℎ are obtained from each group, the maximal cross correlation coefficient (MCC) is employed to select a pair of source signals to form a surrogate signal. Let be the cross correlation function between the source signal and ℎ . Here, the MCC between and ℎ is defined as MCC = .
(7) Figure 2. The PPG signal, its intrinsic mode functions as well as both the reference ECG signal and the reference respiratory signal.

Filtering on the Intrinsic Mode Functions
The normal ranges of the heart rates and the respiratory rates for the young population (including both children between 2 and 18 years old and young adults) are between 45 and 145 beats per minute as well as between 8 and 45 breaths per minute [23,37], respectively. Therefore, this paper chooses the frequency band between 0.1 Hz and 2.55 Hz as the possible signal band. Then, the intrinsic mode functions with the dominating frequencies lying in the range between 0.75 Hz and 2.55 Hz are selected to form a group of the candidate cardiac intrinsic mode functions. On the other hand, the intrinsic mode functions with the dominating frequency lying in the range between 0.1 Hz and 0.75 Hz are selected as a group of the candidate respiratory modes. Hence, the intrinsic mode functions in these frequency bands are categorized into two groups denoted as the cardiac group and the respiratory group, respectively.

Performing both the Independent Component Analysis and the Non-Negative Matrix Factorization on the Intrinsic Mode Functions
To separate the source signals embedding in each group, the set of intrinsic mode functions in the same group are represented as a matrix X with m rows and n columns. Here, each column of X is an intrinsic mode function. For performing the non-negative matrix factorization, the observed matrix is required to be non-negative. Hence, by letting each column of X subtracted from its minimum as the column of a non-negative matrix V, the non-negative matrix is obtained. Now, the source signal separation problem is formulated as two optimization problems defined in Equations (5) and (6), respectively. Then, by applying the FastICA and the non-negative matrix factorization with the sparseness constraints to X and V, respectively, these two rank × n source signal matrices S and H as well as these two m × rank mixed matrices M and W are found. In this paper, rank is set to 2 to separate four source signals. Then, they are mapped to the range between −1 and 1. Finally, they are denoted as s 1 , s 2 , h 1 and h 2 . After the source signals s 1 , s 2 , h 1 and h 2 are obtained from each group, the maximal cross correlation coefficient (MCC) is employed to select a pair of source signals to form a surrogate signal. Let R ij be the cross correlation function between the source signal s i and h j . Here, the MCC between s i and h j is defined as Obviously, the MCC can clearly indicate the similarity between the source signals obtained by different algorithms without affected by the phase difference between these source signals. Therefore, this paper selects a pair of source signal s i and h j corresponding to the maximal MCC as the surrogate signals.

Performing the Principal Component Analysis on the Selected Pairs of Source Signals
In order to both fuse a pair of the selected source signals and retain most of their variations, the first principal components obtained by applying the principal component analysis on the pairs of the source signals from a cardiac group and the respiratory group were used as a surrogate cardiac signal and the respiratory signal, respectively. Therefore, the obtained surrogate cardiac signal and the respiratory signal indicate the cardiac activity and the respiratory activity, respectively.

Estimation of Both the Heart Rate and the Respiratory Rate Using the Surrogate Signals
In order to extract the cardiorespiratory activities such as both the heart rate and the respiratory rate from the corresponding surrogate signals, the FFTs of these surrogate signals are computed.
Let f HR and f RR be the heart rate frequency and the respiratory rate frequency, respectively. They are defined as the peak frequencies of the surrogate heart signal and the surrogate respiratory signal, respectively. Let HR and RR be the estimated heart rate and the estimated respiratory rate, respectively. That is: HR = f HR * 60 (beats/min) (8) and RR = f RR * 60 (breaths/min), (9) respectively.

Database
The Medical Information Mart for Intensive Care database [38] contains 720 complete records of 90 different intensive care unit patients. Every record contains the signals such as the ECG signals, the PPG signals and the respiratory signals with 10 min duration. These signals are acquired simultaneously and sampled at 125 Hz. In this paper, all these 720 records are used to estimate both the heart rate and the respiratory rate. In addition, the obtained results are compared to the reference manually computed heart rate and the reference manually computed respiratory rate, respectively. However, due to the page limits, only four records of four subjects with the subject identification numbers 055m, 212m, 220m and 408m were demonstrated below. These four records are demonstrated because they correspond to different diseases. In particular, 055m is suffered from the respiratory failure, 212m is suffered from the pulmonary edema, 220m is suffered from the brain injury and 408m is suffered from the bleeding. In fact, other records also exhibit the similar results. Figure 3 shows both the surrogate cardiac signal and the surrogate respiratory signal obtained via both the empirical mode decomposition based method [23] and our proposed method. In addition, both the reference ECG signal and the reference respiratory signal are shown in Figure 3. From the figure, it can be seen that the rhythms of both the cardiac activity and the respiratory activity obtained by our proposed method are more obvious compared to those obtained by the empirical mode decomposition based method. This is particularly obvious for the surrogate respiratory signal. In addition, the artifact effect such as the boundary effect obtained by our proposed method is less obvious than that obtained by the empirical mode decomposition based method.

Performances
Sensors 2020, 20, x 8 of 13 proposed method is less obvious than that obtained by the empirical mode decomposition based method. To quantitatively evaluate the performance of our proposed method, both our proposed method and the empirical mode decomposition based method are applied to the signals with the lengths equal to 30 s durations. Here, the signals with only 30 s durations are used for evaluating the numerical performances. This is because it can save the required computational powers. Let be the estimated heart rate, be the manually computed heart rate, be the estimated respiratory rate, be the manually computed respiratory rate and N be the total number of the samples of the entire record within the 30 s duration. Let and be the accuracy of the estimated heart rate and the accuracy of the estimated respiratory rate, respectively. More precisely, they are defined as and The results were shown in Table 1. It can be seen from Table 1 that the digital filtering approach achieves the lowest accuracies on the estimations of both the heart rate and the respiratory rate. On the other hand, the empirical mode decomposition based method achieves the similar accuracy on the estimation of the heart rate compared to that of our proposed method. However, the empirical mode decomposition based method achieves a lower accuracy on the estimation of the respiratory rate compared to that based on our proposed method. In particular, the of the subjects with the subject numbers 055m and 220m achieved by our proposed method are 7.94% and 10.73% higher than those achieved by the empirical mode decomposition based method, respectively. These are the significant results. Overall, both the means and the variances of both the and the over all these 90 subjects achieved by our proposed method, the empirical mode decomposition based method and the digital filtering approach are computed and listed in Table 2. From here, it can be seen that our proposed method achieves the To quantitatively evaluate the performance of our proposed method, both our proposed method and the empirical mode decomposition based method are applied to the signals with the lengths equal to 30 s durations. Here, the signals with only 30 s durations are used for evaluating the numerical performances. This is because it can save the required computational powers. Let HR i be the estimated heart rate, HR i be the manually computed heart rate, RR i be the estimated respiratory rate, RR i be the manually computed respiratory rate and N be the total number of the samples of the entire record within the 30 s duration. Let ACC HR and ACC RR be the accuracy of the estimated heart rate and the accuracy of the estimated respiratory rate, respectively. More precisely, they are defined as and The results were shown in Table 1. It can be seen from Table 1 that the digital filtering approach achieves the lowest accuracies on the estimations of both the heart rate and the respiratory rate. On the other hand, the empirical mode decomposition based method achieves the similar accuracy on the estimation of the heart rate compared to that of our proposed method. However, the empirical mode decomposition based method achieves a lower accuracy on the estimation of the respiratory rate compared to that based on our proposed method. In particular, the ACC RR of the subjects with the subject numbers 055m and 220m achieved by our proposed method are 7.94% and 10.73% higher than those achieved by the empirical mode decomposition based method, respectively. These are the significant results. Overall, both the means and the variances of both the ACC HR and the ACC RR over all these 90 subjects achieved by our proposed method, the empirical mode decomposition based method and the digital filtering approach are computed and listed in Table 2. From here, it can be seen that our proposed method achieves the highest average accuracies compared to both the digital filtering approach and the empirical mode decomposition based method. In addition, our proposed method is more reliable compared to both the digital filtering approach and the empirical mode decomposition based method.  Figure 4a,b show the histograms of the absolute errors of the estimated heart rates obtained by the empirical mode decomposition based method and our proposed method, respectively. Figure 4c,d show the histograms of the absolute errors of the estimated respiratory rates obtained by the empirical mode decomposition based method and our proposed method, respectively. Let AE HR and AE RR be the absolute error of the estimated heart rate and the absolute error of the estimated respiratory rate, respectively. That is: and respectively. It can be seen From Figure 4a,b that the histogram of the absolute errors of the estimated heart rates obtained by the empirical mode decomposition based method is similar to that obtained by our proposed method. On the other hand, it can be seen from Figure 4c,d that there are some large values of AE RR based on the empirical mode decomposition based method (there are some bars in the histogram with small occurrences at the large absolute errors). Whereas, there is no large value of AE RR based on our proposed method. Therefore, our proposed method is more reliable than the empirical mode decomposition based method for computing the surrogate respiratory signal. The results on the estimation of both the heart rate and the respiratory rate of a PPG signal within the 270 s duration are shown in the Figure 5    Finally, the required computational power of the proposed method is evaluated. An Intel(R) Xeon(R) E3-1225 V2 CPU operating at 3.2 GHz with a 16 GB memory is employed for performing the computer numerical simulations. All the algorithms are executed using the Matlab Version 7.11.0.584 (R2010b) operating under the 64 bit Microsoft Windows 7 Version 6.1 with Service Pack 1 and Java 1.6.0_17-b04. It is found that the required computational time for processing a signal with the 30 s duration based on our proposed method is 0.12 s, which is acceptable in most real time applications.

Conclusions
This paper applies the complementary ensemble empirical mode decomposition, the independent component analysis and the non-negative matrix factorization to separate the PPG into two sets of signals as well as the principal component analysis to fuse these two sets of signals to generate two surrogate signals. These two surrogate signals are used to calculate the heart rate and the respiratory rate, respectively. The computer numerical simulation results show that our proposed method could achieve both the more reliable and the more accurate results compared to both the digital filtering approach and the conventional empirical mode decomposition based method.