Neural Comput & Applic DOI 10.1007/s00521-017-2865-3
NEW TRENDS IN DATA PRE-PROCESSING METHODS FOR SIGNAL AND IMAGE CLASSIFICATION
A novel system for automatic detection of K-complexes in sleep EEG ¨ zs¸ en1 • Gu¨lay Tezel2 • Serkan Ku¨c¸c¸u¨ktu¨rk3 • Cu¨neyt Yu¨celbas¸ 1 • S¸ ule Yu¨celbas¸ 2 • Seral O 3 S¸ ebnem Yosunkaya
Received: 27 October 2016 / Accepted: 3 February 2017 The Natural Computing Applications Forum 2017
Abstract Sleep staging process is applied to diagnose sleep-related disorders by sleep experts through analyzing sleep signals such as electroencephalogram (EEG), electrooculogram and electromyogram of subjects and determining the stages in 30-s-length time parts of sleep named as epochs. Subjects enter several stages during the sleep, and N-REM2 is one of them which has also the highest duration among the other stages. Approximately half of the sleep consists of N-REM2. One of the important parameters in determining N-REM2 stage is K-complex (Kc). In this study, some time and frequency analysis methods were used to determine the locations of Kcs, automatically. These are singular value decomposition (SVD), variational mode decomposition and discrete wavelet transform. The performance of them in detecting Kcs was compared.
Furthermore, systems with combinations of these methods were presented with logic AND operations. The EEG recordings of seven subjects were obtained from the Sleep Research Laboratory of Necmettin Erbakan University. A database with total 359 Kcs in 320 epochs was prepared from the recordings. According to the results, the highest average recognition rate was found as 92.29% for the SVD method. Thanks to this study, the sleep experts can find out whether there were Kcs in related epochs and also know their locations in these epochs, automatically. Also, it will help automatic sleep stage classification systems. Keywords DWT Sleep EEG K-complex SVD VMD
1 Introduction & Cu¨neyt Yu¨celbas¸
[email protected] S¸ ule Yu¨celbas¸
[email protected] ¨ zs¸ en Seral O
[email protected] Gu¨lay Tezel
[email protected] Serkan Ku¨c¸c¸u¨ktu¨rk
[email protected] S¸ ebnem Yosunkaya
[email protected] 1
Electrical and Electronics Engineering Department, Selcuk University, 42072 Konya, Turkey
2
Computer Engineering Department, Selcuk University, Konya, Turkey
3
Sleep Laboratory, Faculty of Medicine, Necmettin Erbakan University, Konya, Turkey
One of the important factors affecting people’s lives is sleep quality. Nowadays, sleep stage scoring process has a significant role for detecting the problems associated with sleep. It is an ineluctable method to detect the sleep-related problems. To score sleep, it was benefited from the signals such as electroencephalogram (EEG), electrooculogram (EOG), electromyogram (EMG), and electrocardiogram (ECG) during sleep, and then, these signals are evaluated by sleep experts. There are some rules for staging sleep. The specific rules are defined by [1]. According to the guideline, a night sleep is divided into five stages which are named as Awake (W), Non-REM1, Non-REM2, NonREM3 and REM [2]. K-complex is one of the important parameters to label an epoch as N-REM stage-2. Due to this significance, the determination of Kc in an epoch is extremely important for sleep experts.
123
Neural Comput & Applic
K-complex is defined as temporary waveform which is seen by negative sharp wave and then followed by positive wave [3] and also other definition is exactly the opposite of this study [3] which consists of a positive wave and followed by a negative component [4, 5]. The waveform’s duration is between 0.5 and 1.5 s and mostly bigger than 75 lV amplitude [5, 6]. Also in other study [7], its duration is between 0.5 and 1 s and the minimum amplitude value of the Kcs is 100 lV (peak to peak). Visual identification of Kcs in sleep EEG is difficult, tiring and time-consuming job. Sleep experts usually hesitate about the time duration and amplitude level of Kcs. The features of Kc are generally confused. In the literature, many studies have been performed related to automatic detection of Kcs. Jansen et al. [8] compared two knowledge-based systems to determine the Kcs in sleep EEG. In another study [4], they introduced a feature-based Kc detection system by neural networks that obtained a sensitivity rate of 90%. In study [9], researchers concentrated on the detection of Kcs in sleep EEG using methods based on decomposition of the signal. Researchers in [5] utilized the time–frequency-based detection schemes to identify the Kcs. In [10], a continuous density hidden Markov model was presented for the detection of Kcs from EEG signals. In study [11], an algorithm was presented to detect the Kcs from EEG signals. An automatic method by using feature extraction and fuzzy threshold was introduced for determination of Kcs in study [12]. In other study, the researchers were proposed a two-step Kc detection method which the Kc candidates are determined in first step and they are reduced by a machine learning algorithm for the second [13]. In study [14], a detection algorithm by using wavelet method was presented to determine Kc density in N2-sleep [14]. In [15], various features were extracted from Kc waves in N-REM2 sleep and the features were classified by a machine learning algorithm. In other study, researchers were surveyed the automatic detection of stage 2 sleep in EEG signals [3]. Zamir et al. [16] were introduced 3 convex optimization models that obtain important features of EEG data for detecting Kcs. In other research study [17], usually used automatic determination systems were proposed for sleep spindles and Kcs. Pereira et al. [18] experimented the capabilities of several machine learning techniques. Mostly, researchers developed and presented a procedure to identify Kcs in an epoch by classification. In this study, however, singular value decomposition (SVD), variational mode decomposition (VMD) and discrete wavelet transform (DWT) were utilized to determine the locations of Kcs. Comparison of these methods were done. Also, systems combining these methods by logic AND gates were proposed and used differently from the studies in the literature. A database was prepared by using the EEG waves of seven
123
subjects. The database consists of 320 EEG epochs which are recorded from C4A1 and C3A2 channels, and these epochs include 359 Kcs in total. The locations of 359 Kcs in 320 different epochs were identified by two sleep experts who have professional certificates. According to the results, the highest average recognition rate (ARR) was found as 92.29% by SVD method. Thanks to this study, the sleep experts can find out whether there were Kcs in related epochs and also know their locations in these epochs, automatically. Besides, detecting locations of Kcs in an epoch would help an automatic sleep staging system in a great deal.
2 Methods In this study, the used EEG signals—with a sampling frequency of 128 Hz—of 7 subjects recorded by a PSG equipment (VIASY) in sleep laboratory of Necmettin Erbakan University. For each subject, the data were recorded for about 8 h. And then, the epochs (signal parts in 30 s long) of the data were labeled as W, Non-REM1–2– 3 and REM by two sleep experts. The experimental protocol conformed to the ethical standard. 2.1 Sleep staging and K-complex The sleep staging has being performed according to the R&K rules [1]. SSs and K-complexes are among the decisive the important parameters in determination of N-REM2 [2]. Kc is defined as temporary waveform which is characterized by negative sharp wave and then followed by positive wave [3]. The waveform’s duration is between 0.5 and 1.5 s and mostly bigger than 75 lV amplitude [5, 6]. Also in other study [7], its duration is between 0.5 and 1 s and the minimum amplitude value of the Kcs is 100 lV (peak to peak). The example signals with/without Kcs are presented in Fig. 1. 2.2 Dataset In experiments, although the other signals which are EOG, EMG, ECG, chest, body and leg signals were recorded during the night, we only used C4A1 and C3A2 channels of EEG because Kcs are more pronounced in these channels. Filtering (band-pass filter in 0.5–35 Hz band) and eliminating artifacts were applied to the C3A2 and C4A1 channels of EEG signals. After that, the epochs in 30 s length were prepared according to sampling rate. For this study, we only used epochs belonging to N-REM2 to detect the locations of Kcs. The positions of 359 Kcs in 320 different epochs were labeled and determined by the experts. The subjects used in the study are given in Table 1.
Neural Comput & Applic
Fig. 1 Example EEG signals with/without Kc. a An example epoch with K-complex (Kc) and b an example epoch without K-complex (Kc) Table 1 The subjects in used dataset Subject-1
Subject-2
Subject-3
Subject-4
Subject-5
Subject-6
Subject-7
The total number of epochs
1006
1105
751
963
969
898
916
Time in Bed (min)
503
552.5
375.5
481.5
484.5
449
458
Total sleep time (TST) (min)
493.5
472
335.5
383
454
415
386.5
Percent of sleep efficiency (%)
98
85
89
80
94
92
84
Percent of TST in stage 2 (%)
75.7
70.6
76.8
60.2
72.2
66.5
58.7
Gender (M: male–F: female)
M
M
M
M
M
M
F
Age
46
60
46
45
44
26
60
Height (cm)
162
174
178
174
176
177
170
Weight (kg)
70
80
83
102
86
74
80
The conducted procedure used in this study is demonstrated in Fig. 2. In this study, the ECG artifacts were cleaned from the EEG signal before performing process on EEG signals. To remove the ECG artifacts, R peaks in ECG signals were detected by Pan–Tompkins algorithm [19]. 2.3 Singular value decomposition (SVD) The SVD method is an useful tool which can be used in dimensionality reduction, statistical data analysis, signal processing problems and applications [20]. Singular value decomposition could decompose of an m 9 n matrix A = [a1, a2, a3,…, an] into three matrices [21]: A ¼ UKV T
ð1Þ
where K is an m 9 n diagonal matrix, whose diagonal components are the singular values of matrix Am9n that are arranged in descending order such that r1 C r2
… C rn C 0. The left and right singular vectors of matrix A are the columns of the orthogonal Um9m and Vn9n matrices, respectively [21]. The values indicate the inherent characteristic of the matrix and externalize the energy info included in every subspace. In other words, the greater values have more energy and more info about the structure of samples [21, 22]. In this study, the maximum singular values in each 1-slong EEG signal were taken as features to determine the locations of Kcs as seen in the example application in Fig. 3. 2.4 Variational mode decomposition (VMD) Empirical mode decomposition (EMD) was first acquainted by [23]. This is a recursive technique which can be utilized for nonlinear and non-stationary signals such as EEG and EOG [23]. In this technique, a signal adaptively separated into substatements named as intrinsic mode functions
123
Neural Comput & Applic
than EMD method in tone perception, separation and artifact elimination. On account of this, VMD is more appropriate for using in non-stationary and nonlinear signals [25]. As a brief description of this method, VMD can divide any signal in time domain—x(t)—into K separate number of subsignals or modes—uk—where every mode of the signal is assumed to have intensive frequency around a central frequency—wk [24, 25]. A detailed description of VMD method has given in [24]. ( 2 ) XK j jw t k ot dðtÞ þ min ðtÞ e u k uu k¼1 pt fuk g;fwk g 2 ð2Þ Subject to K X
uk ðtÞ ¼ x
ð3Þ
k¼1
Fig. 2 Conducted procedure
(IMFs). A new method was introduced a non-iterative alternative to the EMD so-called variational mode decomposition (VMD) by [24]. VMD can dynamically divide any signal into subcomponents (modes—uk) that are community of K band-limited intrinsic mode functions (BLIMFs). These experiment results in study [25] have indicated that the VMD performs better
where t is the time domain, d is the Dirac distribution and * signifies convolution. And also, L and k denote Lagrangian and Lagrangian multiplier, respectively. The Lagrangian is given by [24, 25]. 2 K X j jwk t Lðfuk g; fwk g; kÞ ¼ a ot dðtÞ þ pt uk ðtÞ e 2 k¼1 2 * + K K X X þ xðtÞ uk ðtÞ þ kðtÞ; xðtÞ uk ðtÞ k¼1 k¼1 2
ð4Þ in which a defines the balancing parameter. After that, Eq. (3) is solved by the ADMM. The modes in time
Fig. 3 Example SVD implementation and its output. a An example epoch with K-complex and b SVD application of the epoch
123
Neural Comput & Applic
domain—uk(t)—are obtained by the following equations [24, 25]: P ^ x^ðwÞ i6¼k u^i ðwÞ þ kðwÞ 2 nþ1 u^k ðwÞ ¼ ð5Þ 2 1 þ 2aðw wk Þ where x^ðwÞ is the Fourier transform of the signal x(t). The updated equation of wk is shown as [24, 25]: R1 wju^k ðwÞj2 dw nþ1 ^k ¼ R0 1 w ð6Þ ^k ðwÞj2 dw 0 ju when the modes (^ uk ) are analyzed according to frequency components. The primary u^k contains the low frequency info and the last u^k includes higher-frequency information. In this study, the first u^k in modes was operated on because of the locations of Kcs were more prominent and more specific in the signal as shown in Fig. 4.
2.5 Discrete wavelet transform (DWT) The wavelet transform (WT) can be described as a widening version of the classic Fourier transform, but it runs on a multi-scale principle. This multi-scale lets the separation of a signal into a certain number scales. The formula of discrete wavelet transform (DWT) is as follows [26]: 1 DWT ðj; kÞ ¼ pffiffiffiffiffiffiffi j2 j j
Z1
t 2 jk xðtÞ w dt 2j
ð7Þ
1
where w is function and the 2j and 2jk are scrolling factors, respectively. Mallat [27] improved an efficient method for execution this method by using low- and high-pass filter on the signal. In this paper, five-level DWT was used and at the end of the decomposition, and the approximation coefficients (0–2 Hz) in fifth level were utilized. Coefficients in the fifth level are the nearest to the frequency range of Kc. Also, two different wavelet functions (db2 and db4) which are the most similar to Kc were used for this application of DWT method as shown in Figs. 5 and 6. 2.6 Thresholding and specification of Kcs One of the most substantial sections of this paper is thresholding. This step was administrated to identify the locations of Kcs. The thresholding method used in this study was dynamically designed for each epoch. According to this method, after the implementation results obtained for the related epoch, mean and standard deviation values of the obtained waveform in related epoch were computed
separately. And then, the threshold values were calculated by multiplying with specific coefficients places of Kcs that were found by detecting higher-level zones in these obtained new signals. The new signals were obtained by using the SVD, VMD and DWT methods. A sample application for thresholding is shown in Fig. 7. An important parameter that influences the success of thresholding process is the k parameter which plays on important role in detection of the threshold level. In this study, two different calculation methods were experimented to determine the levels as shown in Eqs. 8 and 9. Threshold level 1 ¼ std ðnew signalÞ k
ð8Þ
Threshold level 2 ¼ mean ðnew signalÞ k
ð9Þ
where k is a fixed value which is given as k = 0.1:0.1:10. The threshold level was calculated thought multiplying k by either standard deviation or mean of the new signal for each epoch. Lesser values result many zones surpassing the level, while greater values cause fewer (sometimes none) zones to be identified. As illustrated in Fig. 8 when this value is chosen to be high, only very high signal components will exceed the threshold. The lesser values of this parameter, however, cause so much confusion in that almost all the particles in the signal will exceed this level (Fig. 8). In these implementations, the most appropriate values for k parameter were identified by trial and error. The performance was analyzed for various k values. Also, each calculation method for the threshold level as given in Eqs. 8 and 9 was experimented. In the first place, SVD, VMD and DWT methods were used to detect Kcs separately and the performance of them was compared. 2.7 Application of SVD, VMD and DWT methods The application of SVD, VMD and DWT methods is illustrated in Fig. 9 schematically. As shown in Fig. 9, each of the 30-s-long epochs was initially divided into 30 parts in 1 s length for SVD, VMD and DWT methods. In the next step, these methods were applied to these 1-s-long parts of each epoch. The first singular value (SV) of each 1-s-long windows was obtained in SVD method, and a new signal composed of 30 SVs for an epoch was obtained. In the VMD and DWT methods, however, maximum of absolute values in each 1-s-long window were taken. For example, in VMD method first IMF component of 30-s-long epoch was divided into 1-s-long windows and maximum of absolute values in each windows were taken. Thus, a new signal whose length is 30 data points was comprised by combining these maximum values of each window. In DWT
123
Neural Comput & Applic
Fig. 4 Sample application for VMD method. a An example epoch with K-complex and b VMD application of the epoch Fig. 5 db2 and db4 wavelet functions resembling K-complex as shape
method, this procedure is conducted on fifth-level approximation coefficient instead of first IMF of VMD. Next, the new data were obtained in each of SVD, VMD and DWT methods. And then, they were entitled as ‘SVD new data,’ ‘VMD new data’ and ‘DWT new data.’ After that, the thresholding was applied to detect Kcs with different k parameter values for all methods. Finally, a confusion matrix was formed to calculate sensitivity, specificity and average recognition rates for each method. Evaluation of used methods in this system was carried out over 9600 windows (320 epochs * 30 windows = 9600). Kcs in the epoch were labeled as 1 and 0 as shown in Fig. 10. The 1 and 0 indicate the Yes-Kc and No-Kc in the related epoch, respectively. If there are two adjoining ‘1’ in an epoch, they were taken as only one Kc. Checking of this situation was carried out programmatically.
123
At the end of all these practices, the best method/ methods was/were determined by evaluating the obtained results of this system. After application of SVD, VMD and DWT methods separately, a methodology by combining these methods was proposed with the use of logic AND gate concept. In this system, SVD, VMD and DWT methods were used to label KCs in each 1-s-long windows as the first stage. In the second step, the results obtained from these methods were given as inputs to the logic AND gate system. Thus, it was targeted to increase the sensitivity ratio. The application schema of this system is shown in Fig. 11. As shown in Fig. 11, four different methods were applied in this proposed system. In all methods, the data (‘SVD new data,’ ‘VMD new data’ and ‘DWT new data’) obtained in the proposed first system were used as input. The new data include 0 and 1 values which 0 and 1 indicate
Neural Comput & Applic
Fig. 6 Sample application for DWT method. a An example epoch with K-complex, b DWT application of the epoch for db2 mother wavelet and c DWT application of the epoch for db4 mother wavelet
Fig. 7 Example implementation and thresholding. a An example epoch with K-complex, b locations of the K-complexes that labeled by Expert1 and Expert-2, c application of the example epoch and threshold level and d application result
123
Neural Comput & Applic
Fig. 8 Example for determination of threshold level. a An example epoch with K-complex, b an example application for high threshold level and c an example application for low threshold level
Non-Kc and Kc, respectively. In method 1 of the proposed second system, SVD&VMD new data were given to AND gate as input. SVD&DWT and VMD&DWT new data were used for method 2 and method 3, respectively. In method 4, SVD&VMD&DWT new data were used for input. Finally, a confusion matrix was formed to calculate sensitivity, specificity and average recognition rates for each method. Again, evaluation of used methods in this system was carried out over 9600 windows (320 epochs * 30 windows = 9600). An example application of the proposed second system is shown in Fig. 12. Figure 12 shows that the sensitivity of the system would be increased by using AND implementation. As shown in Fig. 15a, there is only one Kc in the example epoch and in (b), whose location was labeled by two sleep experts as agreed. VMD and DWT application results of the epoch were seen in (c) and (d). According to (c), one location was correctly found and the other was incorrectly labeled as Kc. In (d), the one was correctly and the others were incorrectly labeled as Kc. When the (c) and (d) were given, the proposed AND system as input, and the (e) was obtained. According to the application result of the VMD&DWT logic AND gate system, only one location was determined as Kc, correctly. As a result, the proposed system increased the performances in detecting Kcs.
123
2.8 Test and performance assessment In all used systems, the sensitivity (SN), specificity (SP) and average recognition rate (ARR) measurements were utilized to see performances of the systems in detecting Kcs. A confusion matrix was used for the determination of these parameters (TP, TN, FP, FN) as given in Table 2. The SN and SP performance criteria are determined by using the matrix as follows: SN ¼ TP=ðTP þ FNÞ
ð10Þ
SP ¼ TN=ðTN þ FPÞ
ð11Þ
ARR ¼ ðSN þ SPÞ=2
ð12Þ
Calculation of TP, FP and FN values was simple. Periods of signals passing threshold level were computed with the aid of locations found in thresholding application. The systems searched the Kcs and labeled them if they are really Kc or not, and then, TP and FP values were obtained. FN value was calculated by taking real Kcs and examining whether they were labeled as not-Kcs. Calculation of TN value, however, had not so easy methodology. We took the signal places which didn’t labeled as Kcs by the systems as shown in Fig. 13. As shown in Fig. 13b, the locations that except TP, FP and FN were taken as TN. The TN represents parts in which there are not-Kcs according to both the system and sleep experts. Then, these parts were divided into 1-s-long parts to determine the number of TN.
Neural Comput & Applic
Fig. 9 Experimental procedure applied for SVD, VMD and DWT methods
123
Neural Comput & Applic Fig. 10 Labeling of an epoch as ‘1’ and ‘0’
Fig. 11 Experimental procedure conducted for system-2
123
Neural Comput & Applic
Fig. 12 Example application of the logic AND gate system. a An example epoch with K-complex, b locations of the K-complex that labeled by expert 1 and expert 2, c VMD application result of the
epoch, d DWT application result of the epoch for db2 mother wavelet and e application result of the VMD&DWT logic AND gate method
Table 2 Confusion matrix
obtained results were compared. The first implemented method was SVD. The singular values of each window were obtained to determine the Kcs. Kcs were detected by searching duration of signal components passing the threshold level. Optimum values for k parameter were searched in this method. The application results of this method are given in Table 3. Sensitivity value explains at which degree a system can detect Kcs. However, specify value gives the incorrectly detected Kcs by the system when in fact there were no Kcs. As given from Table 3, Sens. and Spec. rates were obtained as 100% but for different threshold values. The highest ARR was obtained with k value of 1.4 when the mean method was used in threshold calculation. The changes of Sens., Spec. and ARR rates according to the k parameter are shown in Fig. 14. The second implemented method was VMD which is applied as explained in Sect. 2. Again, Kcs were determined by controlling duration of signal windows passing the threshold level. Optimum values for k parameter were
Predicted class Yes (Kc)
No (not-Kc)
Yes (Kc)
True positive (TP)
False negative (FN)
No (not-Kc)
False positive (FP)
True negative (TN)
Actual class
3 Results In this paper, detection of the locations of Kcs in an epoch was aimed. Furthermore, the performances of three different methods in proposed two systems were compared as seen in result tables. 3.1 Comparison of SVD, VMD and DWT methods As a first study step, SVD, VMD and DWT methods were applied as explained in the previous section and the
123
Neural Comput & Applic
Fig. 13 Thresholding and determination of confusion matrix elements. a An example epoch with K-complex and b thresholding and determination of TP, TN, FP AND FN values
Table 3 The results of the SVD method in system 1
Tvcm Mean
Standard deviation (std)
Rates
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%) 57.39
HSenR
0.6
359
7281
1262
0
100
14.77
HSpeR
5.1
6
0
9229
353
1.67
100
50.84
HARR
1.4
337
798
7787
22
93.87
90.70
92.29
HSenR
0.9
359
7601
942
0
100
11.03
55.51
HSpeR
6.9
7
0
9227
352
1.95
100
50.97
HARR
3.6
313
845
7784
46
87.19
90.21
88.70
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate
searched in this method, too. The application results of this method are given in Table 4. As given in Table 4, Sens. and Spec. rates were also obtained as 100% but similarly for different threshold values. The highest ARR was obtained with k value of 1.5 when the mean calculation method was used in determining threshold level. The changes of Sens., Spec. and ARR rates according to the k parameter are shown in Fig. 15. The last implemented method was DWT. Similar to SVD and VMD, Kcs were detected by controlling duration
123
of signal parts passing the threshold level. Optimum values for k parameter were searched in this method. The application results of this method are given in Table 5. In this application, two types of wavelet function—dubhecies 2 (db2) and dubhecies 4 (db4)—were also experimented differently from the previous two methods. Again, 100% Sens. and Spec. values were reached but for different k values as shown in Table 5. It can be found from the table that generally db4 gave better results, but for the best ARR value, db2 obtained higher value than db4. Also, unlike SVD and VMD methods, std calculation
Neural Comput & Applic
Fig. 14 Change of Sens., Spec. and ARR rates according to the k parameter for ‘mean’ (a, c, e) and ‘std’ (b, d, f) in application of the SVD method
Table 4 The results of the VMD method in system 1
Tvcm Mean
Standard deviation (std)
Rates
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%) 62.24
HSenR
0.6
359
6452
2091
0
100
24.48
HSpeR
5.5
4
0
9233
355
1.11
100
50.56
HARR
1.5
323
835
7776
36
89.97
90.30
90.14
359
7528
1015
0
100
11.88
55.94
0
0
9241
359
0
100
50
323
965
7645
36
89.97
88.79
89.38
HSenR
0.7
HSpeR
7
HARR
2.9
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate
method for the threshold value gave better results. The highest ARR rate was obtained with k value of 2.5 and wavelet of db2 when the std calculation method was used in determining threshold level. The changes of Sens., Spec. and ARR rates according to the k parameter are shown in Fig. 16. 3.2 Results of combinations with logic AND gate Signals to be thresholded which was generated by SVD, VMD and DWT methods were given to AND lojik gate
as illustrated in Fig. 11. The purpose of this was to increase the probability of correctly identifying the Kcs. In briefly, it was intended to increase the sensitivity rates. The first implemented method was SVD&VMD. In this method, a two-input AND gate was used (see Fig. 11). The application results of this method are given in Table 6. Whereas 100% was obtained for Sens. and Spec. values, they were obtained for different k values. The ARR rates were recorded as 91.39 and 89.26% for ‘mean’ and ‘std’ threshold calculation methods, respectively.
123
Neural Comput & Applic
Fig. 15 Change of Sens., Spec. and ARR rates with respect to the k parameter for ‘mean’ (a, c, e) and ‘std’ (b, d, f) in application of the VMD method
Table 5 The results of the DWT method in system 1 Tvcm
Rates
Mean
HSenR HSpeR HARR
Standard deviation (std)
WFT
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%)
db2
0.5
359
6768
1775
0
100
20.78
60.39
db4
0.6
359
5902
2641
0
100
30.91
65.46
db2
6.9
2
0
9237
357
0.56
100
50.28
db4
5.8
12
0
9217
347
3.34
100
51.67
db2
1.4
319
1239
7380
40
88.86
85.62
87.24
db4
1.4
318
1210
7407
41
HSenR
db2
0.5
359
7579
964
0
db4
0.7
359
6690
1853
0
100
21.69
60.85
HSpeR
db2
5.8
7
0
9227
352
1.95
100
50.97
db4
5.8
13
0
9215
346
3.62
100
51.81
db2 db4
2.5 2.7
321 304
1260 954
7355 7691
38 55
89.42 84.68
85.37 88.96
87.39 86.82
HARR
88.58
85.96
87.27
100
11.28
55.64
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate; WFT, wavelet function type
3.3 Results of SVD&DWT method The second implemented method was SVD&DWT (see Fig. 11). The application results of this method are given in Table 7.
123
The ARR rates were recorded as 89.48 and 87.28% for ‘mean’ and ‘std’ threshold calculation methods, respectively.
Neural Comput & Applic
Fig. 16 Change of Sens., Spec. and ARR rates with respect to the k parameter for ‘mean’ (a, c, e) and ‘std’ (b, d, f) in application of the DWT method Table 6 The results of the SVD&VMD method in system 2
Tvcm Mean
Standard deviation (std)
Rates
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%) 63.77
HSenR
0.6
359
6190
2353
0
100
27.54
HSpeR
5.1
4
0
9233
355
1.11
100
50.56
HARR
1.2
344
1118
7454
15
95.82
86.96
91.39
HSenR
0.7
359
7478
1065
0
100
12.47
56.23
HSpeR
6.1
12
0
9218
347
3.34
100
51.67
HARR
2.8
321
939
7674
38
89.42
89.10
89.26
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate
3.3.1 Results of VMD&DWT method The next implemented method was VMD&DWT. The application results of this method are given in Table 8. The ARR rates were recorded as 89.47 and 88.21% for ‘mean’ and ‘std’ threshold calculation methods, respectively. 3.3.2 Results of SVD&VMD&DWT method The last implemented method was SVD&VMD&DWT. The application results of this method are given in Table 9. The ARR rates were recorded as 90.06 and 88.26% for ‘mean’ and ‘std’ threshold calculation methods, respectively.
When examining the result tables, the Sens. and Spec. rates were obtained as 100% in both threshold value computation methods (‘mean’ and ‘std’). When assessed in terms of average recognition rates, the most obvious difference between the calculation methods was emerged. For the ARR, the ‘mean’ computation method was more advantageous and, furthermore, had achieved better results. According to these result tables, the best result was obtained in SVD&VMD method in which ARR was obtained as 91.39%. The aim of these combinations was to increase ARR rates by eliminating erroneous decisions but even in the best combination (in which ARR was 91.39%) the result of SVD (which was 92.29%) couldn’t exceeded. In the first place, we had aimed to increase
123
Neural Comput & Applic Table 7 The results of the SVD&DWT method in system 2 Tvcm
Rates
Mean
HSenR
Standard deviation (std)
WFT
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%)
db2
0.5
359
6724
1819
0
100
21.29
60.65
db4
0.6
359
5780
2763
0
100
32.34
66.17
HSpeR
db2
5.1
6
0
9229
353
1.67
100
58.84
db4
5.1
6
0
9229
353
1.67
100
50.84
HARR
db2
1.2
334
1210
7382
25
93.04
85.92
89.48
db4
1.2
332
1204
7388
27
92.48
85.99
89.23
HSenR
db2
0.5
359
7579
964
0
100
11.28
55.64
db4
0.7
359
6686
1857
0
100
21.74
60.87
HSpeR
db2
5.7
9
0
9223
350
2.51
100
51.25
db4
5.8
12
0
9217
347
3.34
100
51.67
HARR
db2 db4
2.3 2.7
330 304
1493 903
7105 7742
29 55
91.92 84.68
82.64 89.55
87.28 87.12
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate; WFT, wavelet function type
Table 8 The results of the VMD&DWT method in system 2 Tvcm
Rates
Mean
HSenR
Standard deviation (std)
WFT
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%)
db2
0.5
359
6328
2215
0
100
25.93
62.96
db4
0.6
359
5258
3285
0
100
38.45
69.23
HSpeR
db2
5.5
4
0
9233
355
1.11
100
50.56
db4
5.5
4
0
9233
355
1.11
100
50.56
HARR
db2
1.2
333
1187
7406
26
92.76
86.19
89.47
db4
1.2
330
1180
7416
29
91.92
86.27
89.10
db2
0.5
359
7486
1057
0
100
12.37
56.19
db4
0.6
359
7015
1528
0
100
17.89
58.94
db2
5.6
11
0
9219
348
3.06
100
51.53
db4
5.7
12
0
9217
347
3.34
100
51.67
db2 db4
2.3 2.3
320 310
1096 1066
7520 7568
39 49
89.14 86.35
87.28 87.65
88.21 87.00
HSenR HSpeR HARR
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate; WFT, wavelet function type
sensitivity by combining different methods with AND operation as shown in Fig. 12. We saw that sensitivity rates were increased by AND operation, but on the other hand, specificity ratios were also decreased because of FP rates were increased by the AND operation. Different strategies to combine the applied methods could be experimented to increase both sensitivity and specificity.
4 Discussion In this study, two systems were presented to automatically determine the locations of Kcs in sleep EEG that taken from seven subjects recorded at sleep laboratory of
123
Necmettin Erbakan University. In the first place, SVD, VMD and DWT methods were used, and then, the numerical values that are obtained from these methods were presented as input to the logic AND gate system. Example outcomes of these methods are shown in Fig. 17. In this paper, the comparison table of exhaustive test results of the all used methods is given in Table 10. The results were given for only the highest average recognition rate in comparison table. When the result tables of systems 1 and 2 were examined, the best sensitivity and specificity rates were obtained as 100% for all methods but at the expense of lower rates of the other. For example, 100% Sens. in SVD application was obtained for k = 0.6 but Spec. value for that k value
Neural Comput & Applic Table 9 The results of the SVD&VMD&DWT method in system 2 Tvcm
Rates
Mean
HSenR
Standard deviation (std)
WFT
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%)
db2
0.5
359
6319
2224
0
100
26.03
63.02
db4
0.6
359
5210
3333
0
100
39.01
69.51
HSpeR
db2
5.1
4
0
9233
355
1.11
100
50.56
db4
5.1
4
0
9233
355
1.11
100
50.56
HARR
db2
1.2
329
991
7610
30
91.64
88.48
90.06
db4
1.2
327
992
7610
32
91.09
88.47
89.78
db2
0.5
359
7486
1057
0
100
12.37
56.19
db4
0.6
359
7015
1528
0
100
17.89
58.94
HSpeR
db2
5.6
11
0
9219
348
3.06
100
51.53
db4
5.7
12
0
9217
347
3.34
100
51.67
HARR
db2 db4
2.3 2.3
320 310
1087 1055
7529 7579
39 49
89.14 86.35
87.38 87.78
88.26 87.07
HSenR
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate; WFT, wavelet function type
Fig. 17 Graphical comparison of the example outcomes for SVD, VMD and DWT methods. a An example epoch with K-complex, b SVD application of the epoch, c VMD application of the epoch and d DWT application of the epoch
123
Neural Comput & Applic Table 10 Comparison table of exhaustive results for methods in system 1 and system 2 Tvcm and rate type Mean and HARR
Standard deviation (std) and HARR
Methods
k
TP
FP
TN
FN
Sens. (%)
Spec. (%)
ARR (%)
SVD
1.4
337
798
7787
22
93.87
90.70
92.29
VMD
1.5
323
835
7776
36
89.97
90.30
90.14
DWT (db2)
1.4
319
1239
7380
40
88.86
85.62
87.24
DWT (db4)
1.4
318
1210
7407
41
88.58
85.96
87.27
SVD&VMD
1.2
344
1118
7454
15
95.82
86.96
91.39
SVD&DWT (db2)
1.2
334
1210
7382
25
93.04
85.92
89.48
SVD&DWT (db4)
1.2
332
1204
7388
27
92.48
85.99
89.23
VMD&DWT (db2)
1.2
333
1187
7406
26
92.76
86.19
89.47
VMD&DWT (db4)
1.2
330
1180
7416
29
91.92
86.27
89.10
SVD&VMD&DWT (db2)
1.2
329
991
7610
30
91.64
88.48
90.06
SVD&VMD&DWT (db4)
1.2
327
992
7610
32
91.09
88.47
89.78
SVD
3.6
313
845
7784
46
87.19
90.21
88.70
VMD
2.9
323
965
7645
36
89.97
88.79
89.38
DWT (db2)
2.5
321
1260
7355
38
89.42
85.37
87.39
DWT (db4)
2.7
304
954
7691
55
84.68
88.96
86.82
SVD&VMD
2.8
321
939
7674
38
89.42
89.10
89.26
SVD&DWT (db2)
2.3
330
1493
7105
29
91.92
82.64
87.28
SVD&DWT (db4)
2.7
304
903
7742
55
84.68
89.55
87.12
VMD&DWT (db2)
2.3
320
1096
7520
39
89.14
87.28
88.21
VMD&DWT (db4)
2.3
310
1066
7568
49
86.35
87.65
87.00
SVD&VMD&DWT (db2)
2.3
320
1087
7529
39
89.14
87.38
88.26
SVD&VMD&DWT (db4)
2.3
310
1055
7579
49
86.35
87.78
87.07
Tvcm, the threshold value computation methods; mean, thres = k*mean(data); std, thres = k*std(data); HSenR, the results for the highest sensitivity rate; HSpeR, the results for the highest specificity rate; HARR, the results for the highest average recognition rate
Fig. 18 Samples of the true K-complex (TKc) and false K-complex (FKc). a An example epoch with K-complex and b example SVD application of the epoch
123
Neural Comput & Applic Table 11 Comparison of experimented methods with the literature References
Methods
The number of subject
The number of K-complex
Sens. (%)
Spec. (%)
ARR (%)
This study
SVD (mean)
7
359
93.87
90.70
92.29
95.82
86.96
91.39
[4]
A feature-based detection approach using neural networks
2
Training set: 200 Non-Kc 200 Kc
90
–
–
Scorer 1: 221
80
–
–
Scorer 2: 280
89
Scorer 3: 240
87
SVD&VMD (mean)
Test set: 49 Non-Kc 51 Kc [11]
Wavelet and Teager energy operator
1
[12]
Features extraction and the use of fuzzy thresholds
5
– (unclear)
61.72
–
–
[13]
One-class classification by two-step methodology
1
279
83
–
84
[14]
Wavelet transform
5
160
74
–
–
[15]
Generalized radial basis function extreme learning machine (MELM-GRBF) algorithm
12
Training set: 814 Kc
61
–
96.15
Test set: 348 Kc [3]
The switching multiple models, an autoregressive hidden Markov model (AR-HMM)
–
6–8 (unclear)
50.02–74.75
–
–
[16]
Convex optimization models, linear least squares optimization model 2 (LLSOM2), WEKA classifiers
–
Training set: 30 Non-Kc
–
–
84
89.19
–
91.40
21 Kc Test set: 9 NonKc 10 Kc [18]
Support vector machines, feature extraction, feature selection
was 14.77%. Thus, the aim was to raise both Spec. and Sens. values and this can be achieved by raising ARR rate. For this reason, we aimed at obtaining highest ARR values in these experimentations. When the results were examined according to the highest values, ARR, the best Sens. rate was obtained as 95.82% for SVD&VMD method. On the other hand, the highest Spec. and ARR rates were recorded as 90.70 and 92.29% for SVD method. Combining SVD with VMD increased the Sens. rates as understood from the results, but Spec. rates decreased in this combination. The reason is that FP rates are a cumbersome in all experimented systems. FP is a general problem of all detectors for Kc detection. Sometimes, it can even be hard for sleep experts to differentiate a Kc from a Non-Kc. An illustration is given in Fig. 18. In this figure, there are some Kcs labeled by the experts as TKc (true K-complex) and one false Kc at the end of the epoch labeled as FKc (false K-complex). It is almost impossible to discriminate FKc from the TKc. This sometimes the case for sleep experts, too. Many unpracticed experts can see these FKcs as TKc in the related signal.
32
111
This study has obtained extremely high sensitivity, specificity and ARR rates according to similar studies. But, we cannot directly compare the results of these systems and methods with the other studies because of the dissimilar signals. Nonetheless, we compared the results of our experimented methods with the literature in Table 11 to give on insight. Also, used methods, datasets, subjects, the number of K-complex and obtained results in related studies are given in Table 11. Feature extraction, feature selection and various classification algorithms were usually used in literature studies. However, novel methods in this area, thresholding and logic AND gate system were utilized in this study. Therefore, this study is different from other studies in the literature. When Table 11 is briefly examined, we can say that our experimented methods performed very well in that all Kcs with their places in the dataset could be successfully found by using SVD and SVD&VMD methods. As a future study, the relationship of Kcs with other sleep parameters in EEG, EOG, EMG and ECG signals will be investigated and these correlations will also be utilized to reduce the number of FP.
123
Neural Comput & Applic
Also, other more advanced but simple methods such as DNA-based computing [28], non-stationary Gaussian process [29] and genetic algorithms can be utilized to reduce the computational complexity for detection of the threshold level in thresholding process and finding the best method. In addition to this, when researchers try to make sense from a nonlinear signal, they can consult and use the return-map or delay-map methods [30]. For instance, let x(t), t = 0,1,…, be a time series of a nonlinear signal. Then, the points (x(t), x(t ? i)), t = 0,1,…, i = 1,2,…, are the return map [30]. The return map is even more useful than the time series while the problem of pattern recognition [30].
Compliance with ethical standards Conflict of interest The authors declare that there is no conflict of interests regarding the publication of this article. Ethical standard All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national Non-invasive Clinical Research Medical Ethics Review Board and with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Informed consent Proper informed consent was obtained from all individual participants included in the study.
References 5 Conclusions Sleep staging process is an indispensable method to diagnose the sleep-related disorders. This process is implemented by using the sleep signals which are EEG, EOG, EMG, ECG and the others of the subjects. And then, the epochs are labeled as Awake, Non-REM stage 1–2–3 and REM by sleep specialists. K-complex is one of the important parameters to label epochs as N-REM2 that includes SS and K-complex. Determining of the N-REM2 within the other sleep stages is crucial for the specialists because of this stage forms about 50–60% of the full sleep time. Thanks to this study, the sleep experts can find out whether there were Kcs in related epochs and also know their locations in these epochs, automatically. Furthermore, it will help automatic sleep stage classification systems. In our study, SVD, VMD and DWT techniques and their combinations have been applied to determine the locations of Kcs in the sleep EEG. According to the results, the highest ARR was found as 92.29% for the SVD method. The results denoted that the systems and methods could be trustworthy utilized in the automatic determination the locations of Kcs. Thanks to this study at least, the sleep specialists will be able to trusty identify the Kcs in which epochs and also know their locations in the related epochs, automatically. Another important benefit of the study that the sleep staging process—exhausting, time-consuming and high margin of error for the specialists—can be implemented in less time and at a higher performance. In future, we plan to rise the ARR rates of these systems using other supporting instruments from signal processing and machine learning area. Acknowledgements This study is supported by the Scientific and Technological Research Council of Turkey (project no. 113E591) and the Scientific Research Projects Coordination Unit of Selcuk University.
123
1. Rechtschaffen A, Kales A (1969) A manual of standardized terminology. Techniques and scoring system for sleep stages of human subjects. US Government Printing Office, Washington 2. Flemons WW, Buysse D, Redline S, Pack A, Strohl K, Wheatley J, Young T, Douglas N, Levy P, McNicholas W, Fleetham J, White D, Schmidt-Nowarra W, Carley D, Romaniuk J, Force AASMT (1999) Sleep-related breathing disorders in adults: recommendations for syndrome definition and measurement techniques in clinical research. Sleep 22(5):667–689 3. Camilleri TA, Camilleri KP, Fabri SG (2014) Automatic detection of spindles and K-complexes in sleep EEG using switching multiple models. Biomed Signal Process Control 10:117–127 4. Bankman IN, Sigillito VG, Wise RA, Smith PL (1992) Featurebased detection of the K-complex wave in the human electroencephalogram using neural networks. IEEE Trans Biomed Eng 39(12):1305–1310 5. Richard C, Lengelle R (1998) Joint time and time-frequency optimal detection of K-complexes in sleep EEG. Comput Biomed Res 31(3):209–229 6. Jansen BH, Desai PR (1994) K-complex detection using multilayer perceptrons and recurrent networks. Int J Biomed Comput 37(3):249–257. doi:10.1016/0020-7101(94)90123-6 7. Bremer G, Smith JR, Karacan I (1970) Automatic detection of the K-complex in sleep electroencephalograms. IEEE Trans Biomed Eng 17(4):314–323 8. Jansen BH, Dawant BM, Meddahi K, Martens W, Griep P, Declecrk AC (1989) AI techniques for K-complex detection in human sleep EEGs. In: Proceedings of the annual international conference of the ieee engineering in engineering in medicine and biology society, 1989. Images of the twenty-first century. IEEE, pp 1806–1807 9. Henry D, Sauter D, Caspary O (1994) Comparison of detection methods: application to K-complex detection in sleep EEG. In: Engineering in medicine and biology society. Engineering advances: new opportunities for biomedical engineers. Proceedings of the 16th annual international conference of the IEEE, vol 1212, pp 1218–1219. doi:10.1109/IEMBS.1994.415401 10. Kam A, Cohen A, Geva AB, Tarasiuk A (2004) Detection of K-complexes in sleep EEG using CD-HMM. In: Engineering in medicine and biology society. IEMBS ‘04. 26th Annual international conference of the IEEE, 1–5 Sept. 2004, pp 33–36. doi:10.1109/IEMBS.2004.1403083 11. Erdamar A, Duman F, Yetkin S (2012) A wavelet and teager energy operator based method for automatic detection of K-complex in sleep EEG. Expert Syst Appl 39(1):1284–1290 12. Devuyst S, Dutoit T, Stenuit P, Kerkhofs M (2010) Automatic K-complexes detection in sleep EEG recordings using likelihood thresholds. In: 2010 Annual international conference of the ieee
Neural Comput & Applic
13.
14.
15.
16.
17.
18.
19. 20.
21.
engineering in medicine and biology, Aug. 31 2010–Sept. 4 2010, pp 4658–4661. doi:10.1109/IEMBS.2010.5626447 Zacharaki EI, Pippa E, Koupparis A, Kokkinos V, Kostopoulos GK, Megalooikonomou V (2013) One-class classification of temporal EEG patterns for K-complex extraction. Conf Proc IEEE Eng Med Biol Soc 2013:5801–5804. doi:10.1109/EMBC. 2013.6610870 Krohne LK, Hansen RB, Christensen JA, Sorensen HB, Jennum P (2014) Detection of K-complexes based on the wavelet transform. Conf Proc IEEE Eng Med Biol Soc 2014:5450–5453. doi:10.1109/EMBC.2014.6944859 Noori SMR, Hekmatmanesh A, Mikaeili M, SadeghniiatHaghighi K (2014) K-complex identification in sleep EEG using MELM-GRBF classifier. In: 2014 21th Iranian conference on biomedical engineering (ICBME), 26–28 Nov 2014, pp 119–123. doi:10.1109/ICBME.2014.7043905 Zamir ZR, Sukhorukova N, Amiel H, Ugon A, Philippe C (2015) Convex optimisation-based methods for K-complex detection. Appl Math Comput 268:947–956 Parekh A, Selesnick IW, Rapoport DM, Ayappa I (2015) Detection of K-complexes and sleep spindles (DETOKS) using sparse optimization. J Neurosci Methods 251:37–46 Hernandez-Pereira E, Bolon-Canedo V, Sanchez-Marono N, Alvarez-Estevez D, Moret-Bonillo V, Alonso-Betanzos A (2016) A comparison of performance of K-complex classification methods using feature selection. Inf Sci 328:1–14 Pan J, Tompkins WJ (1985) A real-time QRS detection algorithm. IEEE Trans Biomed Eng 32(3):230–236 Lee S, Hayes MH (2004) Properties of the singular value decomposition for efficient data clustering. IEEE Signal Proc Lett 11(11):862–866 Xia Y, Zhou W, Li C, Yuan Q, Geng S (2015) Seizure detection approach using S-transform and singular value decomposition. Epilepsy Behav 52:187–193
22. Hassanpour H, Mesbah M, Boashash B (2004) Time-frequency feature extraction of newborn EEG seizure using SVD-based techniques. EURASIP J Adv Signal Process 16:1–11. doi:10. 1155/s1110865704406167 23. Huang NE, Shen Z, Long SR, Wu MLC, Shih HH, Zheng QN, Yen NC, Tung CC, Liu HH (1998) The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis. Proc R Soc A Math Phys 454(1971):903–995 24. Dragomiretskiy K, Zosso D (2014) Variational mode decomposition. IEEE Trans Signal Process 62(3):531–544 25. Liu YY, Yang GL, Li M, Yin HL (2016) Variational mode decomposition denoising combined the detrended fluctuation analysis. Signal Process 125:349–364 26. Ocak H (2009) Automatic detection of epileptic seizures in EEG using discrete wavelet transform and approximate entropy. Expert Syst Appl 36(2):2027–2036 27. Mallat SG (1989) A theory for multiresolution signal decomposition—the wavelet representation. IEEE Trans Pattern Anal 11(7):674–693 28. Ullah AMMS (2010) A DNA-based computing method for solving control chart pattern recognition problems. CIRP J Manuf Sci Technol 3(4):293–303. doi:10.1016/j.cirpj.2011.02.002 29. Ullah AMMS, Harib KH (2010) Simulation of cutting force using nonstationary Gaussian process. J Intell Manuf 21(6):681–691. doi:10.1007/s10845-009-0245-2 30. Ullah A, Harib KH (2006) Knowledge extraction from time series and its application to surface roughness simulation. Inf Knowl Syst Manag 5(2):117–134
123
