SIGNAL PROCESSING METHOD AND APPARATUS

Abstract

A method and apparatus for estimating the frequency of a dominant periodic component in an input signal by modelling the input signal using auto-regressive models of several different orders to generate candidate frequencies for the periodic component, generating synthetic sinusoidal signals of each of the candidate frequencies, and calculating the cross-correlation of the synthetic signals with the original signal. The frequency of whichever of the synthetic signals has the highest cross-correlation with the original signal is taken as the estimate of the frequency for the dominant periodic component of the input signal. The method may be applied to any noisy signal which has a suspected periodic component, for example physiological signals such as photoplethysmogram signals, and in the estimation of heart rate and breathing rate from such physiological signals.

Claims

1. A method of analysing an input signal to estimate the frequency of a periodic signal therein comprising the steps of: receiving the input signal; generating frequency domain candidates for the periodic signal; for each frequency domain candidate generating a synthetic signal of the same frequency; calculating in the time domain a measure of the similarity of each synthetic signal with the input signal; outputting as the estimate of the frequency of the periodic signal the frequency of the synthetic signal with the maximum similarity with the input signal.

2. A method according to claim 1 wherein the measure of similarity is the correlation or the cross-correlation.

3. A method according to claim 1, wherein the frequency domain candidates are generated by performing spectral analysis on the input signal to detect its spectral components.

4. A method according to claim 3 wherein the spectral analysis is by autoregressive modelling.

5. A method according to claim 4 wherein a plurality of autoregressive models of different orders are fitted to the input signal.

6. A method according to claim 5 wherein the plurality of autoregressive models comprise models of order 2 to 20.

7. A method according to claim 4, wherein synthetic signals are generated for each pole in the fitted autoregressive models.

8. A method according to claim 4, wherein synthetic signals are generated for only poles in the fitted autoregressive models which correspond to a predetermined allowed frequency range for the periodic signal.

9. A method according to claim 4, wherein synthetic signals are generated for only poles in the fitted autoregressive models which have a predetermined magnitude.

10. A method according to claim 4, wherein synthetic signals are generated for only the dominant poles in the fitted autoregressive models.

11. A method according to claim 1 wherein the frequency domain candidates are generated by selecting a plurality of different frequencies within a predetermined frequency range.

12. A method according to claim 11 wherein the plurality of different frequencies are evenly spaced from each other over the frequency range.

13. A method according to claim 1 wherein the input signal is temporally windowed, for example from 5 to 60 seconds long, and the steps of calculating the similarity and outputting a frequency estimate are performed for each window.

14. A method according to claim 13 wherein the windows are overlapping windows, for example with an overlap from 0.5 to 10 seconds.

15. A method according to claim 1, wherein the input signal is obtained by detrending and filtering a noisy signal.

16. A method according to claim 15 wherein the input signal is obtained by resampling the detrended and filtered noisy signal.

17. A method according to claim 1, wherein the synthetic signal is sinusoidal.

18. A method according to claim 1, further comprising the step of defining a similarity threshold and inhibiting the outputting step if the similarity of the synthetic signal with the maximum similarity with the input signal is below the threshold.

19. A method according to claim 1, wherein the input signal is a time series of image intensity data.

20. A method according to claim 1 wherein the input signal is video image data comprising a time series of intensity data for each colour channel.

21. A method according to claim 1 wherein the input signal is a measurement of a physiological parameter of a human or animal subject.

22. A method according to claim 21 wherein the periodic signal is the heart rate or respiration rate of the subject.

23. A computer program comprising program code means for executing on a computer system the method of claim 1.

24. Apparatus for analysing an input signal to estimate the frequency of a periodic signal therein comprising: an input for receiving the input signal; a processor for processing input signal; the processor being configured to execute the steps of claim 1; the apparatus further comprising an output to output as the frequency estimate the frequency of the synthetic signal with the maximum cross-correlation with the input signal.

Description

[0022] The invention will be further described by way of example with reference to the accompany drawings in which:—

[0023] FIG. 1 is a flow diagram illustrating the method of one embodiment of the invention;

[0024] FIG. 2 shows the results of applying the embodiment of FIG. 1 to respiratory rate estimation and comparing it with other techniques for a patient with spontaneous breathing;

[0025] FIG. 3 shows the results of applying the embodiment of FIG. 1 to respiratory rate estimation and comparing it with other techniques for a patient with controlled breathing; and

[0026] FIG. 4 schematically illustrates signal processing apparatus according to an embodiment of the invention.

[0027] FIG. 1 is a flow diagram showing the application of an embodiment of the invention to the analysis of a PPG signal. In step 20 a PPG signal, which is the pulsatile waveform produced by a pulse oximeter at one of its two wavelengths (red and infra-red) is received. In step 22 this signal is detrended (e.g. by subtracting the mean or a linear fit over the whole signal) and filtered using a low-pass filter (Kaiser window, f.sub.pass=0.5*cardiac frequency in Hz, f.sub.stop=1.2*cardiac frequency). The signal is resampled, for example at a frequency of 5 Hz. This resampling is required in this embodiment because the intention is to find the respiration rate. The PPG signal itself has a high sampling rate of 100-250 Hz to ensure that the shape and heart rate information are preserved. However at such high sample rates the phase angles corresponding to breathing frequencies are very small and the cardiac-synchronous pulsatile component of the PPG signal (heart beat) is dominant. Downsampling ensures that the cardiac-synchronous pulsatile component is no longer dominant.

[0028] In step 24 the signal is divided into consecutive 32-second windows, each separated by 3 seconds. Longer or shorter windows can be used, and rather than being separated, the windows can overlap.

[0029] In step 26 multiple auto-regressive models with model orders ranging from 2-20 are applied to the windowed signal. In step 28 for each AR model the dominant pole, that it is to say the pole with the highest magnitude is identified. Thus 19 dominant poles, one from each of the model orders, is produced. Although in this embodiment only the dominant pole is taken from each model, it is possible alternatively to take all poles within a certain allowed frequency range, or to take all poles of each model onto the next step. Clearly taking more poles means a greater amount of processing.

[0030] In step 30, for each of the selected poles (in this example the nineteen dominant poles, one from each model order) a synthetic signal which is a pure sinusoidal signal at the frequency represented by the pole's phase angle and whose amplitude is determined by that of the signal that was modelled (for example by setting it equal to the standard deviation of the original signal) is generated. This sinusoid is therefore a time-domain representation of the dominant pole from each model. Then in step 32 the cross-correlation is calculated between each of the synthetic signals and the original modelled signal for all possible phase differences between the synthetic signal and the modelled signal. Each calculation will return a cross-correlation coefficient c. In step 34 whichever of the synthetic signals has the highest cross-correlation with the original signal is found (the maximum value of c) and in step 36 the frequency of this synthetic signal is output as the estimate of the strongest periodic component of the original input signal. In this example, analysing a downsampled PPG signal, the estimated respiration rate in breaths per minute can be obtained by taking the reciprocal of this frequency in Hertz and multiplying by 60.

[0031] FIG. 2 illustrates the results of applying this example of the invention to a PPG recording in comparison with capnometry data as a reference and the auto-regressive modelling technique described in WO-A2-2013/027027. The capnometry reference measurements are shown in dotted, the results of applying the embodiment of the invention are shown in a thin line with ticks, and the results of the prior art auto-regressive modelling technique are shown in a thick line with dots. FIG. 2 is for a patient with spontaneous breathing whereas FIG. 3 shows results for a patient with controlled breathing (in FIG. 3 the reference is shown dashed, the prior art AR model in a solid line, and the embodiment of the invention in a dotted line with crosses). As can be seen the method of this embodiment of the invention tracks the breathing rate much more accurately than the prior art auto-regressive modelling method.

[0032] A specific example of the invention for application to analysis of a PPG signal has been described above. The invention is application, however, to any noisy time series of data where one is seeking a single periodic signal in the time series.

[0033] Thus more generally a given input time-series I(x) may is modelled by:

I(x)=μ(I(x))+Σ.sub.ja.sub.j cos(ω.sub.j+φ.sub.j)+ε(x) [0034] Where μ(I(x)) is the average offset or trend, ε(x) is a random error associated with each data-point, and Σ.sub.j(a.sub.j cos(ω.sub.jx+φ.sub.j)) defines the Fourier series that describes the signal.

[0035] Assuming that this signal has a suspected dominant frequency Ω in a particular frequency range, it is detrended to obtain:

I′(x)=detrend{I(x)}=Σa.sub.j cos(ω.sub.jx+φ.sub.j)+ε(x) [0036] and finally filtered so as to be modelled by:

O(x)=filter{I′(x)}=A cos(Ωx+Φ)+ε(x) [0037] For any given set of frequencies ω.sub.k={ω.sub.1, . . . , ω.sub.n}={2πf.sub.1, . . . , 2πf.sub.n} that are thought to best approximate Ω, a synthetic sinusoidal signal is generated for each frequency S.sub.k(x), such that the generated signal is

S.sub.k(x)=cos(ω.sub.kx+λ) [0038] This signal is cross-correlated with the original detrended and filtered signal

[00001]$O\left(x\right).Math..Math.\mathrm{over}.Math..Math.\mathrm{the}.Math..Math.\mathrm{lag}-\frac{1}{2.Math..Math.f_k}<\lambda \le \frac{1}{2.Math..Math.f_k}$$R_k\left(\lambda \right)=\left(O\left(x\right)*S_k\left(x\right)\right)\left[\lambda \right]$ [0039] The frequency ω.sub.k found to maximise the absolute value of R.sub.k is then selected as the best approximation of the dominant signal frequency Ω, i.e. find ω.sub.k that satisfies

argmax{(O*S.sub.k)[λ]},∀k

[0040] Standard methods of cross-correlation calculation may be used to evaluate R.sub.k. Usually one may express R.sub.k as:

R.sub.k(λ)=Σ.sub.x(O[x]−Ō).Math.(S[x−λ]−S) [0041] Alternatively, a more computational efficient model is to use the convolution theorem, i.e. that:

F{f*g}=(F{f})*.Math.F{g} [0042] where F denotes the Fourier Transform and * denotes the complex conjugate operand. [0043] We may then write the cross correlation of the original signal with the synthetic signal as:

R=(λ)=(O(x)*S.sub.k(x))[λ]=F.sup.1{(F{O})*.Math.F{S}} [0044] Finally, the relative amplitudes of the two signals need to be accounted for. This can be done in a number of ways:
i) Take this into account when generating the synthetic signal, i.e. by setting

S.sub.k(x)=B(x)cos(ω.sub.kx+λ)

where B(x) matches the signal envelope. In its simplest form,

B(x)=std(O(x))

ii) Normalise when cross-correlating, i.e. through:

[00002]$R_k\left(\lambda \right)=\frac{{ℱ}^{-1}.Math.\left\{{\left(ℱ.Math.\left\{0\right\}\right)}^{*}.Math.ℱ.Math.\left\{S\right\}\right\}}{\sqrt{\left(0*0\right)\left[0\right].Math.\left(S*S\right)\left[0\right]}}$

or simply through:

[00003]$R_k\left(\lambda \right)=\frac{{ℱ}^{-1}.Math.\left\{{\left(ℱ.Math.\left\{0\right\}\right)}^{*}.Math.ℱ.Math.\left\{S\right\}\right\}}{\mathrm{std}\left(0.Math.\left(x\right)\right).Math.\mathrm{std}\left(S\left(x\right)\right)}$

[0045] In the above examples the candidates for the periodic component of the input signal are selected using auto-regressive modelling. However it is possible to generate these candidates in other ways. For example, instead of auto-regressive modelling another spectral analysis technique can be used to identify candidate spectral components. For example peaks in the Fourier transform spectrum or a wavelet transform method may be used. Alternatively, if the allowed frequency range for the signal of interest is relatively small, it is possible to generate “all” frequencies in that range (with a certain resolution separating them) and to generate synthetic signals corresponding to each of those frequencies.

[0046] Also in the above example cross-correlation is used as the time domain validation technique but other validation techniques in the time domain can be used. For example calculating the correlation or any other method that maximises the linearity between the synthetic signal and the input signal.

[0047] The invention may be embodied as a computer program for running on a computer system which receives the input time series and outputs an estimate of the dominant periodic component within the time series. Alternatively the invention may be incorporated into a dedicated apparatus, such as a monitor, as schematically illustrated in FIG. 4. Such apparatus comprises an input 40 for receiving the input data, a process 42 for executing the data processing steps described above and an output 44, such as a display, for outputting the estimate of the dominant periodic component, for example as a physiological measurement such as a breathing rate or heart rate.