It is worth noting that for snoring sound analysis, we used 80 ms time windows (M = 820) with 50% overlap to ensure the stationarity assumption.

2) Calculate the zero-mean signal of each segment as:

(9)

3) Multiply the zero-mean signal by the Hanning window, , to control the effect of spectral leakage.

(10)

4) Compute the discrete Fourier transform (DFT) of each segment:

(11)

The raw bispectral estimate () can be calculated as:

(12)

where l, m are the discrete frequencies.

5) The consistent estimate of bispectrum () can be obtained by averaging raw estimates over all segments.

(13)

Consequently, the squared bicoherence can be derived from bispectrum as below [18]:

(14)

The discrete bispectrum has many symmetries in plane. It is only needed to calculate in the non-redundant region or principal domain (D) which is defined as:

[19]. Where are the bifrequencies (in Hertz) correspondent to the normalized bifrequencies.

2.3. Feature Extraction

Suppose that we estimated the bispectrum () in D. This section details on deriving two new features defined in section I: 1) The MBF feature, which is a 2-D feature denoted as and 2) PMBF, which is a 1-D feature denoted as.

2.3.1. MBF Computation

MBF is the bifrequency where the summation of absolute values of becomes half of the summation of absolute value of over all bifrequencies in D. In fact, the procedure looks like:

1) Calculate the summation of at all bifrequencies in D.

2) Set f_{1} = 0.

3) For all bifrequencies satisfying the condition calculate:

4) Check if

If YES, end the algorithm and

If NO, increase f_{1} and go to step 3. (Note that: .)

2.3.2. PMBF Computation

Once the MBF is computed, the PMBF, , can be determined by the projection of onto the line corresponding to the diagonal slice of the bispectrum. Equivalently we have:

2.3.3. Skewness and Kurtosis

Let be a zero-mean random process. Skewness and kurtosis are defined as:

where is the standard deviation of and and are its zero-lag 3^{rd} and 4^{th} order cumulants respectively [20].

2.3.4. 1^{st} Formant Frequency and Energy

Energy (E) and first formant frequency (F_{1}) were obtained from each SS segment. Linear predictive coding (LPC) [21] was used to estimate F_{1}. To meet stationarity assumption, was divided into 80 ms overlapping frames (50% overlap and Hanning window). In each frame, the autoregressive (AR) model of the signal was estimated and the roots of AR model were calculated. To select the AR model order, we used the optimum order model (optimum order = ) suggested in [22]. Therefore, we selected an AR model of 14 to estimate first formant frequency of each frame. F_{1} was estimated by taking median over all frames.

2.3.5. Calculation of Features

As mentioned in the Section 2.1, the number of SS segments is different for each patient. Let us denote SS segment of patient X by First, , , , , for all segment were calculated resulting in a finite number of observations for each feature. Then, the sample median of each feature set was estimated. In fact the followings were calculated for patient X:

This procedure was repeated for all 60 individuals yielding a matrix for MBF feature and a vector for PMBF, energy, 1^{st} formant frequency, skewness, and kurtosis features. The reason we used median instead of mean is the insensitivity of median to outliers; it is known that when the data is not symmetrically distributed, the median outperforms the mean in measuring the middle range of data [23]. Figure 2 shows the sample density of, and estimated by kernel method [24]. As shown in Figure 2, the data is skewed; therefore, the median is a better estimate of the middle range of the data than the mean value in this case.

2.4. Statistical Analysis

To investigate the effect of anthropometric parameters such as age, gender, height, BMI, and AHI on the features, we ran statistical tests assuming the significance level as. Since the distribution of the features deviated from normal distribution, the Kendall’s Tau-b test (nonparametric counterpart of Pearson correlation) [10] was used to measure the correlation among continuous anthropometric parameters and HOS features. The one-way KWAV [11,12] was also used to compare the median of features between men and women.

2.5. Classification

Naïve Bayes classifier [13] was used to evaluate the ability of our feature set to discriminate the subjects to snorers with OSA and snorers without OSA or the so called “simple snorer” groups. Particularly, we were interested to compare the ability of SS features to be used as a signature of OSA when the groups of OSA and non-OSA were matched (Experiment A) and unmatched (Experiment B) in terms of anthropometric parameters.

Therefore, we performed two experiments: Experiment A: We selected a subset of our database including 22 apneic and 6 simple snorers that were matched in terms of gender, BMI, height and AHI. Experiment B: Another subset with the same number of participants (28 including 21 apneic and 7 simple snorers) with unmatched anthropometric parameters was used for classification. Table 2 shows the anthropometric profile of both experiments.

We used the energy, 1^{st} formant frequency, MBF, PMBF, skewness, and kurtosis as our features for linear discriminant analysis. Several combinations of the features were examined and the performance was evaluated using the Leave-One-Out Cross-Validation (LOOCV) technique [25,26]. The LOOCV is a common technique when the number of observations (subjects in this case) is relatively small; it helps to prevent over-fitting. In the LOOCV, one observation is used as testing set and the rest [27] is used as training set. This procedure is repeated for all observations [28] and the average performance is measured. It is worth noting that the Euclidean metric was used to compute the distance.

Figure 2. Kernel density estimate of, , and for a typical subject. Note the asymmetry of the distributions.

3. RESULTS AND DISCUSSION

All SS segments were found to be non-Gausssian, while their linearity varied during the night. In fact, for each patient, there existed some linear SS segments as well as some non-linear ones. We also noticed that the linearity of the SS segments varied among different body positions within each subject. However, this result was not consistent among all subjects. Furthermore, not everybody slept and/or snored in all positions.

It was shown that the body position during sleep changes both duration and intensity of snoring sounds [27]. However, we did not find a consistent and significant relationship between the sleeping position and the change in the linearity of snoring sound generating process.

However, in general, we did not find a consistent and significant relationship between the sleeping position and the change in the linearity of snoring sound generating process; the results in the population of our study were variable between the subjects. However, the lack of finding a general pattern of such relationship could be partially due to lack of snoring data in all different sleeping positions in our study.

It is known that if a signal is non-Gaussian, the 2^{nd} order statistical techniques are only able to extract partial information from the signal [7,8]. Therefore, we used HOS measures to develop new features such as MBF and PMBF from existing data. We also extracted common HOS features such as skewness and kurtosis from the SS segments. It was found that there was a significant relationship between frequency based features such as, , and and all anthropometric parameters except age. As shown in Table 3, four out of five anthropometric parameters (height, BMI, AHI, and gender) significantly affected the HOS features of the SS segments, while gender and BMI were significant parameters affecting energy and first formant features.

The height of individuals was observed to be a significant factor influencing the value of (p < 0.01), , and (p < 0.05). There was a negative relationship between height and these frequency related features. The taller the individuals, the lower frequency components were in their snoring bispectrum. The height has been shown to affect the tracheal sound spectral features [28]. It was reported that the tracheal sounds in children had higher frequency components than in healthy

Table 2. Anthropometric information of two subsets selected for classification.

Table 3. The results of Kendall’s Tau-b and Kruskal-Wallis tests for five anthropometric parameters.

adults. In another study [29], it is shown that the anatomy of the trachea determines the characteristic features of tracheal sounds. However, there was no study confirming the change in the features of SS segments due to the height. Based on our findings, the MBF and PMBF features of the extracted SS segments are negatively related to the height of individuals. Assuming that taller individuals have taller neck, this result implies that the characteristics of SS segments reflect resonances (existing in SS) that depend on the upper airway’s length.

The results of the Kendall’s Tau-b test on BMI groups shows that BMI is significantly associated with the value of (p < 0.05), (p < 0.01) as well as F_{1} (p < 0.05). They were negatively correlated meaning that the higher the BMI, the lower were the values of, , and F_{1}. As known, obesity is a factor strongly associated with the presence of OSA [30]. Obese individuals with sleep apnea have been shown to have more (about 42%) fat in their cervical region than normal subjects as well as non-obese individuals with OSA [31]; thus, resulting in pharyngeal area narrowing. It is also known that higher BMI is associated with increased level of leptin (a hormone produced by the adipose tissue and has also actions on the respiratory centre control) [32]. Therefore, our observed changes in the acoustical properties of the SS segments due to BMI can be explained by both anatomical and hormonal changes of the upper airway.

It was also found that AHI and gender were significantly correlated with energy and frequency-based HOS features of the SS segments. As shown in Table 3, the individuals with higher AHI had lower frequency-based features (and) (p < 0.05) and higher energy (p < 0.01). The female snorers of this study were observed to have higher frequency-based features (and) (p < 0.05) and lower energy feature (p < 0.05) than the male snorers. Although there was no study investigating the gender effect on the snoring sounds, this observation is congruent with findings reported in two studies focused on breath and lung sounds [33,34]. According to those studies, breath and lung sounds in healthy women contain higher frequency components than in men. It has also been shown that men have higher pharyngeal and supraglottic resistances than women [35]. Therefore, given that the size and mechanical properties of pharynx are significantly different between men and women [36], the snoring sounds of women and men can be expected to be significantly different as the results of our study indicate. Moreover, these might be also a reason for greater incidence of OSA in men [35,36].

Two of the frequency-based HOS features (and) were found to be significantly different in snorers with different AHI. This result is congruent with previous studies. In people with OSA, the lateral pharyngeal muscular wall is usually narrower [37]. Therefore, minimum area of the airway has been shown to be significantly smaller in apneic individuals than non-OSA people. The size of airway plays a major role in the frequency components of the sound produced by the flow turbulence in the airway. This explains the change in the frequency based HOS feature of the SS segments between snorers with OSA and simple snorers.

It was observed that for some of the anthropometric parameters (e.g. BMI and gender) two of the frequency-based features (, , and) were significantly correlated, while the third one was not significantly correlated. is linearly related to the summation of and (projection of two). If both have a significant correlation with a parameter, then we expect that would be also significant (as in the case with Height parameter) but having one of them significantly correlated with a parameter, does not necessarily lead to a significant correlation of and that parameter. The reason that only one of the coordinates of MBF is significant depends on the bispectrum of the SS segment. As an example, let us compare MBF for 4 SS segments of participant4 (P4) and participant6 (P6). P4 (BMI = 24.4): = [320,320,360,320] and = [240,40,160,160]. P6 (BMI = 47.1): = [240,240, 200,250] and = [180,160,180,120]. It is clear that BMI significantly changed but not. In fact, the difference between the bispectral information of the two sets of SS segments is well extracted using 2-D MBF feature which is an advantage of bispectral analysis.

One important point is that these frequency changes due to small changes in the airway size may not always be detectable by spectral analysis of the sounds. However, as known, HOS techniques complement the information obtained from 2^{nd} order statistical techniques, i.e. power spectral analysis. Hence, we propose using a combination of HOS techniques and conventional acoustical techniques increases the diagnosis accuracy of OSA. In this paper, we tried to verify this point by applying a simple classifier to our feature set. We partitioned our database into two sets to compare two scenarios, one when the height, gender, and BMI are matched between the two groups of snorers with OSA and simple snorers, and the other one when those parameters are not matched. We observed an increase in the accuracy of classification when the parameters were matched.

Table 4 illustrates the results of classification. Several combinations of features were used as input to the naïve Bayes classifier. Results demonstrate that the highest sensitivity and specificity occurred when a combination of both conventional feature (Energy) and HOS feature set (and skewness) was used. This combination resulted in sensitivity of 93.2% (87.5%) and the specificity of 88.4% (86.3%) for experiment A (B). As shown in Table 4, for experiment A, the sensitivity and specificity

Table 4. Naïve Bayes classification results for different combination of conventional and HOS features.

values for only HOS features were 75.9% - 94.1% and 74.6% - 81.9%, respectively. On the other hand, using only energy and formant frequency resulted in a sensitivity and specificity of 78.2% and 72.1%, respectively.

As expected, overall, the sensitivity and specificity decreased when an unmatched subset was used for classification.

To compare our work with a recently published work [9], we matched the anthropometric parameters of snorers with OSA and simple snorer groups. Moreover, our recordings were performed using a microphone placed over trachea. Therefore, our recorded sounds have a higher signal to noise ratio than those recorded by an ambient microphone. We also improved the sensitivity and specificity of OSA diagnosis among snorers by simultaneous usage of HOS and conventional features. However a major difficulty in our study was to find a larger population with matched anthropometric parameters to validate the results of our analysis. This limitation was partially resolved by using LOOCV which performs well when the population size is relatively small.

4. CONCLUSION

In this study, the relationship between anthropometric parameters of 60 snorers and the 3^{rd} and 4^{th} order statistical features derived from the SS segments were investigated.

In summary, we developed two new frequency-based HOS features from the non-Gaussian SS segments, and investigated statistical correlation of these features along with the zero-lag HOS features with different anthropometric parameters. An important contribution of the statistical investigation is on the application of snoring sound for OSA identification among snorers. Since the common features of snoring sounds used in classification are sensitive to anthropometric parameters, the results of classification is reliable only when the two groups of apneic and controls are matched for those parameters.

When HOS features were used in the classification of the apneic group, the results have shown improvement compared to those of previous studies. The HOS features may also be used (under investigation) to find the site of obstruction in upper airway, and cluster the SS segments based on their production mechanism.

5. ACKNOWLEDGEMENTS

In summary, we developed two new frequency-based HOS features from the non-Gaussian SS segments, and investigated statistical correlation of these features along with the zero-lag HOS features with different anthropometric parameters. An important contribution of the statistical investigation is on the application of snoring sound for OSA identification among snorers. Since the common features of snoring sounds used in classification are sensitive to anthropometric parameters, the results of classification is reliable only when the two groups of apneic and controls are matched for those parameters.

When HOS features were used in the classification of the apneic group, the results have shown improvement compared to those of previous studies. The HOS features may also be used (under investigation) to find the site of obstruction in upper airway, and cluster the SS segments based on their production mechanism.

REFERENCES

- Hoffstein, V. (2002) Apnea and snoring: State of the art and future directions. Acta Otorhinolaryngol Belg, 56, 205-236.
- Lee, B., Hill, P., Osborne, J. and Osman, E. (1999) A simple audio data logger for objective assessment of snoring in the home. Physiological Measurement, 20, 119-127. doi:10.1088/0967-3334/20/2/001
- Dalmasso, F. and Prota, R. (1996) Snoring: Analysis, measurement, clinical implications and applications. European Respiratory Journal, 9, 146-159. doi:10.1183/09031936.96.09010146
- Sola-Soler, J., Jane, R., Fiz, J. and Morera, J. (2003) Spectral envelope analysis in snoring signals from simple snorers and patients with obstructive sleep apnea. IEEEEMBS, Cancun, Mexico, 2527-2530.
- Abeyratne, U., Karunajeewa, A. and Hukins, C. (2007) Mixed-phase modeling in snore sound analysis. Medical and Biological Engineering and Computing, 45, 791-806. doi:10.1007/s11517-007-0186-x
- Beck, R., Odeh, M., Oliven, A. and Gavriely, N. (1995) The acoustic properties of snores. European Respiratory Journal, 8, 2120-2128. doi:10.1183/09031936.95.08122120
- Fackrell. J.W.A. (1996) Bispectral analysis of speech signals. Edinburgh Research Archive, University of Edinburgh.
- Ng, A.K., Wong, K.Y., Tan, C.H. and Koh, T.S. (2007) Bispectral analysis of snore signals for obstructive sleep apnea detection. Proceedings of the 29th Annual International Conference of the IEEE EMBS, 6195-6198.
- Ng, A., Koh, T., Abeyratne, U. and Puvanendran, K. (2009) Investigation of obstructive sleep apnea using nonlinear mode interactions in nonstationary snore signals. Annals of Biomedical Engineering, 37, 1796-1806. doi:10.1007/s10439-009-9744-8
- Kendall, M.G. (1938) A New Measure of Rank Correlation. Biometrika, 30, 81-93. doi:10.2307/2332226
- Kvam, P.H. and Vidakovic, B. (2007) Nonparametric statistics with applications to science and engineering. Wiley-Interscience. doi:10.1002/9780470168707
- Searle, S.R. (1971) Linear models. John Wiley & Sons, Inc, New York.
- Duda, R.O., Hart, P.E. and Stork, D.G. (2000) Pattern Classification. 2nd Edition, Wiley-Interscience.
- Yadollahi, A. and Moussavi, Z. (2009) Acoustic obstructive sleep apnea detection. Engineering in Medicine and Biology Society, EMBC 2009. Annual International Conference of the IEEE, 7110-7113.
- Azarbarzin, A. and Moussavi, Z. (2011) Automatic and unsupervised snore sound extraction from respiratory sound signals. IEEE Transactions on Biomedical Engineering, 58, 1156-1162. doi:10.1109/TBME.2010.2061846
- Brillinger, D.R. (1965) An introduction to polyspectra. The Annals of Mathematical Statistics, 36, 1351-1374. doi:10.1214/aoms/1177699896
- Hinich, M.J. (1982) Testing for gaussianity and linearity of a stationary time series. Journal of Time Series Analysis, 3, 169-176. doi:10.1111/j.1467-9892.1982.tb00339.x
- Brillinger, D.R. and Rosenblatt, M. (1967) Computation and interpretation of k-th order spectra. Spectral Analysis of Time Series, 189-232.
- Hinich, M.J. and Wolinsky, M.A. (1988) A test for aliasing using bispectral analysis. Journal of the American Statistical Association, 83, 499-502.
- Swami, A., Mendel, J.M. and Nikias, C.L. (1998) Higherorder spectral analysis toolbox user’s guide. Version 2.
- Proakis, J.G. and Manolakis, D.K. (2006) Digital signal processing. 4th Edition, Prentice Hall, New York.
- Markel, J. (1972) Digital inverse filtering—A new tool for formant trajectory estimation. IEEE Transactions on Audio and Electroacoustics, 20, 129-137. doi:10.1109/TAU.1972.1162367
- Hogg, R., Craig, A. and Mckean, J. (2004) Introduction to mathematical statistics. Prentice Hall, New York.
- Silverman, B.W. (1986) Density estimation for statistics and data analysis. Chapman and Hall, London.
- Burman, P. (1989) A comparative study of ordinary cross-validation, v-fold cross-validation and the repeated learning-testing methods. Biometrika, 76, 503-514.
- Venter, J.H. and Snyman, J.L.J. (1995) A note on the generalised Cross-Validation Criterion in Linear Model Selection. Biometrika, 82, 215-219. doi:10.1093/biomet/82.1.215
- Nakano, H., Ikeda, T., Hayashi, M., Ohshima, E. and Onizuka, A. (2003) Effects of body position on snoring in apneic and nonapneic snorers. Sleep, 26, 169-172.
- Sanchez, I. and Pasterkamp, H. (1993) Tracheal sound spectra depend on body height. American Journal of Respiratory and Critical Care Medicine, 148, 1083-1087. doi:10.1164/ajrccm/148.4_Pt_1.1083
- Kraman, S., Pasterkamp, H., Kompis, M., Takase, M. and Wodicka, G. (1998) Effects of breathing pathways on tracheal sound spectral features. Respiration Physiology, 111, 295-300. doi:10.1016/S0034-5687(97)00113-8
- Young, T., Palta, M., Dempsey, J., Skatrud, J., Weber, S. and Badr, S. (1993) The occurrence of sleep-disordered breathing among middle-aged adults. New England Journal of Medicine, 328, 1230-1235. doi:10.1056/NEJM199304293281704
- Mortimore, I.L., Marshall, I., Wraith, P.K., Sellar, R.J. and Douglas, N.J. (1998) Neck and total body fat deposition in nonobese and obese patients with sleep apnea compared with that in control subjects. American Journal of Respiratory and Critical Care Medicine, 157, 280-283.
- De Sousa, A.G., Cercato, C., Mancini, M.C. and Halpern, A. (2008) Obesity and obstructive sleep apnea-hypopnea syndrome. Obesity Reviews, 9, 340-354. doi:10.1111/j.1467-789X.2008.00478.x
- Gavriely, N., Nissan, M., Rubin, A.H. and Cugell, D.W. (1995) Spectral characteristics of chest wall breath sounds in normal subjects. Thorax, 50, 1292-1300. doi:10.1136/thx.50.12.1292
- Gross, V., Dittmar, A., Penzel, T., Schuttler, F. and von Wichert, P. (2000) The Relationship between normal lung sounds, age, and gender. American Journal of Respiratory and Critical Care Medicine, 162, 905-909.
- White, D.P., Lombard, R.M., Cadieux, R.J. and Zwillich, C.W. (1985) Pharyngeal resistance in normal humans: influence of gender, age, and obesity. Journal of Applied Physiology, 58, 365-371.
- Brooks, L.J. and Strohl, K.P. (1992) Size and mechanical properties of the pharynx in healthy men and women. American Journal of Respiratory and Critical Care Medicine, 146, 1394-1397.
- Schwab, R., Gupta, K., Gefter, W., Metzger, L., Hoffman, E. and Pack, A. (1995) Upper airway and soft tissue anatomy in normal subjects and patients with sleep-disordered breathing. Significance of the lateral pharyngeal walls. American Journal of Respiratory and Critical Care Medicine, 152, 1673-1689.