Mixed-phase modeling in snore sound analysis

Abstract

Obstructive sleep apnea (OSA) is a highly prevalent disease in which upper airways are collapsed during sleep, leading to serious consequences. The gold standard of diagnosis, called Polysomnography (PSG), requires a full-night hospital stay connected to over 15 channels of measurements requiring physical contact with sensors. PSG is expensive and unsuited for community screening. Snoring is the earliest symptom of OSA, but its potential in OSA diagnosis is not fully recognized yet. In this paper, we propose a novel model for SRS as the response of a mixed-phase system (total airways response, TAR) to a source excitation at the input. The TAR/source model is similar to the vocal tract/source model in speech synthesis, and is capable of capturing acoustical changes brought about by the collapsing upper airways in OSA. We propose an algorithm based on higher-order-spectra (HOS) to jointly estimate the source and TAR, preserving the true phase characteristics of the latter. Working on a clinical database of signals, we show that TAR is indeed a mixed-phased signal and second-order statistics cannot fully characterize it. Night-time speech sounds can corrupt snore recordings and pose a challenge to snore based OSA diagnosis. We show that the TAR could be used to detect speech segments embedded in snores, and derive features to diagnose OSA via non-contact, low-cost instrumentation holding potential for a community screening device.

Appendices

Appendix A

Appendix A describes a novel method to jointly estimate the pitch period and the TAR (or VTR) using an exhaustive search process.

The approach is to assume a set of test-values (q) for the pitch of x _v (n) and then determine the best value q _opt of q based on a performance measure (see Eq. (14)). The quantity q _opt is then considered as the pitch of x _v (n). Steps (S1)–(S8) below describe the procedure (Please see Fig. 15 for a block diagram).

• (S1) Initialization: Set q = 1; set the window parameter l = 1.

• (S2) Estimate ${\tilde{h}}_{v} (n)$ via (5)–(8) and (9)–(12) for the assumed values of q and l; generate a test data block ${\tilde{s}}_{v} (n)$ via ${\tilde{s}}_{v} (n) = {\tilde{x}}_{v} (n)* {\tilde{h}}_{v} (n),$ where ${\tilde{x}}_{v} (n)$ is the unit impulse train with the inter-impulse pitch set to the assumed q.

• (S3) Normalize the original data block s _v (n) and the test block ${\tilde{s}}_{v} (n)$ each to have rms values equal to 1; let the Fast Fourier Transform (FFT) of the blocks, respectively be S _v (j) and ${\tilde{S}}_{v} (j), \quad j=0,1,2,\ldots, N-1,$ where N is the length of the FFT.

• (S4) Calculate the distance measure e(l,q) defined as:

  $$ e(l,q) = {\sqrt \frac{\sum_{k = 0}^{N - 1} {\left({\vert S_{v} (j) \vert - \vert {\tilde{S}}_{v} (j) \vert} \right)}^{2}} {N}}$$

  (14)

• (S5) Repeat steps (S2)–(S4) increasing q to q + 1 at each repetition until q reaches an arbitrarily large number M _q . Write each e(l,q),(q = 1,2, ... M _q ) in the error matrix E[ ] M _lx M _q row l column q,  (q = 1,2, ... M _q ).

• (S6) Repeat steps (S1)–(S5) increasing l to l +1 at each repetition until l reaches an arbitrarily large number M _l . Write each e(l, q), (q = 1,2, ... M _q ) in the error matrix E[ ] M _lxM _q row l,(l = 1, 2, ... M _l ).

• (S7) Find the global minimum entry, e _min(l,q) of E[ ] M _lxM _q .

• (S8) Define a threshold value ɛ for e _min(l,q). If e _min(l,q) < ɛ, the data block is periodic and thus contains a voiced-segment. Obtain the corresponding value of q as the optimum pitch value q _opt. Use q _opt as the pitch of source excitation x _v (n), and consider the corresponding ${\tilde{h}}_{v} (n)$ to be the TAR. If e _min(l,q) > ɛ the data block is treated as from an unvoiced-segment, and no pitch value is recorded.

Fig. 15

Block diagram of the exhaustive search process in Appendix A

Appendix B

1.1 Minimum-phase/all-pass decomposition of the TAR

Any mixed-phase signal can be represented as the convolution between a unique minimum-phase equivalent and an all-pass component [17]. For the reason ${\tilde{h}}(n)$ is modeled as a mixed phase sequence and estimated preserving the true phase characteristics, it is amenable for the min-phase/all-pass decomposition as given by:

$$ {\tilde{h}}(n) = {\tilde{h}}_{m} (n)* {\tilde{h}}_{a} (n). $$

(15)

In (15) ${\tilde{h}}_{m} (n)$ is the minimum phase equivalent of ${\tilde{h}}(n)$ such that their Fourier transform magnitudes, |H(e ^jω)| = |H _m (e ^jω)| and ${\tilde{h}}_{a} (n)$ is the all pass sequence of ${\tilde{h}}(n).$

Then,

$$ {\left\vert {H_{a} (e^{{j\omega}})} \right\vert} = \frac{{{\left\vert {H(e^{{j\omega}})} \right\vert}}}{{{\left\vert {H_{m} (e^{{j\omega}})} \right\vert}}} = 1, $$

(16)

And

$$ \angle H_{a} (e^{{j\omega}}) = \angle H(e^{{j\omega}}) - \angle H_{m} (e^{{j\omega}}), $$

(17)

where H(e ^jω),   H _m (e ^jω) and H _a (e ^jω) are the Fourier transforms of ${\tilde{h}}(n),\, {\tilde{h}}_{m} (n)$ and ${\tilde{h}}_{a}(n),$ respectively.

Note that in the special case when ${\tilde{h}}(n)$ is a pure minimum phase signal |H _a (e ^jω)| = 1 and $\angle H_{a} (e^{{j\omega}}) = 0.$

Figure 16 shows the process of decomposition of ${\tilde{h}}(n)$ into the minimum-phase and all pass components.

Fig. 16

Minimum-phase/all-pass decomposition of the TAR

The real cepstrum ${\hat{c}}_{h} (n)$ of ${\tilde{h}}(n),$ a nonlinear operator used in speech analysis, is frequency-invariant linear filtered [17] to obtain the complex cepstrum, ${\hat{h}}_{m} (n)$ of the minimum phase sequence, ${\tilde{h}}_{m} (n),$ via;

$$ {\hat{h}}_{m} (n) = {\hat{c}}_{h} (n)\,w_{m} (n), $$

(18)

Where w _m (n) = 2u(n) − δ(n).

Now we obtain the minimum phase equivalent signal ${\tilde{h}}_{m} (n)$ of ${\tilde{h}}(n)$ from the inverse complex operation in (11).

Transforming (18) in to complex cepstrum domain where convolution becomes an addition, followed by a simple algebraic operation, the complex cepstrum ${\hat{h}}_{a} (n)$ of the all pass sequence ${\tilde{h}}_{a} (n)$ is obtained as:

$$ {\hat{h}}_{a} (n) = {\hat{h}}(n) - {\hat{h}}_{m} (n). $$

(19)

The all pass sequence ${\tilde{h}}_{a} (n),$ can now be obtained from the inverse complex operation defined in (11).