Introducing the hybrid “K-means, RLS” learning for the RBF network in obstructive apnea disease detection using Dual-tree complex wavelet transform based features

Signal segmentation with the time duration of 60 s is the first stage of preprocessing. It has been shown in [25] that the 60 seconds duration is the most suitable period for apnea classification. After this, we have the noise cancelation procedure where the baseline wandering and power line interference of the signals are canceled using a Chebyshev type II band-pass filter with the low and high cut-off frequencies of 0.5 Hz and 48 Hz, respectively. Fig. 3, shows only 3 seconds of a model record of the Physionet apnea database [24].

Fig. 3

After we performed filtering, the weight calculation approach is applied for deleting the noisy segments. In [1] a simple method is proposed for the automatic cancelation of the noisy parts. In this method, a weight (W) is calculated for each segment based on the similarity of its Autocorrelation Function (ACF) with other segments ACF, by taking into account the cosine pair-wise similarity as the metric. The similarity values are then given as: (2) dst=(Xs−X¯s)(Xt−X¯t)′(Xs−X¯s)(Xs−X¯s)′(Xt−X¯t)(Xt−X¯t)′ d_{st} = {{(X_s - \bar X_s )(X_t - \bar X_t )^\prime} \over {\sqrt {(X_s - \bar X_s )(X_s - \bar X_s )^\prime} \sqrt {(X_t - \bar X_t )(X_t - \bar X_t )^\prime} }} where, d_st is the correlation distance, and X_s and X_t are the ACF of two different segments. X¯s \bar X_s and X¯t \bar X_t are the mean value of the X_s and X_t respectively. The output of (1) is a vector showing the similarity of each segment. W of each segment is calculated with the normalized summation of all the d_st values. For our method, all the segments with a calculated weight lower than 0.8 were pointed out as the noisy segments. As the results of this paper are presented for the segment-by-segment (minute-by-minute) case rather than the subject-by-subject case, the selection of these segments is very important.

In [20], the usage of DT-CWT in ECG feature extraction is proposed. The main deficiency of DWT based feature extraction in analyzing 1D ECG signals is the lack of shift invariance. It means that the amplitude of the wavelet coefficients varies substantially as the input signal is shifted a little. This happens because of the down sampling operation at each level. A better way of achieving shift invariance is to implement the undecimated form of the dyadic filter tree. However, this method has heavy computation demands and high redundancy in the output.

The DT-CWT tackles this problem with a redundancy factor for 1D signal, which is significantly lower than the undecimated DWT. In [20], the authors have explained the shift invariance property of DT-CWT in detail. The DT-CWT implements two trees of real filters DTCWT of a signal x(n) (ECG) is implemented using two critically sampled DWTs in parallel to the same data. The filters are (Tree A and Tree B) as shown in Fig. 1. The two trees correspond to the real and complex part of the complex wavelet transform. The designed so that the sub band signals of the upper DWT can be interpreted as the real part of a complex wavelet transform and sub band signals of the lower DWT can be interpreted as the imaginary part. When the transform is designed in this manner, the DT-DWT is approximately shift invariant, unlike the critically sampled DWT. The filters implemented in each stage are of length 10. The sets of filter coefficients (H) used in this transform are given in [20]. The selected transform coefficients are x_1a, x_01a, x_001a, x_000a, x_1b, x_01b, x_001b and x_000b.

Fig. 4

In Fig. 5 and 6 we depicted the sub band signals for three levels of the tree A and B respectively. It is important to mention that all of these signals are depicted for the same model record of the Fig. 3 of the Physionet database. It is important to mention that for the following figures from Fig. 5 to 8, the horizontal axis depicts the number of signal samples and the vertical axis shows the amplitude of the signal.

Fig. 5

Fig. 6

The absolute energy of the signal x_000a is depicted in Fig. 7, and the absolute energy of x_000b is depicted in Fig. 8.

Fig. 7

Fig. 8

After extracting the subbands of the DT-CWT from the selected ECG segments, we calculate some non-linear features based on the extracted transform coefficients. In [1] it has been shown that the ApEn, FE, IQR, RP and Poincare plot features make large differences among the two classes (Apnea and Normal). These features are collected in Table 2 and as they are explained in [1], we do not present their theoretical calculations here.

Using these 7 feature extraction methods and with the 8 DT-CWT coefficients that are explained in part III, we have 56 features for each ECG to be fed to the classifier. However, we used the feature reduction to reduce these features as much as possible.

After the extraction of the final features from the DT-CWT coefficients, it is time to select the most suitable ones to reduce the volume of the features and the amount of computation that is needed for processing them. In [1], the sequential forward feature selection (SFFS) method has been proposed for this task that is based on the detection results. Here we propose a data based feature reduction method. SRDA is one of the most efficient methods for feature reduction. We implemented this technique in our proposed method. To start the SRDA algorithm, suppose we have a set of data points x₁, ... x_m ∈ ℝ^N that belong to Nc different classes and m_k denotes the number of training samples of kth class (Σk=1Ncmk=m) (\Sigma _{k = 1}^{Nc} m_k = m) . The steps of SRDA can be summarized as [23]:

① Let (3) yk=[0,…,0︸∑i=1k−1mi,1,…,1︸mk,0,…,0︸∑i=k+1Ncmi]T k=1,…,Nc {\bf\it{y}}_k = \left[ {\underbrace {0, \ldots ,0}_{\sum\nolimits_{i = 1}^{k - 1} {m_i } },\;\underbrace {1, \ldots ,1}_{m_k },\;\underbrace {0, \ldots ,0}_{\sum\nolimits_{i = k + 1}^{Nc} {m_i } }} \right]^T \;\;\;\;\;\;\;\;\;k = 1, \ldots ,Nc and y₀ = [1,1, … ,1]^Trepresents the vector of ones. The Gram-Schmidt process [23] is utilized for orthogonalization of {y_k}. Since y₀ is in the subspace described by {y_k}, the Nc − 1 vectors are obtained as: (4) {y¯k}k=1Nc,(y¯iTy0=0 where y¯iTy¯j=0,i≠j) \left\{ {{\bf\it{\bar y}}_k } \right\}_{k = 1}^{Nc} ,\left( {{\bf\it{\bar y}}_i^T {\bf\it{y}}_0 = 0\;where\;{\bf\it{\bar y}}_i^T {\bf\it{\bar y}}_j = 0,i \ne j} \right)

② In this step, a new entry “l” is appended to each x_i which will also be shown as x_i. Then, Nc − 1 verctors {ak}k=1Nc−1∈ℝN+1 \left\{ {{\bf\it{a}}_k } \right\}_{k = 1}^{Nc - 1} \in {\mathbb R}^{N + 1} are constructed, where a_k is the solution of the regularized least squares problem given in (5): (5) ak=arg⁡ min⁡a(∑i=1m(aTxi−y¯ik)2+α‖a‖2) {\bf\it{a}}_k = \mathop {\arg \;\min }\limits_{\bf\it{a}} \left( {\sum\nolimits_{i = 1}^m {\left( {{\bf\it{a}}^T {\bf\it{x}}_i - \bar y_i^k } \right)^2 + \alpha \left\| {\bf\it{a}} \right\|} ^2 } \right) where y¯ik \bar y_i^k is the ith element of y¯k {\bf\it{\bar y}}_k and a ≥ 0 is a parameter to control the amount of reduction.

③ The Nc − 1 verctors {a_k} are the basis vectors of SRDA. Let A = [a₁, ... , a_Nc−1], which is a (N + 1) × (Nc − 1) transformation matrix. The x can be embedded into z in the (Nc − 1) dimension subspace by: (6) z=AT[x1] {\bf\it{z}} = A^T \left[ {\matrix{ {\bf\it{x}} \hfill \cr 1 \hfill \cr } } \right] The SRDA method reduces the 56 features that are mentioned in part B, based on the value of α which is between [0, 1].

The SVMs are the most prevalently used classifiers in the field of the disease detection. The RBF networks, on the other hand, are not used as much as SVMs. The hybrid RBF network is the solution for this, because, they can rival the SVMs. The hybrid RBF [21,22] consists of three layers and the middle and the output layers work with the K-means and the RLS algorithms, respectively and for this reason, the hybrid adjective is attributed to them.

In this part, we describe the RBF classifier with hybrid learning that is our proposed classifying tool in this paper. We call the proposed classifier as hybrid RBF because it has a hybrid learning procedure with two stages as follows:

Stage 1: implements the K-means clustering algorithm to train the hidden layer in an unsupervised scheme. Usually, the number of clusters and the computational units in the hidden layer is notably smaller than the size of the train sample.

Stage 2: implements the recursive least-squares (RLS) algorithm (or another adaptive algorithm) to determine the weight vector of the linear output layer.

The two-stage design procedure has some desirable features such as low computational complexity and fast convergence.

The RBF network consists of 3 layers as in Fig. 9. Here we describe them briefly [21]:

Fig. 9

1. Input layer, which contains the source nodes that connected the network to its inputs. The inputs of the network for classification are the features vectors.

2. The second layer, consisting of hidden units, implements a nonlinear transformation from the input space to hidden (feature) space. For most applications, the dimensionality of the only hidden layer of the network is high; this layer is trained in an unsupervised manner using stage 1 of the hybrid learning scheme. Each unit in the hidden layer is described mathematically by a radial basis function: (7) φj(x)=φ(‖x−xj‖) j=1,2,…,N \varphi _j ({\bf\it{x}}) = \varphi \left( {\left\| {{\bf\it{x}} - {\bf\it{x}}_j } \right\|} \right)\;\;\;\;\;j = 1,2, \ldots ,N The jth input data point x_j identifies the center of the radial basis function and the vector x is the signal (pattern) applied to the input layer. Therefore, the links connecting the source nodes to the hidden units are direct connections with no weights. There are multiple radial basis functions for using in the hidden layer but we implement the Gaussian function for the sake of comparison between SVM and RBF in [1].

3. The output layer is linear and provides the response of the network to the activation pattern implemented to the input layer; this layer is trained in a supervised fashion using stage 2 of the hybrid scheme. There is no limitation on the size of the output layer, except that typically, the size of the output layer is much smaller than that of the hidden layer.

Here we describe learning algorithms of RBF:

K-means is a method that utilizes distances for clustering with two steps [21, 22]:

Step 1: The total cluster variance is minimized with respect to the assigned set of cluster means {μ^}j=1K \left\{ {{\bf\it{\hat \mu }}} \right\}_{j = 1}^K , the following minimization must be performed: (8) min⁡{μ^}j=1K∑j=1K∑C(i)=j‖xi−μ^j‖2 for a given C \mathop {\it min }\limits_{\left\{ {{\bf\it{\hat \mu }}} \right\}_{j = 1}^K } \sum\nolimits_{j = 1}^K {\sum\nolimits_{C(i) = j} {\left\| {{\bf\it{x}}_i - {\bf\it{\hat \mu }}_j } \right\|^2 } \;\;\;\;\;for\;a\;given\;C}

Step 2: After computing the optimized cluster means {μ^}j=1K \left\{ {{\bf\it{\hat \mu }}} \right\}_{j = 1}^K , we optimize the encoder as follows: (9) C(i)=arg⁡ min⁡1≤j≤K‖xi−μ^j‖2 C(i) = \it arg \;\mathop {\it min }\limits_{1 \le j \le K} \left\| {{\bf\it{x}}_i - {\bf\it{\hat \mu }}_j } \right\|^2

Adaptive algorithms have been designed to converge to certain weights. These weights in the RBF network are adjusted in the learning phase. The RLS algorithm is one of the most powerful adaptive algorithms. In this section, we explain the role of RLS in the output layer of the RBF network [21]. Let the K × 1 vector: (10) Φ(xi)=[φ(xi,μ1)φ(xi,μ2)⋮φ(xi,μK)] {\bf\it{\Phi }}({\bf\it{x}}_i ) = \left[ {\matrix{ {\varphi ({\bf\it{x}}_{i,} {\bf\it{\mu }}_1 )} \cr {\varphi ({\bf\it{x}}_{i,} {\bf\it{\mu }}_2 )} \cr \vdots \cr {\varphi ({\bf\it{x}}_{i,} {\bf\it{\mu }}_K )} \cr } } \right] represent the outputs of the K units in the hidden layer. This vector is constructed to respond to the stimulus x_i, i = 1,2, … , N.Thus, insofar as the supervised training of the output layer is concerned, the training sample is defined by {Φ(i),d(i)}i=1N \left\{ {{\bf\it{\Phi }}(i),d(i)} \right\}_{i = 1}^N , where d_i is the desired response at the overall output of the RBF network for input x_i. This training is implemented by the RLS algorithm described as below [21]:

Given the training sample {Φ(i),d(i)}i=1N \left\{ {{\bf\it{\Phi }}(i),d(i)} \right\}_{i = 1}^N , do the following calculations for iterations n = 1,2, … , N: (11) P(n)=P(n−1)−P(n−1)Φ(n)ΦT(n)P(n−1)1+ΦT(n)P(n−1)Φ(n) {\bf\it{P}}(n) = {\bf\it{P}}(n - 1) - {{{\bf\it{P}}(n - 1){\bf\it{\Phi }}(n){\bf\it{\Phi }}^{\bf\it{T}} (n){\bf\it{P}}(n - 1)} \over {1 + {\bf\it{\Phi }}^{\bf\it{T}} (n){\bf\it{P}}(n - 1){\bf\it{\Phi }}(n)}} (12) g(n)=P(n)Φ(n) {\bf\it{g}}(n) = {\bf\it{P}}(n){\bf\it{\Phi }}(n) (13) α(n)=d(n)−w^T(n−1)Φ(n) \alpha (n) = d(n) - {\bf\it{\hat w}}^{\bf\it{T}} (n - 1){\bf\it{\Phi }}(n) (14) w^(n)=w^(n−1)+g(n)α(n) {\bf\it{\hat w}}(n) = {\bf\it{\hat w}}(n - 1) + {\bf\it{g}}(n)\alpha (n)

To initialize the algorithm, we have w^(0)=0 {\bf\it{\hat w}}(0) = {\bf\it{0}} and P(0) = λ⁻¹I where λ is a small positive constant.

In [222] a complete analysis was made to show the superiority of hybrid RBF to the SVM classifier both computationally and with respect to accuracy. Also, at least a 30% percent time saving is guaranteed using RBF in comparison with SVM [21]. This is important because we compared the results of our proposed method with that of the reference [1].

The conducted research is not related to either human or animal use.