METHOD FOR DETECTING AND DISCRIMINATING BREATHING PATTERNS FROM RESPIRATORY SIGNALS

Abstract

A Cheyne-Stokes (CS) diagnosis system classifies periods of CS-like breathing by examining a signal indicative of a respiratory parameter. For example, nasal flow data is processed to classify it as unambiguously CS breathing or nearly so and to display the classification Processing may detect and display: apnoeas, hypopnoeas, flow-limitation and snore. The signal may be split into equal length epochs and event features are extracted. Statistics are applied to these primary feature(s) to produce secondary feature(s) representing the entire epoch. Each secondary feature is grouped with other feature(s) extracted from the entire epoch rather than from the epoch events. This final group of features is the epoch pattern. The epoch pattern is classified to produce a probability for possible event classes (e.g., Cheyne-Stokes breathing, OSA, etc.). The epoch is assigned to the class with the highest probability, which may both be reported as an indication of disease state.

Claims

1. A method for a controlled ventilatory apparatus to automatically determine an absence or a presence of Cheyne-Stokes breathing in a person comprising: accessing a signal indicative of respiration of a person derived from a sensor; repeatedly in a processor, determining a presence of an apnea or an hypopnea from said signal and, for the determined presence of an apnea or hypopnea, calculating a maximum increase from multiple increases of said signal of a hyperpnoea following the determined presence of an apnea or an hypopnea, wherein each of the multiple increases is determined from a predetermined time interval in a period of the hyperpnoea between its beginning and maximum, wherein a plurality of said maximum increases are determined from a plurality of hyperpnoeas, the plurality of maximum increases including a first maximum increase and a second maximum increase such that the first maximum increase is from one hyperpnoea and the second maximum increase is from another hyperpnoea; determining any peak in a frequency spectrum of said signal; and in the processor, determining the absence of Cheyne-Stokes breathing based on the first maximum increase, and determining the presence of Cheyne-Stokes breathing based on a combination of a peak in the frequency spectrum and the second maximum increase, the first maximum increase being greater than the second maximum increase such that the first maximum increase characterizes a hyperpnoea typical of obstructive sleep apnea and the second maximum increase characterizes a hyperpnoea typical of Cheyne-Stokes breathing, wherein the processor is configured to operate as a classifier with a discriminant function that generates a probability estimate for a class representing presence of Cheyne-Stokes breathing and a class representing absence of Cheyne-Stokes breathing, and generates output that indicates the determined absence or presence of Cheyne-Stokes breathing based on a highest probability from the discriminant function.

2. The method of claim 1 wherein said signal indicative of respiration is a ventilation signal.

3. The method of claim 1 wherein said signal indicative of respiration comprises oxygen saturation.

4. The method of claim 1 wherein the determining any peak in a frequency spectrum of said signal comprises calculating a spectrogram of said signal and determining any peak from said spectrogram.

5. The method of claim 1 further comprising determining whether the signal has a peak in a frequency range of 0 Hz to 0.075 Hz.

6. The method of claim 1 wherein the determining said peak in a frequency spectrum of said signal includes calculating a Fourier transform of the signal.

7. The method of claim 1 further comprising calculating a shape of said signal and utilizing said shape as an indication of central apnea.

8. The method of claim 1 wherein the processor is configured to generate the output as a display on an epoch-by-epoch basis that represents the determined absence and/or presence of Cheyne-Stokes breathing.

9. The method of claim 2 further comprising processing said ventilation signal to derive an envelope of said ventilation signal during a hyperpnoea.

10. The method of claim 9 wherein the envelope is derived by a droopy-peak detector.

11. The method of claim 10 further comprising interpolating the envelope.

12. The method of claim 11 wherein the envelope comprises the predetermined time intervals.

13. The method of claim 9 further comprising scaling a maximum difference value by a mean of the ventilation signal.

14. The method of claim 1, wherein the predetermined time intervals are on the order of time of a breath.

15. The method of claim 1, wherein the predetermined time interval is a two second time interval.

16. The method of claim 2, wherein the ventilation signal is an absolute value of a nasal flow signal.

17. A ventilatory apparatus for automatically detecting an absence or a presence of Cheyne-Stokes breathing in a person comprising: a sensor for sensing a signal indicative of respiration of a person; a detector for repeatedly determining a presence of apnea or hypopnea in said signal; and a processor configured to: (a) determine any peak in a frequency spectrum of said signal, (b) repeatedly, for the determined presence of apnea or hypopnea in said signal, calculate a maximum increase from multiple increases in said signal with respect to time of a hyperpnoea following the determined apnea or hypopnea, wherein each of the multiple increases in said signal with respect to time is determined from a predetermined time interval in a period of the hyperpnoea between its beginning and its maximum, wherein a plurality of said maximum increases are determined from a plurality of hyperpnoeas, the plurality of maximum increases including a first maximum increase and a second maximum increase such that the first maximum increase is from one hyperpnoea and the second maximum increase is from another hyperpnoea, and (c) determine the absence of Cheyne-Stokes breathing from the first maximum increase, and the presence of Cheyne-Stokes breathing from a combination of a peak in the frequency spectrum and the second maximum increase, the first maximum increase being greater than the second maximum increase such that the first maximum increase characterizes a hyperpnoea typical of obstructive sleep apnea and the second maximum increase characterizes a hyperpnoea typical of Cheyne-Stokes breathing, wherein the processor is configured to operate as a classifier with a discriminant function that generates a probability estimate for a class representing presence of Cheyne-Stokes breathing and a class representing absence of Cheyne-Stokes breathing, and to generate output that indicates the determined absence or presence of Cheyne-Stokes breathing based on a highest probability from the discriminant function.

18. The apparatus of claim 17 wherein said signal indicative of respiration is a ventilation signal.

19. The apparatus of claim 17 wherein said signal indicative of respiration comprises oxygen saturation detected by a pulse oximeter.

20. The apparatus of claim 17 wherein said processor is configured to determine a peak in a frequency spectrum of said signal by calculating a spectrogram of said signal.

21. The apparatus of claim 17 wherein said processor is configured to determine whether the signal has a peak in a frequency range of 0 Hz to 0.075 Hz.

22. The apparatus of claim 17 wherein said processor is configured to determine a peak in a frequency spectrum of said signal by calculating a Fourier transform of said signal.

23. The apparatus of claim 17 wherein said processor is configured to calculate a shape of said signal and utilize said shape as an indication of central apnea.

24. The apparatus of claim 18 wherein the processor is further configured to process said ventilation signal to derive an envelope of said ventilation signal during a hyperpnoea.

25. The apparatus of claim 24 wherein the envelope is derived by a droopy-peak detector of the processor.

26. The apparatus of claim 25 wherein the processor is configured to interpolate the envelope.

27. The apparatus of claim 26 wherein the envelope comprises the predetermined time interval.

28. The apparatus of claim 24 wherein the processor is further configured to scale a maximum difference value by a mean of the ventilation signal.

29. The apparatus of claim 17, wherein the predetermined time intervals are on the order of time of a breath.

30. The apparatus of claim 17, wherein the predetermined time interval is a two second time interval.

31. The apparatus of claim 18, wherein the ventilation signal is an absolute value of a nasal flow signal.

32. The apparatus of claim 17 wherein the processor is configured to generate the output as a display on an epoch-by-epoch basis that represents the determined absence and/or presence of Cheyne-Stokes breathing.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0024] FIG. 1 is a block diagram of the signal processing pathway leading from the patient's respiratory signal through pre-processing, feature extraction based on epochs, through to classification.

[0025] FIG. 2 shows a typical respiratory signal epoch including a number of “events” (in this case apnea-hyperpnoea sequences). Several primary features are either shown explicitly (2.1 apnea/hypopnoea length, 2.2 hyperpnoea length, 2.3 cycle length) or inferred (2.4 shape features of the hyperpnoea, 2.5 a feature representing the initial “jump” at the beginning of the hyperpnoea).

[0026] FIG. 3 shows details of the calculation of primary features. 3.1 shows the basis functions used in the determination of the hyperpnoea shape features. 3.2 shows an hyperpnoea typical of OSA together with the calculation of the jump feature displayed graphically. 3.3 shows a similar depiction of an hyperpnoea more typical of CS breathing. In both cases the calculated jump features and shape features are tabled.

[0027] FIG. 4 shows examples of epoch classification, e.g., using bar charts, seen in reports 4.1 and 4.2 of the figure.

[0028] FIG. 5 shows the distribution of the normalized max jump in a hyperpnoea signal.

[0029] FIG. 6 shows a cluster analysis of CS and OSA patients.

[0030] FIG. 7 shows results from a LD (linear) classifier.

[0031] FIG. 8 shows results from a QD (quadratic) classifier.

[0032] FIG. 9 shows the correction of data for baseline offset.

[0033] FIG. 10 shows a Cheyne-Stokes flow waveform, the long-term ventilation and the left-shifted 10-second ventilation for a typical patient.

[0034] FIG. 11 depicts those parts of the flow waveform that correspond to hyperpnoeas.

[0035] FIG. 12 shows a Cheyne-Stokes patient's nasal flow signal over about 15 minutes.

[0036] FIG. 13 shows a patient's Spo2 signal (saturation) and ventilation signal (low-pass filtered absolute value of flow).

DETAILED DESCRIPTION

Process Description

[0037] FIG. 1 shows one embodiment of the classification process. While the following may be explained in terms of a sequential process, it is understood that the process can be carried out using a non-linear, non-sequential, or non-staged process, or the order of the process may be changed. Also, while FIG. 1 describes an entire process, aspects of the invention may relate to only a subset of that process. A signal representative of respiration is first recorded from a patient using a logging device which includes a data-acquisition system and a memory. The respiratory signal is then processed either on-board by the recording device or off-line using a computer.

[0038] Preferably, the signal is initially pre-processed. For example, the signal is filtered to remove unwanted noise and, where appropriate, the baseline is zeroed. The signal may also be linearised depending on the transducer used to detect the respiration.

[0039] In the next stage the signal is divided into n epochs of equal length. The epoch length can be as long as the entire record or as short as is practicable to enable detection of respiratory patterns. In one preferred embodiment the epoch length is 30 minutes.

[0040] FIG. 2 shows a typical epoch recording in a patient with Cheyne-Stokes breathing. The shape of the curve is reminiscent of the shape of the Sydney Harbor Bridge and is sometimes referred to by that name. See also FIG. 10. The recording consists of eight “events”, each event consisting of a hypopnoea (in this case also an apnea) followed by a hyperpnoea. For each event an algorithm is used to detect the beginning and end points such that event lengths can be calculated: e.g., apnea/hypopnoea length and hyperpnoea length. A further algorithm may be used to reject events if they do not follow the correct sequence of hypopnoea/apnea-hyperpnoea. Another further algorithm may be used to reject events that fall outside sensible length scale limits.

Determination of Shape Features

[0041] Each hyperpnoea is further processed to derive four so-called “shape features”. These features indicate different shaped hyperpnoeas (bell-shaped versus triangle-shaped for example). The shape features are calculated using singular value decomposition of the hyperpnoea ventilation signal as follows: First, the hyperpnoea is extracted from the respiratory signal and the absolute value is taken of the respiratory signal, giving a ventilation signal. The ventilation signal is scaled by its mean value to give a vector of values V.sub.hyperp. For mathematical convenience the time base of the hyperpnoea [0 . . . T], where T is the end of the hyperpnoea, is mapped to the interval [0 . . . 2π]. A set of four orthogonal functions are calculated and arranged as a 4×m matrix (where m is the number of values in the hyperpnoea signal). A convenient set of orthonormal function are:

[00001]$M_{\mathrm{Basis}}=\left(\begin{array}{c}\frac{1}{\sqrt{\pi }}\sin \left(\frac{t}{2}\right)\\ \frac{1}{\sqrt{\pi }}\cos \left(\frac{t}{2}\right)\\ \frac{\left(3\pi \sin \left(t\right)-8\cos \left(\frac{t}{2}\right)\right)}{\sqrt{\pi \left(9{\pi }^2-64\right)}}\\ \frac{\left(3\pi \cos \left(t\right)+4\sin \left(\frac{t}{2}\right)\right)}{\sqrt{\pi \left(9{\pi }^2-16\right)}}\end{array}\right)$

where t is the time base over the hyperpnoea from 0 to 2π. The basis functions are shown plotted in FIG. 3.1. The four shape features are then calculated as:

F.sub.P(14)=V.sub.hyperp×PseudoInverse(M.sub.Basis),

and are normalized by:

[00002]$F_{P\left(1-4\right)}=\frac{F_{P\left(1-4\right)}}{L_2.Math.F_{P\left(1-4\right)}.Math.}$

where L.sub.2∥ is the L2 or Euclidean norm,

[00003]$\sqrt{\underset{i}{\overset{n}{.Math.}}{x}_i^2}$

[0042] The pseudoinverse M.sup.+ of a matrix M is a generalization of a matrix inverse, and exists for any (m,n) matrix M (for convenience assume m>n). If such a matrix M has full rank (n) one defines: M.sup.+=(M.sup.TM).sup.−1M.sup.T. The solution of Mx=b is then x=M.sup.+b. (Pseudoinverses are useful because of a general theorem stating that F=M.sup.+v is the shortest length least squares solution to the problem MF=v.)

Jump Determination

[0043] Since sudden jumps in the ventilation/flow at the beginning of an hyperpnoea are characteristic of OSA, (see FIG. 2) each hyperpnoea is further processed to derive the so-called “jump” feature, indicative of the extent of any sudden increase in flow, as follows: Again, the hyperpnoea is extracted from the respiratory signal, the absolute value is taken of the respiratory signal, giving a ventilation signal, a droopy peak-detector is used to approximate the envelope of the ventilation signal:

[00004]$e\left[1\right]=v\left[1\right]$$\mathrm{for}i=2\mathrm{.Math.}m$$\mathrm{if}v\left[i\right]\geqq e\left[i-1\right]$$\mathrm{ }e\left[i\right]=v\left[i\right]$$\mathrm{else}$$\mathrm{ }e\left[i\right]=e\left[i-1\right]+\frac{1}{2.5f_s}\left(v\left[i\right]-e\left[i-1\right]\right)\mathrm{end}$

where e[i] is the approximate envelope, f.sub.s is the sampling frequency and v[i] is the ventilation signal. The envelope is interpolated over a new two-second time base (chosen to be roughly the time-length of a breath) to give ei (between non-breathing intervals). The maximum positive difference e.sub.1−e.sub.1(i-1) (over the two second interval) is found in the interpolated signal in the interval between the beginning of the envelope and the point at which the envelope attains its maximum value. Finally, the maximum difference is scaled by the mean value of the ventilation signal to give the “jump feature”. FIGS. 3.2 and 3.3 show this process graphically for two representative hyperpnoeas.

Secondary Feature Determination

[0044] Secondary features are calculated from primary features using the (measure of variation) statistics detailed below. (Note log denotes the logarithm to base e.) First we define the standard deviation as:

[00005]$\mathrm{STD}\left(F_P\right)=\sqrt{\frac{1}{n-1}\underset{i=1}{\overset{n}{.Math.}}{\left(F_{\mathrm{Pi}}-\stackrel{\_}{F_P}\right)}^2},\mathrm{where}\stackrel{\_}{F_P}=\frac{1}{n}\underset{i=1}{\overset{n}{.Math.}}{F}_{\mathrm{Pi}}$

For length measures (e.g. hypopnoea length) and the jump feature the four features are:

[00006]$1.\frac{1}{n}\underset{i=1}{\overset{n}{.Math.}}\log \left(F_{\mathrm{Pi}}\right)$$2.\log \left(\mathrm{STD}\left(\log \left(F_P\right)\right)\right)$$3.\log \left(\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{.Math.}}{\left(F_{\mathrm{Pi}}-F_{P\left(i-1\right)}\right)}^2}\right)\left(\mathrm{To}\mathrm{get}a\mathrm{normed}\mathrm{deviation}\right)$$4.\log \left(\mathrm{STD}\left(F_{\mathrm{Pi}}-F_{P\left(i-1\right)}\right)\right)$

[0045] For hyperpnoea shape features the four features are:

[00007]$1.\frac{1}{n}\underset{1=1}{\overset{n}{.Math.}}\left(F_{Pi}\right)$$2.\log \left(STD\left(\left(F_P\right)\right)\right)$$3.\log \left(\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{.Math.}}{\left(F_{\mathrm{Pi}}-F_{P\left(i-1\right)}\right)}^2}\right)$$4.\log \left(STD\left(F_{\mathrm{Pi}}-F_{P\left(i-1\right)}\right)\right)$

Additional Feature Determination

[0046] Additional features can be calculated using the entire (e.g. 30 minute) epoch signal. One such feature is derived from the spectrogram of the epoch signal and determining that Cheyne-Stokes breathing is present if the spectrogram indicates that the signal has a peak. This feature is calculated as follows: First, the mean of the respiratory signal is calculated and subtracted from the respiratory signal and the resulting signal is chopped into n slices which overlap each other by exactly half the slice length. Each slice is next windowed, preferably using a Hanning window (to reduce edge effects).

[0047] The use of a Hanning window to prepare the data for a FFT is as follows: The FFT function treats the N samples that it receives as though they formed the basic unit of a repetitive waveform: It assumes that if one took more samples they would repeat exactly, with the (N+1) sample being identical to the first sample, and so on. The usual case is that if one's N samples start at one point in a cycle, they end at some other point, so that if one really did play these back in a loop one would get a discontinuity between the last sample and the first. Hence one tapers both ends of the set of samples down to zero, so they always line up perfectly if looped. The formal name for this process is “windowing”, and the “window function” is the shape that we multiply the data by. When the window function is the “raised cosine” 1+cos t the window is termed a Hanning window. Other periodic functions can be used, yielding other windows.

[0048] Next, since CS data appears periodic, a fast Fourier transform is applied to each windowed slice, yielding a complex vector result for each slice. The absolute value is taken of each complex result yielding a real valued vector per slice. The mean is taken of the resulting vectors to yield one vector. The natural log is taken of the subsequent vector and the values in the frequency range 0 Hz to 0.075 Hz are extracted to form a sub-vector, which is then de-trended. Cheyne-Stokes behavior is present if the spectrogram indicates the signal has a peak in the range 0 Hz to 0.075 Hz.

[0049] Briefly, the method of detrended fluctuation analysis is useful in revealing the extent of long-range correlations in time series, where the time series is a vector of data pairs (t.sub.i, x.sub.i), where t represents time and x represents the variable being measured. De-trending consists of subtracting from the x values, values that have been calculated using a polynomial of order n that has been fitted to the data. For example, for order zero the polynomial is simply the mean of all the x values, and that mean is subtracted from the original values. For order one, the polynomial is simply a linear fit to the data (t.sub.i, x.sub.i). Values calculated using the best linear fit are then subtracted from the original values (so removing any linear “trend”). For order two the fitted polynomial is a quadratic, for order three a cubic etc.

[0050] The feature is then calculated as the maximum minus the mean of the de-trended vector. Alternatively one could calculate the entropy of the FFT instead of its peak.

[0051] Additional features can be derived by applying wavelet analysis to each epoch. In this case wavelet coefficients or statistics derived from wavelet coefficients are used as features for the epoch. This yields the location of the peak in time. In wavelet analysis a wave packet, of finite duration and with a specific frequency, is used as a window function for an analysis of variance. This “wavelet” has the advantage of incorporating a wave of a certain period, as well as being finite in extent. A suitable wavelet (called the Morlet wavelet) is a sine wave multiplied by a Gaussian envelope.

Classification

[0052] A subset of features is then selected for use by the classifier. It is known that a particular subset of features can provide more accurate classification than the full set of features. This is caused in part by the so-called “curse of dimensionality”, whereby the required number of training samples increases with the number of features used. The curse of dimensionality causes networks with lots of irrelevant inputs to behave relatively badly: Where the dimension of the input space is high, the network uses almost all its resources to represent irrelevant portions of the space.

[0053] An algorithm is employed to select the best subset based on the training data. Ideally every subset of features should be tested for accuracy and the best subset chosen. The number of subsets is 2.sup.n-1 where n is the number of features. Unless there is a small number of features the exploration of all subsets is impractical and, in any case, accuracy measures tend to be noisy which further hampers the search for the best subset. Alternative algorithms that enable selection of “good” feature subsets include “best first”, “remove worst”, “random start with add and remove”, “simulated annealing” and genetic algorithms.

[0054] A method often used to measure accuracy is 10-fold cross-validation. The training data are split into ten groups or folds and ten different accuracy tests are performed. In each case 9 tenths of the folds are used for training and the resulting classifier is tested for accuracy on the remaining tenth. Statistics are performed on the 10 results to give a measure of accuracy.

Training the Classifier

[0055] Once a feature subset is chosen, the classifier is trained using the entire training data set. A number of classifier types are available including: Baysean maximum likelihood linear and quadratic discriminants, neural networks and support vector machines. In each case a discriminant function is calculated which, when applied to features calculated from new data to be classified, provides probability estimates for different classes. The data (epoch) is assigned to the class with the highest probability.

[0056] In one particular embodiment the discriminant function includes or preferably consists of two weight vectors (of the same length as the feature subset) and two constants. When the desired feature subset has been extracted from the respiratory epoch, the discriminant functions and probability are calculated as follows:

[00008]$d_1=W_1•F+C_1$$d_2=W_2•F+C_2$$\mathrm{probability}=\frac{e^{\left(d1-d2\right)}}{1+e^{\left(d1-d2\right)}}$

[0057] where W.sub.1, W.sub.2 are vectors and C.sub.1, C.sub.2 are constants.

[0058] The probability cutoff may be set at 0.5 in which case a probability of 1.0 would equate to class A and a probability of 0.0 to class B. The cutoff can be adjusted to suit the desired sensitivity and specificity. This is a two-way classification. With suitable training data, a three-way classification is also possible as are even higher n-way classifications.

[0059] In one particular embodiment the classification of each epoch could be displayed in a bar chart as in FIG. 4. Frame 4.1 shows a bar chart where many epochs show a high probability of a class of respiration (in this case CS-like breathing). This provides an “at-a-glance” indication of a patient record. Frame 4.2 shows a bar chart where only a single epoch displays strong CS-like tendency. This provides an indication of where in the patient's record a more detailed investigation is warranted.

Cheyne-Stokes Classifier Based on a Flow Signal or an Spo2 Signal or Both

[0060] The ApneaLink device is capable of measuring an estimate of a patient's flow signal which can be used as an input to the algorithm described herein. Equally there are similar portable devices that can measure and log SpO2, the saturation of oxyhemoglobin in the blood as estimated by pulse oximetry. Pulse oximetry is a simple non-invasive method of monitoring the percentage of haemoglobin (Hb) which is saturated with oxygen. The pulse oximeter consists of a probe attached to the patient's finger or ear lobe which is linked to a computerised unit.

[0061] SpO2 is typically reported as a percentage, values above 95% being normal and those below 95% indicating some degree of hypoxia (lack of oxygen in the blood). Should a patient undergo an apnoea or hypopnoea, it is usual for the SpO2 signal to fall concomitantly with the ventilation, albeit after some delay. During Cheyne-Stokes breathing the SpO2 signal will undergo the classic waxing and waning pattern also characteristic of the ventilation.

[0062] Hence, it is conceivable that the algorithm described herein might use a flow signal estimate (ventilation) or an SpO2 signal or both signals to classify breathing patterns as being typical of Cheyne-Stokes, OSA, mixed apnoeas etc.

[0063] FIG. 12 shows a Cheyne-Stokes patient's nasal flow signal over about 15 minutes. FIG. 13 shows the same patient's Spo2 signal (saturation) and ventilation signal (low-pass filtered absolute value of flow). The signals have been normalized and shifted to display them n the same graph. The same pattern recognition techniques may be applied to both signals. For example: segment the signal into hypopnoeas/hyperpnoeas; analyze the shape of the hypopnoeas; determine the cycle lengths and space ratios; perform a spectrogram (average of absolute value of a number of FFTs); determine peaks in the spectrogram at the CS frequency; determine a morphologic feature in both signals such as the jump feature; perform a continuous wavelet transform on both signals and use ridge finding techniques t follow any CS frequency component over time.

EXAMPLE 1

[0064] A set of data for testing the ability of flow data to be classified into OSA and CS instances consisted of 90 patient studies of approximately 8 hours each. For purposes of the test, both nasal pressure, flow and two effort signals (abdomen & thorax) were recorded, making a confirming diagnosis of the underlying disease possible. The set was divided into 3 groups of 30 patients: OSA (obstructive apnoea), CS and Mixed. The data were further classified (initially) on a 30-minute time-bin basis. The time periods were classified into the following categories: No apnoeas or hypopnoeas (<5) within the time period; Primarily CS breathing (>90%); Primarily OSA (>90%); Primarily (>90%) apnoeas and hypopnoeas of mixed type (i.e. having a central component followed by a number of obstructed breaths); A combination of different events, typically brief periods of CS or mixed apnoeas interspersed between OSA; Patient is moving and the signal is of too low a quality to be useable.

[0065] Typically if CS disease is present, CS breathing will occur in large blocks of at least 20-30 minutes. The data set contained very few periods of “pure” mixed apnoeas. Rather, the mixed group of 30 patients contained periods of OSA, CS breathing or a mixed picture.

Feature Analysis

[0066] All features were analyzed by calculating distributions for the different groups (OSA, Mixed, and CS). The distribution was normalized by application of an appropriate function, for example FIG. 5 shows the distribution of the normalized max jump in hyperpnoea signal between beginning of hyperpnoea and time of peak flow after application of log-to-base-e. The leftmost curves represent a “normal” or Gaussian distribution. It can be seen that the application of the log function has normalized the distributions and, further, that this feature shows good separation between the CS (left) and OSA (rightmost) groups.

Cluster Analysis

[0067] Both k-means and fuzzy k-means clustering techniques were utilized to visualize feature separation power. The features were first averaged on a per-patient basis and then cluster analysis was used to demonstrate a natural clustering into correct groups. FIG. 6 shows such an analysis. The Euclidean 2-norm distance from each of two cluster centers is plotted one against the other. The CS and OSA patients naturally fall into two groups except for one CS patient. The Mixed patients fall into one group or the other depending on the length of time spent during the night in different breathing patterns. The separating diagonal in the figure represents a naive classifier suitable for per-patient grouping. What such a classifier cannot do is find a short period of CS breathing from amongst an otherwise OSA-dominated night.

Feature Temporal Averaging

[0068] The training of a classifier using patterns assigned to individual events is problematical. Temporal averaging was used to reduce the amount of calculation, while also (potentially) increasing statistical power. A 30-minute time-bin was chosen as a best first-guess. After temporal averaging, a new set of per-time-bin patterns is created. The raw features used (visible separation of groups) were: hypopnoea length; hyperpnoea length; 1.sup.st Fourier shape feature; 2.sup.nd Fourier shape feature; and normalized max jump. The time-averaged 30-minute bin features tested were (std=standard deviation, meansq=mean of square of values, sqrt=square root, shift=allows calculation of temporal difference).

Classifier Training and Testing

[0069] Once the data had been processed and the “expert” diagnosis made, a group of 1440 30-minute bins was available for classifier training (90 patients×16 bins).

Classifier Selection

[0070] Numerous statistical methods exist for the training of a classifier from n-dimensional data, e.g.: nearest neighbor, neural nets, cluster analysis. However, because the data “appeared” linearly separable, Bayesian decision theory was used. This theory (which relies on underlying normal probability density functions) uses the minimization of the Bayes error to calculate a discriminant surface. Such a surface separates the data into one of n categories (in this case 2). Both linear and quadratic discriminant functions were utilized. The former separates the data with a hyperplane in m dimensions (where m is the number of features) while the latter separates the data with a hyperquadric. A hyperplane discriminant is always preferred (assuming accuracy of the same order), as it will tend to be well behaved in areas of minimal data coverage.

Over-Optimistic Train and Test

[0071] The classifier was trained using the training set and then the classifier was tested using the same data. This results in over-optimistic values of sensitivity and specificity, as one would intuitively expect. However, again this is an insightful process and one can use a minimal features set (≤3 features) in order to visualize the result. FIG. 7 shows an LD classifier (plane shows equi-probability surface). FIG. 8 shows a QD classifier (quadric shows equi-probability surface).

Results

[0072] During each test the accuracy, sensitivity and specificity were noted as was the current features set (or group of feature sets) with the best accuracy. Estimates of accuracy, sensitivity and specificity resulted of the order of 91%, 91% and 96% respectively.

EXAMPLE 2

Flow Filtering

[0073] The flow is filtered first to remove unwanted and uninteresting high-frequency content. The filter used is a digital FIR (finite impulse response) filter designed using the Fourier method using a rectangular window. The filter has a pass-band from 0 to 0.9 Hz, a transition band from 0.9 to 1.1 Hz and a stop band above 1.1 Hz. The number of terms in the filter varies with sampling frequency. The flow signal is filtered by convolving the time series point-wise with a filter vector.

Ventilation Calculation

[0074] A long-term ventilation signal y.sub.long is calculated using a simple (first order) low-pass filter applied to the flow signal. A time constant of 200 seconds is used (longer than the longest possible cycle of Cheyne-Stokes breathing). In order to measure ventilation (and not mean flow), the filter is applied to the square of the flow signal and the square root is taken of the filter output. Next, a ten-second ventilation y.sub.10 is calculated (a more “recent” measure). This measure is created by convolving the square of the flow signal with a 10-second square wave with an area of one, i.e. a 10-second-long moving average, and then taking the square-root of the result. This filter will have a five second delay constant over the frequency range of interest. For this reason the signal is shifted left by five seconds so that it “lines up” with the original signal for timing purposes. FIG. 10 shows a Cheyne-Stokes flow waveform (large amplitude rapid varying curve), the long-term ventilation (low amplitude slowly varying curve) and the left-shifted 10-second ventilation (moderately varying curve) for a typical patient.

Event Detection from Ventilation Signals

[0075] The 10-second ventilation signal is used to create low and high thresholds for detection of events (hypopnoea-hyperpnoea sequences). The thresholds are:

thresh.sub.low=0.5×y.sub.long;

thresh.sub.high=0.75×y.sub.long;

The timing of events is calculated using the following algorithm:
1. Find all points where y.sub.10<thresh.sub.low.
2. Find all contiguous sections of the above points. These are provisional hypopnoeas.
3. Find all points where y10>thresh.sub.high.
4. Iterate over all of the hypopnoeas identified in step 2. If no points identified in step 3 (hyperpnoeas) fall between hypopnoea n and hypopnoea n+1, then the hypopnoeas n & n+1 are joined together (because no hyperpnoea has been identified between them) to form one hypopnoea. Repeat for all iterations.
5. The hypopnoeas are now confirmed. All inter-hypopnoea regions are considered hyperpnoeas. Each hypopnoea-hyperpnoea “event” constitutes one possible Cheyne-Stokes cycle. E.g. in FIG. 10 there are five cycles shown.

Calculate Event Timings

[0076] Event timings are calculated for each event as follows:

τ.sub.hypopnoea=t(end_of_hypopnoea)−t(beginning_of_hypopnoea);

τ.sub.cycle=t(beginning_of_next_hypopnoea)−t(beginning_of_hypopnoea);

τ.sub.hyperpnoea=t.sub.cycle−t.sub.hypopnoea;

[0077] Obviously the above events will include some unwanted “garbage”. For example, a one-hour-long period of normal breathing bracketed on each side by Cheyne-Stokes breathing will look like a one-hour-long hyperpnoea! (y10 always greater than threshold). Hence, the following sensible limits are applied to the events:

τ.sub.hypopnoea: minimum=10 seconds, maximum=100 seconds,

τ.sub.cycle: minimum=15 seconds, maximum=250 seconds,

τ.sub.hyperpnoea: minimum=5 seconds.

[0078] All events outside these limits are rejected and not processed. We now have event timings and the ability to extract parts of the flow waveform for further analysis. For example, we can iterate over all the events and select out those parts of the flow waveform that correspond to hyperpnoeas. FIG. 11 is an example where we have selected out an hyperpnoea from the above sequence and plotted it separately. In all further processing it is the 1 Hz filtered flow signal that is used for feature extraction.

[0079] While the invention has been described in connection with what are presently considered to be the most practical and preferred embodiments, it is to be understood that the invention is not to be limited to the disclosed embodiments, but on the contrary, is intended to cover various modifications and equivalent arrangements included within the spirit and scope of the invention.