WO2009127799A1

WO2009127799A1 - Method and apparatus for measuring breathing rate

Info

Publication number: WO2009127799A1
Application number: PCT/GB2009/000428
Authority: WO
Inventors: Lionel Tarassenko; Susannah Fleming; David Andrew Clifton
Original assignee: Oxford Biosignals Limited
Priority date: 2008-04-17
Filing date: 2009-02-18
Publication date: 2009-10-22
Also published as: GB0807039D0

Abstract

A method and apparatus for measuring the breathing rate from the electrocardiograph signal (ECG) by measuring the frequency of modulation of the R-R intervals using auto-regressive modelling. The R peaks are detected using a standard QRS detection technique and the R-R intervals are measured to produce a time series of values. These are resampled evenly and the resulting series is windowed in overlapping windows of typically 30 seconds' length, overlapping by 25 seconds, to obtain discrete sections which are AR modelled using an all- pole auto-regressive (AR) model. The AR model allows identification of the dominant frequencies in the signal and the pole corresponding to the breathing rate is identified- by considering its magnitude and the breathing rate it represents. Each 30 second window gives a breathing rate estimate and use of successive windows displaced by 5 seconds results in a breathing rate estimate every 5 seconds. The time series of breathing rate estimates can be Kalman filtered to reject measurements which have a large change in magnitude or represent a large change in breathing rate. The measurements may also be fused with measurements from another sensor.

Description

METHOD AND APPARATUS FOR MEASURING BREATHING RATE

The present invention relates to a method and apparatus for measuring the breathing rate of a human or animal subject.

Measurement of the breathing rate, also known as respiratory rate, is an important aspect of health monitoring, particularly in subjects who are at risk from respiratory suppression. However existing methods of measuring breathing rate have a number of difficulties associated with them. For example, electrical impedance pneumography, which is based on the measurement of the change in electrical impedance across the chest between inspiration and expiration, is prone to artefact, particularly when a patient moves. Airflow sensors, for instance nasal thermistors, which measure the breathing rate by monitoring the difference in temperature between inhaled and exhaled air, are uncomfortable for patients to use over a long period of time. Acoustic monitoring of the throat-has also been proposed, but again this is subject to artefact.

Proposals have been made for deriving a breathing rate signal from the photoplethysmogram (hereinafter referred to as PPG)₅ this being a signal normally used to measure oxygen saturation in the_. blood and the heart rate (also referred to as the pulse rate). The PPG signal is an oscillatory signal at the heart rate, but with a longer period modulation at the breathing rate. Extracting the respiratory waveform from the PPG is not easy and most techniques have not proved entirely satisfactory in practice. Further, the PPG is not always obtained by multi- parameter patient monitoring devices.

Accordingly the present invention provides a method of measuring the breathing rate of a subject, comprising obtaining an electrocardiogram from the subject, measuring the R-peak to R-peak time intervals in the electrocardiogram to obtain a time series of R-R intervals, and further comprising the steps of: modelling each of a plurality of discrete time periods of the time series of R-R intervals using respective autoregressive models, finding poles of the autoregressive models of the discrete time periods of the time series of R-R intervals, and calculating from the poles for each of the discrete time periods the breathing rate for each of the discrete time periods. Thus the present invention is based on the observation that the R-peak to R-peak interval in an electrocardiogram (hereinafter referred to as an ECG) is modulated at the breathing rate. It measures the R-R intervals in the ECG to produce a time series of R-R intervals, and then it uses autoregressive (AR) modelling of successive sections or windows of the time series to measure the frequency content of the modulation of the R-R intervals in each section. Each pole of the AR model represents a particular frequency component, and the breathing rate is obtained by identifying the pole which represents the respiratory waveform. Finding this pole in each successive section of the signal allows calculation of the breathing rate for each of the successive sections.

The invention therefore allows the estimation of breathing rate from the ECG signal. This is highly advantageous as the ECG signal is available from many existing types of patient monitor, including multi-parameter monitors, which do not themselves provide for breathing rate or oxygen saturation measurement. Such an ECG signal is often available too at locations remote from the patient, either relatively local to the patient such as a central nurse station in a ward for a hospitalised patient, or more remotely in the case of telemonitoring, e.g. of patients at home. The invention has the further advantage of allowing subsequent (i.e. non-real time) estimation of breathing rate from patient data sets where the breathing rate was not originally measured.

Preferably the discrete time periods are successive overlapping windows. The amount of overlap and the length of window can be chosen according to the amount of delay acceptable before outputting the breathing rate for a particular time period and depending on the level of accuracy required (longer time periods allow greater accuracy of the AR modelling, but require a longer wait before the result is output). For example the windows may be successive windows from 60 to 300 seconds long, each displaced from the other by 5 to 15 seconds. Ih one example the windows are 120 seconds long each being displaced by 5 seconds. This results in output of a breathing rate every 5 seconds, this corresponding to the average breathing rate over the preceding 120 seconds.

The method may include at least one of the following steps before the step of autoregressive modelling: (a) evenly resampling the time series of R-R intervals

(b) envelope normalisation of the time series of R-R intervals

(c) detrending the time series of R-R intervals before the step of autoregressive modelling (d) high pass filtering the time series of R-R intervals before the step of autoregressive modelling.

Any combination of these four steps may be used. They are effective to emphasize the spectral power of the breathing rate signal and simplify the modelling process.

To identify the pole corresponding to the breathing rate in the AR model, poles are identified which correspond to frequencies in a range of interest, e.g. between 0.066 and 0.36 Hz (4 - 22 breaths per minute), and the pole with the highest magnitude in this range is identified as the dominant pole representing the breathing rate.

The order of the autoregressive model is preferably from 8 to 13, more preferably 12. The higher the order of the model the more accurate the modelling, but the greater the processing time and the greater the number of poles from which the breathing rate pole needs to be identified.

The method above produces a series of breathing rate measurements, one for each window, i.e. one every five seconds in the example above. Preferably this series of measurements is further processed to improve the breathing rate estimate. For example, because it would be unusual for the breathing rate to change greatly from one reading to the next, the series of measurements may be smoothed, e.g. using a Kalman filter. This is effective to reject any outlier measurement outside a predetermined change in breathing rate.

The breathing rate measurements produced by the measurements above may also be combined with breathing rate measurements obtained from another sensor, for example, impedence pneumography, nasal thermistor, PPG etc., and the two measurements may be combined with respective weights derived from the level of confidence in each of the measurements as described in WO 03/051198. Another aspect of the invention provides apparatus for measuring the breathing rate in accordance with the method above. Such apparatus accepts the ECG input and, optionally, a breathing rate input by another sensor, processes the signals as described above and displays the breathing rate on a display.

The invention may be embodied in software and thus extends to a computer program comprising program code means for executing the processing steps in the method, and to a computer-readable storage medium carrying the program.

The invention will be further described by way of example with reference to the accompanying drawings in which:-

Figure 1 schematically illustrates a breathing rate (BR) measurement apparatus according to one embodiment of the present invention;

Figure 2 is a flow diagram of the processing according to an embodiment of the present invention;

Figure 3 (a) illustrates schematically a raw ECG signal;

Figure 3(b) illustrates schematically the R-R interval values plotted against time and with the cubic spline interpolation dashed;

Figure 3(c) illustrates schematically the resampled points on the cubic spline interpolation; and

Figure 3(d) illustrates the poles of the AR model calculated therefrom using a 12 order AR model with a resampled frequency of 2Hz;

Figure 4 illustrates example breathing rate measurements using an embodiment of the invention against a reference breathing rate;

Figure 5 illustrates example breathing rate measurements using an embodiment of the invention against a reference breathing rate; Figure 6(a) and (b) illustrate respectively measurements of the R-R interval for a patient over 120 seconds and the corresponding breathing rate calculated in accordance with an embodiment of the invention;

Figure 7 illustrates the results of measuring the breathing rate in accordance with one embodiment of the invention together with the reference breathing rate and a Kalman filtered version of the estimated breathing rate;

Figure 8 illustrates a time series of R-R intervals with the detected peaks and troughs joined;

Figure 9 illustrates the power spectral density of the time series of Figure 8;

Figure 10 illustrates the time series of R-R intervals of Figure 8 after envelope normalisation; and

Figure 11 illustrates the power spectral density of the envelope normalised time series of Figure 10.

Figure 1 schematically illustrates a breathing rate measurement apparatus in accordance with one embodiment of the invention. The modelling and filtering processes which produce the breathing rate measurement are conducted in a processor 1 and a breathing rate is displayed on a display 2 as a number or a time series or both. The processor 1 receives an ECG input 3 from a standard electrocardiogram and optionally another breathing rate input 4 from a different type of sensor such as a nasal thermistor or an electrical impedance pneumograph.

Before describing the signal processing in this embodiment in detail it may be useful here to give a brief explanation of the general principles of autoregressive (AR) modelling, though AR modelling is well-known, for example, in the field of speech analysis.

AR modelling can be formulated as a linear prediction problem where the current value x(ή) of the signal can be modelled as a linearly weighted sum of the preceding p values. Parameter p is the model order, which is usually much smaller than the length N of the sequence of values forming the signal. Thus:-

P x(ⁿ) = -∑ ^aΛⁿ ~k) + ein) (1)

4-1

The value of the output x(n) is therefore a linear regression on itself, with an error e(ή), which is assumed to be normally distributed with zero mean and a variance of σ². The model can alternatively be visualised in terms of a system with input e(ή), and output x(ή), in which case the transfer function H can be formulated as shown below:

1 z^p

As shown in Equation 2, the denominator of H{z) can be factorised into p terms. Each of these terms defines a root z_;- of the denominator of H(z), corresponding to a pole of H(z). Since H(z) has no finite zeros, the AR model is an all-pole model. The poles occur in complex-conjugate pairs and define spectral peaks in the power spectrum of the signal. They can be visualised in the complex plane as having a magnitude (distance from the origin) and phase angle (angle with the real axis). Higher magnitude poles corresponding to higher magnitude spectral peaks and the resonant frequency of each spectral peak is given by the phase angle of the corresponding pole. The phase angle θ corresponding to a given frequency f, is defined by Equation 3 which shows that it is also dependent on the sampling interval Δt (reciprocal of the sampling frequency):

Thus fitting a suitable order AR model to a signal reveals the spectral composition of the signal. As will be explained below, the pole in an AR model of the ECG signal which corresponds to the breathing rate can be identified from a search of the poles with phase angles within a range defined by the expected breathing frequencies for a normal signal.

To find the poles, the Burg algorithm (J. P. Burg, "Maximum entropy spectral analysis," in Proc. 37th Meeting Soc. Exploration Geophys., 1967) is used. This returns the values of the poles. As each pole z* can be written as a complex number Xk + Wk, the frequency represented by that pole can be calculated from the phase angle of that pole in the upper half of the complex plane:

θ = \2xf^ι ylx = 2πf_k.llf_s (4)

where^ is the sampling frequency,

and the magnitude r is r = (x²+y²)^m.

With that background in mind, Figure 2 illustrates the processing of the ECG signal in accordance with one embodiment of the invention. The ECG signal is received in step 10, this typically being a signal sampled at a high sampling frequency f_s such as 256Hz.

In step 12 the R peaks in the ECG waveform are identified using a QRS detector based on a standard algorithm such as Pan-Tompkins or Hamilton-Tompkins and then in step 14 the intervals from one R peak to the next are measured. Figure 3(a) illustrates schematically an original ECG waveform with the R peaks labelled R and the interval between them Δt. The measured intervals Δt are plotted against time in Figure 3(b) and it can be seen that they are, of course, unevenly spaced (as they are separated on the time axis by Δt which varies). The modulation of the R-R interval Δt with breathing rate is schematically shown in Figure 3(b).

In step 16 the data is windowed into overlapping windows with window lengths of 120 seconds and with the consecutive windows each starting 5 seconds after its predecessor. Longer or shorter window lengths can be used, and different offsets, if desired. Each breathing rate measurement will be, in essence, an average over the window, so although the use of a longer window can increase the precision of the modelling, it also introduces a longer delay and means that the measurement obtained is less representative of the breathing rate at the current instant

In step 18 the unevenly spaced time series of Δt values is re-sampled evenly at a chosen sampling frequency. This can be done by cubic spine interpolation of the measured Δt values (the cubic spine curve is represented dotted in Figure 3(b)), and in this embodiment a re- sampling frequency of 2 Hz is used, as represented by the crosses in Fig. 3(c), though other re-sampling frequencies are possible, for example 4 Hz. The chosen sampling frequency determines the highest frequency that can be measured, as well as affecting the number of data points which have to be processed. A sampling frequency of 2 Hz has been found to give satisfactory results.

In step 20 the signal is de-trended to reduce the influence of any DC offset and to increase the stability of the model.

An AR model is then fitted to the re-sampled de-trended time series in step 22 using a standard technique such as the Burg algorithm and the poles are found The order chosen for the AR model can significantly affect the accuracy of the breathing pole placement. In general it is found that model orders from about 8 to 13 give an acceptable accuracy (error) of about 1 to 2 breaths per minute. In the illustrated embodiment and for the results illustrated later a model order of 12 is used.

In step 26 the dominant pole representing the breathing rate is then found as explained below, and the frequency (breathing rate) that it represents is calculated using Equation 4. To find the breathing pole, poles in the range of 0.066 - 0.36 Hz (4 - 22 breaths per minute) are examined, the pole with the highest magnitude in this range is selected as the breathing pole, and its frequency is used to calculate the breathing rate. If the model does not contain any poles in the range 4-22 breaths per minute, the same process is repeated for the range 0-40 breaths per minute. If this fails to identify any poles the breathing rate is taken to be zero. The calculated breathing rate is then displayed at step 32.

Other ways of identifying the breathing pole are possible: for example by including constraints based on prior knowledge of the expected breathing rate, using a threshold on the magnitude and so on. For example, there is a^'positive correlation between the heart rate and breathing rate which can be used, or past estimates of the breathing rate can be used, possibly over long time windows e.g. up to 5 minutes.

Figure 3(d) illustrates poles obtained using an 11^th order AR model with a downsampling frequency of 1 Hz. Each of the poles is one of a complex conjugate pair and so to identify the pole corresponding to the breathing rate it is necessary only to consider the upper half of the diagram. To identify the pole representing breathing rate both the magnitude of the pole (distance r from the origin) and frequency (phase angle θ) are considered. The pole representing the breathing frequency should have a high magnitude and be in the region of 4.8 - 42 breaths per minute. Therefore in the diagram a sector of interest is defined from 0.08 to 0.7 Hz (4.8 to 42 cycles per minute) and the lowest frequency pole within that sector that has a magnitude of more than 95% of the pole with the highest magnitude within the same sector is identified (labelled BR in the diagram), hi Figure 3(d) the phase angle of the breathing pole labelled BR corresponds to a breathing frequency of 0.27 Hz (16 breaths per minute).

Although the method above gives a good estimate of the breathing rate, which could be directly displayed in step 32, further steps can be taken to reduce the effect of artefacts in the ECG signal. This is particularly important if short length windows are used, hi particular artefacts and other similar signal quality problems can lead into sudden changes in the estimate of breathing rate and reductions in the magnitude of the breathing pole, hi this embodiment of the invention, therefore, a Kalman filter can optionally be used in step 28 to smooth the time series of breathing rate measurements and reject those representing too large a change in breathing rate or pole magnitude.

As is well-known, a Kalman filter uses probabilistic reasoning to estimate the state x of a system based on measurements z. In this case, a 1 -dimensional Kalman filter is used, with x and z each being scalars corresponding to the breathing. The evolution of these quantities is described by Equations 5.1 and 5.2 below.

_X;- = Ax,_{1 + W} w « N(o,Q)

The square matrix A describes how state x evolves over discrete time intervals (i-1 is followed by i), and is known as the "state transition matrix". In this case we do not expect a particular trend in breathing rate over time, so A is the identity matrix. The vector w is the "process noise" describing the noise in the true state, and is assumed to be normally distributed with zero mean and covariance Q. Q can be estimated by observing the statistics of a known breathing rate. z. = Hx, + v

(5.2) v » N(O₅R)

The current value of observation z is assumed to be a function of the current state x with added noise v. The square matrix H describes how the measurement relates to the state, and is known as the "observation matrix". In this case, we measure the state directly, so H is the identity matrix. The "measurement noise", v, is assumed to be normally distributed with zero mean and covariance R. R can be estimated by observing the statistics of measurement errors of a known breathing rate.

Before running the filter, it is necessary to supply initial values for the state estimate X₀ and its covariance matrix P₀. Again, these can be calculated by finding the mean and variance of a true breathing rate. The filter first estimates the new values of x and P using Equations 5.3 and 5.4. x, = Ax_M (5.3) P₍. = AP_wA^r +Q (5.4)

From these, the Kalman gain K and the innovation I can be calculated using Equations 5.5 and 5.6. These will be used to update the prediction of the state x using the current value of the measurement z.

K, = P₁H(HP₄H^2" + R)-¹ (5.5)

I₁ = z, - Hx, (5.6)

If the innovation is too great (e.g. greater than three standard deviations), the measurement is flagged as invalid (rejected), but the update is still carried out.

X₁ = X, + K₁I, (5.7) .

P, = P, -K,HP. (5.8)

The initial values used in this embodiment of the invention are as follows: X₀ = Il P₀ = 20 Q = 0.2 R = 20

Figures 4 and 5 illustrate examples of measurements of the breathing rate of patients using an embodiment of the invention. In Figure 4 120 second windows are used with 4 Hz re- sampling. The AR based measurements are shown in a plot labelled "ar" against a reference labelled "ref . It can be seen that the AR-based measurements follow the reference reasonably well with an average error of 1.34 breaths per minute. Figure 5 illustrates results for a different patient using the same 120 second windowing and 4 Hz re-sampling. In this case the average error is 2.12 breaths per minute, which is higher because of the greater variability in breathing rate.

Figure 6 illustrates the results of applying an embodiment of the invention to a patient being monitored in a high-dependency care unit. Figure 6(a) illustrates the R-R intervals measured from the ECG, and these were windowed into 120 second windows and re-sampled at 2Hz. Figure 6(b) illustrates the breathing rate calculated from these R-R intervals and the reference breathing rate is also plotted. It can be seen that the reference breathing rate is absolutely constant at 20 breaths per minute, suggesting that this patient was being artificially ventilated. The method of the invention gives an average error of 1.4 rate breaths per minute.

Figure 7 illustrates the results of Kahnan filtering of the AR-derived breathing rate estimate. The reference breathing rate is labelled "ref and the basic AR estimated breathing rate labelled "ar" and shown dotted. The Kahnan filtered version of this is labelled "ar-K". This shows that Kahnan filtering does successfully remove any artefactual readings, but that a series of incorrect readings close together cause the signal to be lost as can be seen at 1,250 and 2,900 seconds into the recording. Although unfortunate, this behaviour is both expected and desired as this is also the way that the filter detects and tracks true changes in breathing rate, as seen at 1,000 and 2,000 seconds into the recording.

In a second embodiment further preprocessing steps may be applied after the detrending of step 20 and before fitting of the AR model (step 22). In particular, the ability of the algorithm to select the correct breathing rate can be hampered by the presence of large peaks in the time- series of R-R intervals. For example, the time-series of R-R intervals in Figure 8 has a corresponding power spectral density shown in Figure 9, where it may be seen that a large quantity of power exists at a frequency corresponding to approximately 5 breaths per minute.

I l The actual breathing rate for this time-series of R-R intervals was approximately 25 breaths per minute. Thus, when breathing rate is extracted from this signal using frequency-based techniques, which select the frequency component with dominant power, an incorrectly low breathing rate may be obtained.

In order to obtain correct breathing rates, the time-series of R-R intervals may be processed in step 21b such that its envelope is made to cover the range [-1 I]. This is performed by detecting the peaks and troughs of the time-series of R-R intervals, as shown in Figure 8, where detected peaks and troughs are shown by circles. The envelope of the signal is thus estimated by joining the detected peaks, and joining together the detected troughs (as shown in Figure 8 by lines joining the circles). This envelope is then scaled to the desired [-1 1] range, as shown in Figure 10.

The power spectral density of the resulting envelope-normalised signal is shown in Figure 11. It may be seen that the power associated with the dominant large peaks in the original time- series of R-R intervals has been reduced, allowing the power associated with the actual breathing rate to be clearly visible, at 25 breaths per minute in the figure. Frequency-based breathing rate estimation will thus correctly select the breathing rate, because it corresponds to the frequency component with the largest amount of power.

It is also possible, as illustrated by step 21a, to high-pass filter the time series of R-R intervals, e.g. at 0.1 Hz (which corresponds to a breathing rate of 6 breaths per minute), before performing the envelope normalisation.

In some circumstances, despite such steps, the accuracy of the AR breathing pole may be low. In order to allow estimation of the breathing rate even in such circumstances the AR breathing rate estimate may be fused in step 30 with an estimate from another breathing rate sensor, for example a conventional sensor. The two measurements may be fused using the technique described in WO 03/051198 which uses a one-dimensional Kalman filter on each series of breathing rate measurements to calculate a confidence value associated with that measurement, this confidence value then being used as a weight in combining the two breathing rate estimates.

Claims

1. A method of measuring the breathing rate of a subject, comprising obtaining an electrocardiogram from the subject, measuring the R-peak to R-peak time intervals in the electrocardiogram to obtain a time series of R-R intervals, and further comprising the steps of: modelling each of a plurality of discrete time periods of the time series of R-R intervals using respective autoregressive models, finding poles of the autoregressive models of the discrete time periods of the time series of R-R intervals, and calculating from the poles for each of the discrete time periods the breathing rate for each of the discrete time periods.

2. A method according to claim 1 wherein the time periods are successive overlapping windows.

3. A method according to claim 2 wherein the start points of the successive overlapping windows are displaced from each other by 5 to 15 seconds and the windows are from 60 to 300 seconds long.

4. A method according to claim 3 wherein the start points of the successive overlapping windows are displaced from each other by 5 seconds and the windows are 120 seconds long.

5. A method according to any one of the preceding claims further comprisingat least one of the following steps before the step of autoregressive modelling: (a) evenly resampling the time series of R-R intervals (b) envelope normalisation of the time series of R-R intervals;

(c) detrending the time series of R-R intervals before the step of autoregressive modelling

(d) high pass filtering the time series of R-R intervals before the step of autoregressive modelling.

6. A method according to any one of the preceding claims further comprising the step of identifying as a pole representing the breathing rate a pole which has the greatest magnitude within a predetermined frequency range.

7. A method according to claim 6 wherein the predetermined frequency range corresponds to a predefined allowable range of breathing rates.

8. A method according to any one of the preceding claims wherein the order of the model

9. A method according to any one of the preceding claims wherein the order of the model is 12.

10. A method according to any one of the preceding claims further comprising the step of filtering the breathing rate measurements for the plurality of discrete time periods to produce a smoothed breathing rate measurement.

11. A method according to claim 10 wherein said filtering step comprises Kalman filtering.

12. A method according to claim 10 or 11 further comprising rejecting measurements which represent greater than a predetermined change in breathing rate.

13. A method according to any one of the preceding claims further comprising the step of combining the measurement of breathing rate with a breathing rate measurement from another sensor.

14. A method according to claim 13 wherein the combination is a weighted combination with weights derived from confidence in the measurements.

15. Apparatus for measuring the breathing rate of a subject, comprising an input for receiving an electrocardiogram from the subject, a data processor for executing the steps of measuring the R-peak to R-peak time intervals in the electrocardiogram to obtain a time series of R-R intervals, modelling each of a plurality of discrete time periods of the time series of R- R intervals using respective autoregressive models, finding the poles of the autoregressive models of the discrete time periods of the time series of R-R intervals, and calculating from the poles for each period the breathing rate for each time period as defined in any one of claims 1 to 15.