WO2017029506A1

WO2017029506A1 - Method and apparatus for determining a repeat period in periodic data

Info

Publication number: WO2017029506A1
Application number: PCT/GB2016/052557
Authority: WO
Inventors: Philip ASTON; Mark Christie; Gary CHAFFEY; Manasi NANDI
Original assignee: University Of Surrey; King's College London
Priority date: 2015-08-18
Filing date: 2016-08-18
Publication date: 2017-02-23
Also published as: GB201514637D0

Abstract

A method and apparatus for determining a repeat period in periodic data are disclosed. The repeat period is determined by obtaining a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window, and obtaining a plurality of transformed data points by transforming each of the vectors of delay coordinates into a coordinate system comprising one or more predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points. The repeat period can then be determined based on the transformed data points. In some embodiments an optimised value of a time delay is determined, for which the spread of the transformed data points around a known axis is minimised, and the repeat period is then determined from the optimised value of the time delay and a known relationship between the optimised time delay and the repeat period. In other embodiments, the repeat period is determined as the time taken to complete one loop in an attractor defined by the transformed data. In some embodiments a measure of cardiac contractility is determined based on a peak in the density of a predefined vector as a measure of time, where the predefined vector is defined such that a transformed data point for an exactly periodic signal lies on the predefined vector when a certain relationship between the time delay and the repeat period is satisfied. Furthermore, in some embodiments two or more repeat periods within a signal may be determined, if the signal under analysis contains more than one frequency component.

Description

Method and Apparatus for Determining a Repeat Period in Periodic Data

Technical Field

The present invention relates to a method and apparatus for determining a repeat period in periodic data. More particularly, but not exclusively, the present invention relates to determining a repeat period in pseudo-periodic data, for example

determining a heart rate from blood pressure data. Background of the Invention

There are many applications where it is necessary to determine the repeat period of periodic data. For example, in a medical context, a measurement of the repeat period in an approximately periodic physiological signal such as blood pressure can provide valuable diagnostic information. Many methods have been proposed for determining cycle lengths from approximately periodic data such as a blood pressure trace or an Electrocardiography (ECG) signal. The simplest approach is to determine a point of interest on each cycle such as a maximum or minimum point. However, a drawback of this approach is that irregularities within the data can result in extra maxima or minima in each cycle, or can result in a point of interest such as the peak of the R wave for an ECG being missing in some cycles. Accordingly, the repeat period may be incorrectly calculated.

The invention is made in this context. Summary of the Invention

According to a first aspect of the present invention, there is provided apparatus for determining a repeat period in periodic data, the apparatus comprising: a delay coordinate unit configured to obtain a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window; a coordinate transforming unit configured to obtain a plurality of transformed data points by transforming each of the vectors into a coordinate system comprising a plurality of predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points; and a repeat period calculating unit configured to determine the repeat period based on the transformed data points. In some embodiments according to the first aspect, the delay coordinate unit is configured to obtain each vector of delay coordinates based on a time delay τ, the predefined vectors are selected such that a transformed data point for an exactly periodic signal lies on one of the predefined vectors( v ) and the magnitude of another one of the predefined vectors ( w ) is zero when a relationship between τ and the repeat period is satisfied, and the apparatus further comprises: a time delay optimising unit configured to determine an optimised value of τ for which the spread of the

transformed data points around said one of the predefined vectors ( ) is minimised, by finding the value of τ which minimises a norm of the values of said other one of the predefined vectors ( w ) for the plurality of transformed data points, wherein the repeat period calculating unit is configured to determine an average repeat period in the approximately periodic data based on the optimised value of τ. For example, the periodic data may be blood pressure data, and the apparatus can be further configured to determine a density in said one of the predefined vectors ( v ) as a measure of time, over a plurality of repeat periods in the periodic data, and determine a measure of cardiac contractility based on the area under a peak in the density, or based on the magnitude of the peak. Furthermore, in some embodiments according to the first aspect a measure of the overall periodicity maybe derived, for example the standard deviation of the periods, in order to give an indication as to how regular the signal is. Such information may have diagnostic value, for example in enabling early detection of failure of a component in an engineering context.

In some embodiments according to the first aspect, the delay coordinate unit is configured to obtain each vector of delay coordinates based on a first time delay τΐ, and is further configured to obtain a second delay coordinate from one of said one or more predefined vectors based on a second time delay τ2, the time delay optimising unit is further configured to determine an optimised value of the second time delay τ2 which minimises a function related to the difference between said one of the predefined vectors and the second delay coordinate, and the repeat period calculating unit is further configured to determine an average repeat period in the approximately periodic data from the optimised value of the second time delay τ2. The function related to the difference between said one of the predefined vectors and the second delay coordinate may be a scale invariant function. In some embodiments according to the first aspect, the plurality of transformed data points define a projection of an attractor of the periodic data along an axis normal to the plurality of predefined vectors, and the repeat period calculating unit is configured to determine the repeat period by determining a time taken to complete one loop of the obtained attractor in the projection along said axis. The repeat period calculating unit may further be configured to determine a first time point as being a time at which a trace of the periodic data crosses a section arranged to cut one side of the attractor in the projection along said axis, determine a second time point as being a time at which the trace of the periodic data next crosses the section, and determine the time taken to complete one loop as the difference in time between the first time point and the second time point. The repeat period calculating unit may be configured to determine a time at which the trace of the periodic data crosses the section by interpolating between two consecutive data points falling on opposite sides of the section, and/or maybe configured to arrange the section to be perpendicular to said one side of the attractor. The delay coordinate unit can be configured to obtain three-dimensional vectors of delay coordinates [x, y, z], and the axis normal to the predefined vectors is x = y = z.

In some embodiments according to the first aspect, the delay coordinate unit is configured to obtain first and second vectors of delay coordinates for each one of the plurality of samples of the periodic data in the time window, the first vector of delay coordinates being obtained using a first time delay and the second vector of delay coordinates being obtained using a second time delay longer than the first time delay, the apparatus further comprising: a time delay optimising unit configured to obtain an optimised value of the first time delay while holding a value of the second time delay constant, and subsequently obtain an optimised value of the second time delay using the optimised value of the first time delay, wherein the repeat period calculating unit is configured to calculate a first repeat period from the optimised value of the first time delay and to calculate a second repeat period from the optimised value of the second time delay. In such embodiments, the time delay optimising unit may be configured to obtain the optimised value of the second time delay by applying a time delay

transformation to the plurality of samples of the periodic data to obtain a plurality of transformed signal samples, each transformed signal sample being obtained from a difference in values between one of the plurality of samples of the periodic data at a first point in time and one of the plurality of samples of the periodic data at a second point in time, the time delay optimising unit being further configured to set the difference between the first and second points in time based on the optimised value of the first time delay, wherein the delay coordinate unit can be configured to recalculate the first and second vectors of delay coordinates for the transformed signal samples, wherein the coordinate transforming unit can be configured to transform each of the recalculated vectors into the coordinate system comprising the plurality of predefined vectors, and wherein the time delay optimising unit can be configured to obtain the optimised value of the second time delay from said plurality of predefined vectors.

According to a second aspect of the present invention, there is provided a method of determining a repeat period in periodic data, the method comprising: obtaining a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window; obtaining a plurality of transformed data points by transforming each of the vectors of delay coordinates into a coordinate system comprising one or more predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points; and determining the repeat period based on the transformed data points. It should be understood that 'periodic' does not mean that a signal represented by the data is perfectly periodic. Periodic data may also be referred to as 'pseudo-periodic' or 'approximately periodic' data.

In an embodiment according to the second aspect, each vector of delay coordinates is obtained based on a time delay τ, the one or more predefined vectors are arranged such that a transformed data point lies on a predefined axis when a relationship between τ and the repeat period is satisfied, and the repeat period is determined by: determining an optimised value of τ for which the spread of the transformed data points around the predefined axis is minimised; and determining the repeat period from the optimised value of τ. The predefined axis may be one of the predefined vectors, such that the magnitude of another one of the predefined vectors is zero when the relationship between τ and the repeat period is satisfied, and the optimised value of τ may be determined by finding the value of τ which minimises the norm of the values of said other one of the predefined vectors for the plurality of transformed data points.

The periodic data may be physiological data, for example blood pressure data, and the method can further comprise: determining a density in said one of the predefined vectors as a measure of time, over a plurality of repeat periods in the periodic data; and determining a measure of cardiac contractility based on the area under a peak in the density, or based on the magnitude of the peak. In another embodiment according to the second aspect, each vector of delay

coordinates is obtained based on a first time delay τΐ, and the method further comprises: obtaining a second delay coordinate from one of said one or more predefined vectors based on a second time delay τ2; determining an optimised value of the second time delay τ2 which minimises the difference between said one of the predefined vectors and the second delay coordinate; and determining an average repeat period in the approximately periodic data from the optimised value of the second time delay τ2. In another embodiment according to the second aspect, the plurality of transformed data points define a projection of an attractor of the periodic data along an axis orthogonal to the one or more predefined vectors, where movement along said axis corresponds to baseline variation in the periodic data, and the repeat period is determined by determining a time taken to complete one loop of the obtained attractor in the projection along said axis. Determining the time taken to complete one loop of the obtained attractor may comprise: determining a first time point as being a time at which a trace of the periodic data crosses a section arranged to cut one side of the attractor in the projection along said axis; determining a second time point as being a time at which the trace of the periodic data next crosses the section; and determining the time taken to complete one loop as the difference in time between the first time point and the second time point. A time at which the trace of the periodic data crosses the section can be determined by interpolating between two consecutive data points falling on opposite sides of the section. In some embodiments, the section is arranged to be perpendicular to said one side of the attractor.

The obtained vectors of delay coordinates may be three-dimensional vectors [x, y, z], and the axis normal to the predefined vectors may be x = y = z.

In some embodiments according to the second aspect, obtaining the vector of delay coordinates comprises obtaining first and second vectors of delay coordinates for each one of the plurality of samples of the periodic data in the time window, the first vector of delay coordinates being obtained using a first time delay and the second vector of delay coordinates being obtained using a second time delay longer than the first time delay, and the method further comprises: obtaining an optimised value of the first time delay while holding a value of the second time delay constant; and obtaining an optimised value of the second time delay using the optimised value of the first time delay, wherein a first repeat period is determined from the optimised value of the first time delay and a second repeat period is determined from the optimised value of the second time delay. In some embodiments according to the second aspect, the optimised value of the second time delay is obtained by: applying a time delay transformation to the plurality of samples of the periodic data to obtain a plurality of transformed signal samples, each transformed signal sample being obtained from a difference in values between one of the plurality of samples of the periodic data at a first point in time and one of the plurality of samples of the periodic data at a second point in time, a difference between the first and second points in time being set based on the optimised value of the first time delay; recalculating the first and second vectors of delay coordinates for the transformed signal samples; transforming each of the recalculated vectors into the coordinate system comprising the plurality of predefined vectors; and obtaining the optimised value of the second time delay from said plurality of predefined vectors.

The periodic data may be physiological data.

According to a third aspect of the present invention, there is provided a non-transitory computer-readable storage medium arranged to store software instructions which, when executed on one or more processors, perform a method according to the second aspect.

According to a fourth aspect of the present invention, there is provided apparatus for determining a repeat period in periodic data, the apparatus comprising: processing means for executing computer program instructions; and computer-readable memory arranged to store computer program instructions executable by the processing means, wherein the stored computer program instructions include instructions for performing a method according to the second aspect.

According to a fifth aspect of the present invention, there is provided apparatus for determining a repeat period in periodic data, the apparatus comprising: processing means for executing computer program instructions; and computer-readable memory arranged to store computer program instructions executable by the processing means, wherein the stored computer program instructions include instructions for: obtaining a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window; obtaining a plurality of transformed data points by transforming each of the vectors of delay coordinates into a coordinate system comprising one or more predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points; and determining the repeat period based on the transformed data points.

In an embodiment according to the fifth aspect, the computer program instructions may be further configured to perform any of the methods according to the second aspect.

Brief Description of the Drawings

Embodiments of the present invention will now be described, by way of example only, with reference to the accompanying drawings, in which:

Figure ι illustrates blood pressure traces from a healthy mouse;

Figure 2 is a flowchart showing a method of determining a repeat period in periodic data, according to an embodiment of the present invention;

Figure 3 illustrates a trajectory of a blood pressure signal in three-dimensional reconstructed phase space derived from the data in the lower panel of Fig. 1, according to an embodiment of the present invention;

Figure 4 is a flowchart showing a method of determining a repeat period in periodic data by determining a time taken to complete one loop in a projection of an attractor, according to an embodiment of the present invention;

Figure 5 illustrates a projection of the trajectory of Fig. 3 when viewed along the axis x = y = z, according to an embodiment of the present invention;

Figure 6 is a flowchart showing a method of determining an average repeat period in periodic data by obtaining an optimised value of the time delay, according to an embodiment of the present invention;

Figure 7 illustrates transformed data points used when determining the average repeat period in the method of Fig. 6, according to an embodiment of the present invention; Figure 8 schematically illustrates apparatus for determining a repeat period in periodic data, according to an embodiment of the present invention; and

Figure 9 illustrates an idealised, piecewise linear signal, according to an embodiment of the present invention;

Figure 10 is a graph plotting the value of v as a function of time (t) during one repeat period, for an idealised signal; Figure 11 is a graph plotting the magnitude of the delta functions at the end points of a density;

Figure 12 illustrates the x points at which v₂ = o and dv₂/ dt < o in the blood pressure data from the top panel in Fig. 1;

Figure 13 is a graph plotting the individual cycle lengths in the data from the top panel in Fig. 1, as determined by the method shown in Fig. 4;

Figure 14 illustrates the x points at which the trajectory in the (v,w) plane intersects a Poincare section similar to the one shown in Fig. 5, in the blood pressure data from the top panel in Fig. 1;

Figure 15 is a graph plotting the individual cycle lengths in the data from the top panel in Fig. 1, as determined by the method shown in Fig. 4;

Figure 16 is a flowchart showing a method of determining first and second repeat periods from a compound signal, according to an embodiment of the present invention; and

Figure 17 is a flowchart showing a method of determining an optimised value of the second time delay in the method of Fig. 16, according to an embodiment of the present invention.

Detailed Description

Embodiments of the present invention can be used to analyse periodic data and determine a repeat period within the periodic data. Depending on the embodiment, the repeat period that is determined may be an individual repeat period for a single cycle within the periodic data, may be an average value of the repeat period over a plurality of cycles, or may quantify the variability of the value of the repeat period. For example, in some embodiments the method can be used to determine the Heart Rate Variability (HRV) from heart rate data. The periodic data may relate to an approximately periodic signal, that is, a signal which comprises repeating similar waveforms but which is not perfectly periodic. As such, it should be understood that 'periodic' in this context does not mean that a signal represented by the data is perfectly periodic. The observed periodicity can vary slightly from one repeat period to the next, and may vary significantly over longer time periods. This is particularly true of physiological signals, but also applies to any type of pseudo-periodic data in other contexts, for example in engineering applications. Similarly, any references herein to a 'repeating waveform' in the periodic data should be understood as meaning a waveform that is similar, but not necessarily identical, from one repeat period to the next. An example of data which can be analysed in an embodiment of the invention is shown in Fig. l, which illustrates blood pressure (BP) traces for a healthy mouse. The top panel of Fig. ι illustrates a ι second sample of BP data from a healthy, conscious mouse. As shown in the top panel of Fig. l, the BP data constitutes a pseudo-periodic signal comprising a waveform that is almost identical, to a close approximation, from one repeat period to the next. The blood pressure data shows cyclic behaviour due to the regular heartbeat. The beat-to-beat interval for this data is approximately 120 milliseconds (ms), corresponding to a heart rate of 500 beats per minute, and so there are around 120 data points per cycle. However, the blood pressure signal is not strictly periodic in a mathematical sense, as beat-to-beat variability in both period and waveform shape results in variation in this signal over time.

As shown in the bottom panel, which illustrates a 10 second sample of BP data, baseline variation over longer time periods can complicate the determination of the repeat period when a simple thresholding technique is used, since the maximum and minimum values of the signal within each repeat period can vary significantly over time.

The blood pressure data illustrated in Fig. 1 was collected using an implanted radiotelemetry device at a sample rate of 1000 Hertz (Hz) for a 24 hour period, resulting in a series of large datasets. In the present embodiment, blood pressure data is recorded as a digital signal comprising a plurality of samples. Each sample represents the magnitude of the quantity being measured, in this case blood pressure, at a point in time. In other embodiments, samples can be obtained from a received analogue signal. The invention is not limited to a sample rate of 1000 Hz, and in other embodiments any suitable sample rate can be used, depending on the typical timescale for the signal type being measured.

A method of analysing periodic data to determine a repeat period will now be described with reference to Figs. 2 and 3, in the context of the BP data shown in Fig. 1. Figure 2 is a flowchart showing a method of determining a repeat period in periodic data, according to an embodiment of the present invention, and Figure 3 illustrates a trajectory of a blood pressure signal in three-dimensional reconstructed phase space. The trajectory in n-dimensional reconstructed phase space can also be referred to as the 'attractor' of the signal, and is obtained by plotting a vector of delay coordinates for each data point. To determine the repeat period in periodic data, such as that shown in Fig. 1, a vector of delay coordinates is first obtained for each one of a plurality of samples of the periodic data in a time window, in step S201. The delay coordinates can be obtained based on Takens' theorem, which allows an attractor to be reconstructed in an n-dimensional 'phase space' for n≥2 from a single signal x(t) by using a vector of delay coordinates as follows:

[x{t), x{t - T), x{t - 2r),..., x{t - (n - l)r)] where x(t) is the signal value at time t, and τ > o is a fixed time delay.

When using Takens' delay coordinates to reconstruct an attractor, it is necessary to choose appropriate values for two key parameters, namely (i) the embedding dimension n, and (ii) the value of the time delay τ.

With regard to the embedding dimension n, Takens showed that an m-dimensional manifold requires 2m+i delay coordinates for there to be a continuous, injective mapping from the manifold to the reconstructed phase space. However, the dimension of the underlying model or of the attractor of that model is generally not known, particularly in the case of complex physiological signals. Taking the present example of a blood pressure signal, there are many factors that influence the blood pressure in a conscious animal, including the sympathetic nervous system, parasympathetic nervous system, respiratory system and motor activity. Accordingly, a full mathematical model including all possible factors would be highly complex. Even if such a model could be constructed, the dimension of the attractor would still not be known. Thus, a theoretical value for the embedding dimension n is not available. However, there are a variety of methods in the literature for estimating a minimum embedding dimension, for example using a singular value analysis or the method of false nearest neighbours.

In the present embodiment, the value n = 3 is chosen for convenience, without attempting to find the optimal embedding dimension. However, in other embodiments different values can be chosen for the embedding dimension n. For example, n = 2 could be used in some embodiments, or a value higher than 3 could be chosen. With regard to the time delay τ, if τ is chosen to be very small then there will be only a small difference between the variables. As a result, the trajectory will always lie close to the axis in the phase space given by i = x₂ = ... = x_n, where Xi(t) = x(t - (i - ι)τ) for i = 1, n. On the other hand, if τ is chosen to be very large, then there may be little correlation between each of the variables.

Investigations by the inventors have shown that when τ is equal to 1/3 or 2/3 of the period Tof a perfectly periodic signal, the closed loop attractor has Z₃ cyclic symmetry when viewed along the axis x = y = z. For physiological data, the value of τ is chosen to make the density of the projected attractor as close to being Z₃ symmetric as possible by minimising a symmetry measure. It can be shown that the optimal time delay is then related to the average cycle length of the signal in the time window. In other embodiments, a fixed value of τ could be used or a different measure optimised. In the present embodiment, the trajectory shown in Fig. 3 is plotted using a time delay τ of 31 milliseconds (ms), for the 10-second trace shown in the lower panel of Fig. 1. Examples of methods that can be used for choosing the time delay τ in other embodiments of the present invention include a mutual information based method, and a continuity statistic based method. The continuity statistic based method can be used to determine both the optimal time delay and embedding dimension simultaneously.

As described above, the embedding dimension n = 3 is used in the present

embodiment. Accordingly, two new variables y(t) and z(t) are defined, allowing a three-dimensional vector of delay coordinates to be obtained for each sample of the signal x(t) within a time window as follows:

where

As described above, the trajectory shown in Fig. 3 is obtained by plotting the obtained vectors for data within the defined time window, which in the present example is 10 seconds. In other embodiments, data can be plotted for a time period of longer or shorter duration than 10 seconds.

After obtaining the vectors of delay coordinates for a plurality of samples, in step S202 a plurality of transformed data points are obtained by transforming each of the vectors of delay coordinates into a coordinate system comprising one or more predefined vectors. Specifically, the predefined vectors are arranged so as to factor out baseline variation in the periodic data from the transformed data points. For example, when three-dimensional vectors of delay coordinates are used, baseline variation in the signal will result in movement along the axis x = y = z. In general, for an n-dimensional vector comprising n variables x x₂, ... x_n, where n≥2, the axis ¾ = x₂ = ... = x_n can be referred to as the central axis. The predefined vectors can be chosen to be orthogonal to the central axis in order to factor out the baseline variation. As a result of the coordinate transformation in step S202, the repeat period can be determined in step S203 without being affected by baseline variation in the signal. Examples of methods for performing the coordinate transformation in step S202 and determining the repeat period in step S203 will now be described, with reference to Figs. 4 to 7.

Figure 4 is a flowchart showing a method of determining a repeat period in periodic data by determining a time taken to complete one loop in a projection of an attractor, according to an embodiment of the present invention. First, in step S401 a vector of delay coordinates is obtained for each sample of the periodic data within the time window, as described above. The vectors of delay coordinates define an attractor of the periodic signal in n-dimensional reconstructed phase space.

Then, in step S402 a coordinate transformation is performed to obtain a projection of the attractor onto a plane orthogonal to the central axis, to factor out baseline variation in the signal. A projection of the attractor from Fig. 3 is illustrated in Fig. 5. As shown in Fig. 5, the trace 501 of the attractor on the plane orthogonal to the central axis is approximately triangular in the present embodiment.

An explanation of the coordinate transformation method used in step S402 of Fig. 4 is as follows. For a reconstructed attractor in a three-dimensional space, such as the one shown in Fig. 3, if the time delay is chosen to be the trivial case of τ = o then x = y = z, and so the trajectory in the reconstructed phase space simply moves up and down the line in the phase space defined by x = y = z, which is the central axis. For small τ, the trajectory will stay near to this axis, since the three points will be close together, and so the dominant motion will again be up and down this central axis. However, for larger τ, the trajectory will start to move away from this axis. As the blood pressure rises and falls, the motion projected onto the central axis will also rise and fall. Therefore a new coordinate system is defined, consisting of the central axis together with two axes that are orthogonal to this axis, and each other. We then normalise each of these vectors to give an orthonormal basis.

Although in the present embodiment the new coordinate system comprises three orthogonal unit vectors, in other embodiments the new vectors may not be unit vectors, in which case the step of normalising the vectors to give an orthonormal basis can be omitted. Also, in some embodiments the vectors may not be orthogonal, in which case an alternative central axis can be defined other than x = y = z.

In more detail, a unit vector in the direction of the central axis x = y = z is given by

V! = (1, 1, i)^T/V3. In order to project the attractor onto the plane normal to the central axis, the remaining two basis vectors must be orthogonal to this one and to each other. In the present embodiment we choose two other predefined unit vectors as v₂ = (1, i,-2)^T/V6 and v₃ = (1,-1, o)^T/V2, however, any rotation of these two vectors in the plane orthogonal to the central axis is also possible. Thus, the matrix M, which has columns v_1; v₂ and v₃, is an orthogonal matrix. If we have coordinates (u, v, w) with respect to the new basis vectors, then the old and new coordinates are related by:

or equivalently x = u where x = (x, y, z)^rand u = (u, v, w)^T. Thus, the new coordinates are defined by u = ^rx since Mis an orthogonal matrix, or equivalently

1

v {x + y - 2z)

From this, we see that u is almost the mean of the three original variables x, y and z, and in some embodiments u can be redefined to be the mean.

The new variables v and w are the coordinates of a point in the three-dimensional phase space projected onto the plane orthogonal to the line x = y = z. The trajectory shown in Fig. 3 projected onto this plane is shown in Fig. 5. The new variable u(f) captures the vertical motion (baseline variation) of the blood pressure signal, whilst the other two variables v(t) and w(f) are not affected by this motion.

Next, in steps S403 to S405 the repeat period is determined as being the time taken to complete one loop of the obtained attractor in the projection onto the (y, w) plane. In the present embodiment, the time taken to complete one loop of the attractor is determined by defining a Poincare section 502 arranged to cut one side of the attractor in the (y, w) plane, and to be substantially perpendicular to the side of the attractor, as shown in Fig. 5. The time taken to complete one loop is then given by the difference in time between a first time point at which a trace of the periodic data crosses the section, and a second time point at which the trace of the periodic data next crosses the section.

In scenarios where one side of the attractor is approximately horizontal and the attractor describes loops around the origin in the (y, w) plane, such as in the case of an attractor for blood pressure data as shown in Fig. 5, the Poincare section 502 can be defined as the line v = o with w < o. In the original (x, y, z) coordinates, these conditions are equivalent to: x + y - 2z - 0 , with x < y

The first condition can also be expressed as:

z = ^ {x + y)

That is to say, z is the average of the two other coordinates. There are generally two points in each cycle when this occurs, but the one of interest (i.e. where the trace intersects the Poincare section) is defined by the additional inequality that x < y . By recording consecutive time points at which the trajectory intersects the Poincare section, the cycle length for a single repeat period can then be found as the difference between consecutive time points.

Alternative choices for the Poincare section in the present embodiment could be lines that bisect the other two sides of the triangular attractor, given by a rotation of the line v = o with x < y by 2π/3 or 471 /3. A rotation by 2π/3 in the (y, w) plane corresponds to a cyclic permutation of the original three variables. Thus the alternative Poincare sections are defined by: x - ^ (y + z), with y < z

y = ^ (x + z), with z < x

When a data point does not coincide with the section, the time at which the trace of the periodic data crosses the section can be determined by interpolating between two consecutive data points falling on opposite sides of the section. By making the Poincare section 502 perpendicular to the side of the attractor which is cut by the section, the repeat period can still be accurately determined even when the course of the trace deviates substantially from the ideal attractor during a single repeat period, causing the trace to cut the Poincare section 502 at different points in the (u, w) plane. If the Poincare section 502 was sharply inclined to the edge of the attractor, such deviation could result in the cycle length being determined incorrectly. Although in the present embodiment the embedding dimension 71=3 is used, in other embodiments a higher embedding dimension may be used (n≥4). In such

embodiments, it is still possible to restrict to well-defined two dimensional subspaces, defined by eigenvectors associated with a complex conjugate pair of eigenvalues that are roots of unity, on which there is approximate cyclic symmetry. A similar approach can be used in this case also by defining a half line that will generally, but not necessarily, go through the origin that cuts through the attractor in the given plane. The time points can again be found as successive intersections of the trajectory with this line, from which a repeat period can be determined.

Furthermore, in some embodiments an embedding dimension of n=2 may be used. In such embodiments, new coordinates can be defined as:

* - y) where v₂ is the coordinate that factors out baseline variation in the signal. This variable will simply move up and down the line in the interval -h≤v₂≤ h, where h is the pulse pressure. Time points could again be identified by the condition that v₂ = o with dv₂/dt > o (or < o), with the difference between consecutive such points defining the individual cycle lengths. Figures 12 and 13 illustrate an example of results obtained by applying this n=2 method to the mouse blood pressure data from the top panel of Fig. 1. Figure 12 illustrates the x points at which v₂ = o and dv₂/ dt < o, and Fig. 13 illustrates the individual cycle lengths in the data as determined by this method.

Although in the present embodiment the repeat period is determined by detecting points of intersection with a Poincare section, in other embodiments other methods can be used to determine the period. For example, the time intervals between the maximum or minimum points of w{t) or v{t) could be used to determine the period, although this approach could suffer from the same drawbacks as finding the maximum or minimum points of the original signal as described earlier, and may therefore be less robust than the Poincare based approach described above. In the method shown in Fig. 4, an initial value of the time delay τ to be used when calculating the vectors of delay coordinates can be selected in various ways. For maximum accuracy, the value of the time delay τ should be chosen so as to maximise the symmetry of the attractor in the projection onto the plane orthogonal to the central axis. For example, in one embodiment a starting value of τ is chosen based on a symmetry measure which provides an indication of how symmetric the attractor is. Various definitions of the symmetry measure are possible, depending on the desired shape of the attractor. In the present embodiment, it is desired to find a value of the time delay τ which gives an attractor that is as close to being Z₃ symmetric as possible when projected onto the (v, w) plane. For an ideal periodic signal, this attractor will have this threefold rotational symmetry around the x = y = z axis, and so a matrix ¾(τ) can be defined as:

D_S {T) = ^ {D₁ {T) + D₂ {T) + D₃ {T))

where the first density matrix Α(τ) = D (τ), the second density matrix Α>(τ) is a density matrix obtained for the data rotated in the (v, w) plane by 2π/3, and the third density matrix Ι¾(τ) is a density matrix obtained for the data rotated in the (v, w) plane by 4^/3·

A symmetry measure Six) is then defined as:

In the ideal case of a perfectly repeating waveform with a constant repeat period T, and a time delay τ equal to (i/3)Tor (2/3)Γ, ¾(τ) will be equal to D (τ) and will be equal to zero. The symmetry measure can therefore be used to select a suitable value of the time delay τ, by computing S( ) for different values of τ to find a value of the time delay τ which minimises the function S( ).

However, in other embodiments this procedure can be omitted, and different methods can be used to select the time delay τ. For example, the mean peak-to-peak distance in the physiological data can be used to provide an initial estimate of the repeat period T, or the repeat period Tcan be found using autocorrelation, and an initial value of the time delay τ can be chosen as 1/3 or 2/3 of this repeat period T. Alternatively, the time delay τ could be fixed or could simply be chosen by trial-and-error.

Once an initial value of the time delay τ has been selected, which may be defined based on an average or approximate value of the repeat period T, the method shown in Fig. 4 can then be used to accurately determine the length of each individual cycle within the time window. Furthermore, in some embodiments the value of the time delay τ may be periodically updated during ongoing analysis of the periodic data, to account for variation in the repeat period Tover time.

Although in the present embodiment the attractor is constructed in three dimensions, and then projected onto a plane orthogonal to a chosen axis to give a two-dimensional projection, the invention is not limited to these numbers of dimensions. In other embodiments, any number of dimensions greater than or equal to two can be used for the delay coordinate vectors. For example, to extend the delay coordinate vectors described above to n dimensions, new variables ; can be defined as follows: x_i = x(t - (i - ί)τ) where i = 1, .... , n. The central axis can then be defined as the axis

^X\ ~ ^X2 ~ ^Xn and additional vectors v„ orthogonal to the central axis can be defined as necessary. For example, additional mutually orthogonal vectors can be defined as follows:

v

₂ = (i -i,o,...,o)^r

v

₃ = (l,l,-2,0,...,0)^r „ = (l,...,l,-(n - l))^r

Coordinates of the vectors v₂,...v_n give a projection along the central axis in this case. An alternative method of performing a coordinate transformation and determining the average repeat period will now be described, with reference to Figs. 6 and 7. Figure 6 is a flowchart showing a method of determining a repeat period in periodic data by obtaining an optimised value of the time delay, and Fig. 7 illustrates transformed data points configured to factor out baseline variation in the periodic data, according to an embodiment of the present invention. In the method of Fig. 6, the repeat period is determined based on a known relationship between the repeat period and an optimised value of the time delay used to obtain the vector coordinates. First, in step S601 a vector of delay coordinates is obtained for each sample of the periodic data within the time window, as described above. The vectors of delay coordinates define an attractor of the periodic signal in n-dimensional reconstructed phase space. Then, in step S602 a coordinate transformation is performed. In the present embodiment, new coordinates are defined such that a transformed data point lies on a known axis when the time delay τ is a predefined fraction of the repeat period T.

Specifically, in the present embodiment new coordinates v and w are defined such that when τ = Γ/2, the transformed data points for a perfectly periodic signal lie entirely on the v axis.

Next, in step S603 an optimised value of the time delay τ is determined for which the spread of the transformed data points around the v axis is minimised. Then, in step S604 the repeat period is determined from the optimised value of the time delay τ, based on the known relationship between the time delay τ and the repeat period T. Specifically, in the present embodiment, once the optimised value of τ has been determined then the repeat period T = 2τ.

A method of defining the new coordinates v and w is as follows. In the present embodiment, we consider the case of the variables v and w described above with reference to Figs. 4 and 5, when the time delay τ = Γ/2. In the case of an exactly periodic signal, taking τ = Γ/2 gives z(t) = x(t - 2τ) = x(t - T) = x(t). Replacing z with x in the definition of the variable v gives:

or equivalently, w -— J3v . Thus, when τ = Γ/2 the orbit of the attractor collapses onto a line through the origin with slope - 3 . For convenience, in the present embodiment the new variables v and w are defined based on an anticlockwise rotation of the (v, w) plane, such that this line corresponds to the horizontal axis:

sinf- cos 3 J

Relating these back to the original delay coordinates x, y, z, gives:

w

When τ = Γ/2, z = x and w - 0 , meaning that the periodic orbit of the attractor lies entirely along the v axis. For approximately periodic data, in which the repeat period may vary from one cycle to the next, the transformed (v , vv) data points are spread out in the w direction as shown in Fig. 7, and will be most concentrated around the v axis when the time delay τ is equal to half of the average cycle length T. This value of the time delay τ can be referred to as the optimised time delay, and can be found by minimising a suitable measure of the spread in the w direction.

In the present embodiment, the optimised value of the time delay τ is determined by finding the value of τ for which the taxicab norm of the w values, , is minimised.

An alternative would be to minimise , or the Euclidean norm , or any other norm of w . Another alternative would be to minimise w / p: where x(t) is the signal.

The value of ||wfl/| | is scale invariant, meaning that the value will remain constant when the signal is scaled by some constant a. Minimising the value of a scale invariant function such as ||w||/||xj| ensures that the calculation of the optimised time delay will not be affected due to scaling of the signal, for example if the signal is attenuated due to losses. Yet another alternative would be to minimise the value of a density function, since a density function is also scale invariant. For example, a first density Α(τ) can be derived from the (v , w) data, and a second density Α(τ) from the reflected data

(v ,-vv) . The density will be as symmetric as possible relative to a reflection in the v axis if we minimise the function:

I (T) - ( (r ) + D₂ = \ D_x (r ) - D₂ (r|

In comparison to methods based on density functions, minimising the taxicab or Euclidean norm offers the advantages that there are no densities or matrix norms to compute, and that more weighting is given to points further away from the v axis, meaning that the minimisation process tries to concentrate the mass as close to the v axis as possible. By using a method such as the one shown in Fig. 6, an average value of the repeat period for a plurality of data samples within a time window can be calculated in a computationally efficient manner.

The method described above minimises the norm of the difference between x(t) and x(t-2x). If there is significant baseline variation in the signal, then this can result in a difference between x(t) and x(t-2x) for some values of t even if τ is correctly chosen as τ=Γ/2. Thus, the method can also be applied to a signal in which the baseline variation has been removed, such as v. In this case, 1 1 v(t)-v (t-20) 1 1 can be minimised to find Q=T/ 2. The new parameter Θ could be chosen to be the same as τ or it could be different. Any other norm could also be used, and the variable w could also be used instead of v, or other similarly defined variables. This variation of the method is not affected by baseline variation, allowing the repeat period to be reliably calculated even when there is substantial variation in the signal level over time, for example as shown in the bottom trace of Fig. 1.

In the method described above with reference to Figs. 6 and 7, the repeat period is obtained by finding an optimised value of the time delay τ which was used when deriving the original vector of delay coordinates [x,y,z]. In a related embodiment, a second time delay τ₂ can be used to obtain a second delay coordinate from one of the variables obtained following the coordinate transformation. For example, for the variable v we can define a new delay coordinate v₂ as follows: v₂ t) = v t - T₂ )

The variable v constitutes a periodic signal with the same repeat period as the original input signal. The average repeat period can then be determined by finding the optimised value of the second time delay τ₂ which minimises the difference between v and v₂ . For a plurality of samples of v and v₂ , the second time delay τ₂ can be optimised by minimising a norm of the quantity ( v - v₂ ).

Referring now to Fig. 8, apparatus for determining a repeat period in periodic data is schematically illustrated, according to an embodiment of the present invention. As shown in Fig. 8, the apparatus 800 comprises a delay coordinate unit 801, a coordinate transforming unit 802, and a repeat period calculating unit 803. The delay coordinate unit 801 is configured to obtain a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window, and may comprise any suitable means for obtaining the vector of delay coordinates. The coordinate transforming unit 802 is configured to obtain a plurality of transformed data points by transforming each of the vectors received from the delay coordinate unit 801 into a coordinate system

comprising a plurality of predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points, and may comprise any suitable means for obtaining the plurality of transformed data points. The repeat period calculating unit 803 is configured to determine the repeat period based on transformed data points received from the coordinate transforming unit 802, and may comprise any suitable means for determining the repeat period. In the present embodiment, the apparatus 800 further comprises an output unit 805 configured to output the determined repeat period. The output unit 805 may, for example, be configured to display information about the repeat period on a display unit, or may comprise a network interface configured to transmit information about the repeat period to another device. In the present embodiment, the apparatus 800 is configured to implement a method such as the one shown in Fig. 6, and accordingly further comprises a time delay optimising unit 804 configured to determine an optimised value of τ, as described above with reference to Fig. 6. The repeat period calculating unit 803 can then calculate the repeat period from the optimised value of τ, based on the known relationship between τ and the repeat period. In other embodiments, when a different method of determining the repeat period is used, the time delay optimising unit 804 may be omitted.

The various elements of the apparatus 800 can be implemented using hardware, software, or a combination of both. In software-implemented embodiments, one or more elements can be implemented as software instructions in one or more computer programs stored on a non-transitory computer-readable storage medium. When executed on one or more processors, the software instructions cause the relevant steps of the methods disclosed herein to be performed. In the present embodiment, the delay coordinate unit 801, coordinate transforming unit 802, repeat period calculating unit 803, and time delay optimising unit 804 are implemented in software, and the apparatus 800 further comprises computer-readable memory 806 configured to store computer program instructions which, when executed, perform the functions associated with each of these elements. In embodiments where the necessary functionality is implemented in dedicated hardware, the computer-readable memory 806 may be omitted.

In some embodiments, a method such as the one shown in Fig. 6 can be used to obtain information relating to cardiac contractility. Contractility is related to the inherent strength and vigour of the heart's contraction during systole, and a measurement of contractility can be of diagnostic value. Changes in contractility result in movement of the breakpoint within a single repeat period of a blood pressure signal. Specifically, if each waveform in a blood pressure signal is considered a piecewise linear signal consisting of two linear segments (systole and diastole, respectively) with the break between them occurring at approximately 1/3 of the period, it can be shown that when the break point occurs earlier in the cycle, there is a characteristic increase in the density at the corners of a projected triangular attractor such as the one shown in Fig. 5. Conversely, if the break point occurs later in the cycle, this will result in a reduction in the density and rounding of the corners of the attractor derived from blood pressure data. These changes are characteristic of increases and reductions in cardiac contractility, respectively and may be used for the diagnosis of cardiac changes induced by sepsis, shock, cardiotoxic drugs, or pharmacological agents with a therapeutic action on the heart.

In principle therefore a measure of the density at the corners of a projected triangular attractor for a blood pressure signal will contain useful information about the contractility. However, there are difficulties associated with determining contractility directly from a projected triangular attractor. If the attractor consists of a thick band then the magnitude of the peaks in the density at the corners will be reduced, in comparison to a case in which the attractor is a thin band. Therefore the magnitude of the peak can be influenced by the thickness of the bands along each side of a projected triangular attractor, which may obscure the effect of changes in contractility. Also, if the breakpoint of the periodic piecewise linear idealised signal occurs at aX, where Tis the period of the signal, then there is only a peak in the density for o≤ a≤ 1/3. If a > 1/3 then it becomes difficult to distinguish different cases.

In some embodiments of the present invention, a method based on the one shown in Fig. 6 can be adapted to overcome these problems and provide a more reliable indication of cardiac contractility. The method can be understood by considering the theoretical example of an idealised, piecewise linear signal as shown in Fig. 9. If x(t) is periodic with period T and τ = Γ/2, x(t) = At + B for some range of t and y(t) also lies on this line, then it can be shown that v {t) - -^/ (v - x) = -AT ΐ β and so is constant. In the example idealised signal shown in Fig. 9, the breakpoint between the two lines occurs at t = aXfor some <≡ . This function is defined by:

Assuming that the optimised value of the time delay, τ, which has been obtained using the method shown in Fig. 6 is equal to T/ 2, then for (a+V2 T≤ t≤ T, x and y are both on the downstroke of the signal which is a straight line with slope -h/((i- a)T) and so v = v₀ , where:

Since the pulse pressure h and the breakpoint a are constant during one repeat period, v is constant during the time interval (a+V2)T≤ t≤ T. Similarly, for oT≤ t≤ T/2, x and y are again both on the downstroke but their order is reversed due to the periodicity of the signal. Therefore, in this case v = -v₀ and again is constant during this time interval.

During the time interval Γ/2≤ t≤ (a+V2)T, x is on the downstroke and y is on the upstroke. During this interval, v is given by:

which is a linear function of time. It can be shown that:

This clearly constitutes motion at constant speed from left to right between the two end points ± v₀. Similarly, when o≤ t≤ aT, there is motion at constant speed between the two end points from right to left.

Thus, the overall motion, starting from t = o, consists of uniform motion from v₀ to - v₀ , followed by no movement for a time of ( /2-a)T, followed by uniform motion back from - v₀ to v₀ , finally followed by no movement for another time interval of length ( /2-a)T. Figure 10 is a graph plotting the value of v as a function of time (t) during one repeat period, for an idealised signal. It follows from Fig. 10 that the density in v₀ as a measure of time is constant between - v₀ and v₀ with value oT/ v₀ (twice the reciprocal of the velocity as the interval is traversed twice) together with two delta function spikes at the two end points, each of magnitude ( /2-a)T.

In the present embodiment, a measure of contractility is obtained as follows. After optimising the time delay τ using a method such as the one shown in Fig. 6, such that the magnitude of the vector w is zero for a perfectly periodic signal when τ=Τ/ 2, the density in v as a measure of time is determined over a plurality of repeat periods in the periodic data. The integral of the density is then the total time for one cycle, namely T, and in the present embodiment the total density is normalised to one by dividing by T. The normalised density therefore has two delta functions at the end points of magnitude Vi-a. The magnitude of the delta function in the normalised density, as a function of a, is plotted in Fig. 11. The value of a, which represents a measure of contractility, can therefore be obtained by determining the magnitude of the peak in the density, that is, the magnitude of the delta functions at the end-points. In other embodiments, the step of normalisation can be omitted, in which case the magnitude of the delta functions will be given by (V2-a)T, where T is the average repeat period over the time interval being studied.

From this idealised example, several observations can be made:

• The spikes in the density vary for <≡ , which is better than the method based on a projected triangular attractor which only gives peaks in the density for or e (θ, {].

• The one-dimensional density in v is symmetric about the origin and is constant between the end points.

• The magnitude of the delta function spikes at the end points depends only on a.

For smaller values of a, the spikes are larger and vice versa. The spikes vanish when a = Vi.

• The constant value of the density in the middle region depends on the

breakpoint a and the pulse pressure h. In terms of a, the maximum value is obtained when a = Vi. • v₀ , which defines the two end points of the density interval, also depends on a and h. When a = o, v₀ = h/ [6 and when a = Vi, v₀ = 2h/ 6 . Clearly these two values differ by a factor of 2.

• The clearest measure of contractility is therefore the magnitude of the spikes at the end points of the interval, which depends only on a and not on h.

When Vi < a < 1, a similar analysis shows that the density is very similar to the previous case and consists of two delta functions located at v = ±v_l where v_l = h/ j6 ) of magnitude a - V2 and a constant density in between these end points with value (l - a)v_l - 6a(l - a)/ h . The results for a≤ Vi and a > Vi can be combined. Thus, in general, the density has two delta functions located at v = ±v_d where

- (A/6 (ja - -_j| + -_j)) of magnitude |γ - cr| and a constant density in between the end points with value 6a(l - a)/ h . The magnitude of the delta functions is shown in Fig. 11 as a function of a.

An approximately periodic signal that is similar in nature to the idealised signal shown in Fig. 9 should therefore generate a density in v that is qualitatively similar to the one described above. Thus, we would expect two peaks at either end with a flatter region in the middle which is symmetric about zero. The best measure of contractility (a) is either the height of the peaks or the area under the peaks, since the other features also depend on the pulse pressure (h). By determining the contractility based on the magnitude of peaks or the area under the peaks in the density in one of the predefined vectors, in this case v , embodiments of the present invention can therefore provide a reliable measure of contractility that is independent of the pulse pressure or other factors.

Examples of results obtained from using the above-described methods to analyse the mouse blood pressure data shown in the top panel of Fig. 1 are illustrated in Figs. 12 to 15. Figures 12 and 13 relate to the n=2 method as described above.

Figures 14 and 15 relate to the method described above with reference to Fig. 4. Figure 14 illustrates the x points at which the trajectory in the (y,w) plane intersects a Poincare section similar to the one shown in Fig. 5, and Fig. 15 illustrates the individual cycle lengths in the data as determined by this method.

Methods as described above in Figs. 4 and 6 have also been applied by the inventors to analyse the noisy data shown in the lower panel in Fig. 1, and the results are shown below in Table 1. The first row gives the mean, standard deviation, maximum, and minimum across all 109 cycles, measured in ms, as determined by the method of Fig. 4 with embedding dimension n = 2. The second row gives the mean, standard deviation, maximum, and minimum across all 109 cycles, measured in ms, as determined by the method of Fig. 4 with embedding dimension n = 3. Both methods correctly identify all 109 cycles in the data.

Table 1 From a comparison of the individual cycle lengths plotted in Figs. 13 and 15, and the data shown in Table 1, it can be seen that the results obtained using both methods are in close agreement, indicating that both methods are able to reliably detect and measure the length of individual cycles in pseudo-periodic data. Embodiments of the invention have been described in which blood pressure data, which is a type of pseudo-periodic signal, is analysed. However, in other embodiments of the invention different types of physiological signals can be analysed, including but not limited to arterial blood pressure, pulse oximetry, central venous pressure,

Electrocardiography (ECG), intracranial pressure, respiration data, Left Ventricular Pressure (LVP), Electroencephalography (EEG), Electromyography (EMG), PPG,

Doppler blood flow, accelerometry (tremor), gait monitoring, continuous blood glucose, and acoustic (speech) signals. Furthermore, in other embodiments the data being analysed may not be physiological. For example, in some embodiments data obtained in an engineering context may be analysed, such as vibration data from machinery.

In embodiments of the invention, a vector of delay coordinates is obtained for each one of a plurality of samples of periodic data. Depending on the embodiment, the periodic data may be the original signal itself, for example heart-rate data, or may be a signal derived from the original. For example, in some embodiments the above-described methods can be applied to a periodic data such as υ, ιυ, ν or w , which are themselves periodic signals with the same repeat period as the original data.

Furthermore, other embodiments of the invention may be used to determine two or more repeat periods within a signal, if the signal under analysis contains more than one frequency component. A signal which contains multiple frequencies may be referred to as a compound signal. Examples of types of compound signals include but are not limited to Amplitude Modulation (AM), Baseline Wander (BW), and Frequency

Modulation (FM) signals. In some embodiments of the invention, a compound signal maybe analysed using a delay coordinate method to determine the repeat period, and hence the frequency, of each signal component within the compound signal. Referring now to Fig. 16, a flowchart is illustrated showing a method of determining first and second repeat periods from a compound signal, according to an embodiment of the present invention. In the present embodiment, a compound signal which includes two component signals of different frequencies is analysed using a first vector of delay coordinates (¾, y and a second vector of delay coordinates (x₂, y₂, ¾). First, in step S1601 each vector of delay coordinates is calculated using a different time delay τ, and τ₂, as follows: x_x (t) - x(t)

y_l {t) = x{t - *i )

z_l t) = x{t - 2r_l )

and

x₂ (t) - x(t

y₂ {t) = x{t

z₂ = where τ, < τ₂. It can be seen from the above equations that the second vector of delay coordinates is related to the first vector of the delay coordinates by the second time delay, τ₂. That is, the second vector of delay coordinates at a first point in time t takes value equal to a value of the first vector of delay coordinates at a second point in time (t-T₂), where the difference between the first and second points in time is equal to the second time delay τ₂. A signal defined by the second vector of delay coordinates can therefore be seen as a time-delayed copy of a signal defined by the first vector of delay coordinates.

As with the embodiments described previously, in step S1602 the plurality of transformed data points are then obtained by transforming each of the first and second vectors of delay coordinates into a coordinate system comprising one or more predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points. Here, the predefined vectors are each a function of both the first and second delay coordinates. In the present embodiment, the vectors Vi(t) and v₂(f) are defined as shown below in order to factor out baseline variation in compound BW signals which exhibit linear baseline drift. In other embodiments different vectors may be defined in order to factor out a different type of baseline drift, for example, suitable vectors may be defined which are invariant for signals with quadratic drift.

L (2 _! (f ) - y, (f ) - z, (f ) - 2x₂ (f ) + y₂ (f ) + z₂ (f )),

In the present embodiment, the smaller time delay ii is obtained first by choosing a fixed value for τ₂. In step S1603, a suitable fixed value for the second time delay τ₂ is selected to lie within a range of expected values for the type of signal under analysis. Then, in step S1604 the first time delay Ti is found by minimising the value of the mean of

(t) for varying values of the first time delay τ, while keeping the second time delay τ₂ constant. In other embodiments a different function may be optimised in order to obtain the optimised value of the first time delay Tj. For example, any of the above-described methods for determining the optimum value of a time delay can be used to optimise the first time delay T whilst holding the second time delay τ₂ constant. If a minima is not found, the process can be repeated using a different fixed value of the second time delay τ₂.

In some embodiments, if it is not known whether the input signal contains more than one frequency component, a system such as the one shown in Fig. 8 may be configured to repeatedly attempt to find a value of the first time delay using different fixed values of the second time delay. For example, the delay coordinate unit 801 can calculate an initial set of delay coordinates using a first fixed value of τ₂, and the repeat period calculating unit 803 can then attempt to find a value of τ,. If after a certain number of attempts it has not been possible to determine the first time delay, for example after ten unsuccessful attempts, the repeat period calculating unit 803 can conclude that the signal only contains one frequency component and the system 800 can revert to one of the above-described methods.

Once a value for the first time delay τ, has been determined, the second time delay τ₂ is then determined. In the present embodiment, a different method of determining the second time delay τ₂ is used if the signal type is known in advance. Hence in step S1605 it is checked whether the type of compound signal, for example BW, AM or FM, is known. If the signal type is known, then in step S1606 a suitable algorithm can be selected for determining the second time delay τ₂. On the other hand, if the signal type is not known, then in step S1607 a general-purpose algorithm can be used to determine the second time delay τ₂.

Once the second time delay τ₂ has been determined, the first and second repeat periods are then calculated in step S1608 based on a known relationship between the first repeat period T, and the optimised first time delay τ,, and a known relationship between the second repeat period T₂ and the optimised second time delay τ₂. In the present embodiment, the first repeat period T_t is equal to twice the optimised first time delay Τι (i.e.

and the second repeat period T₂ is equal to twice the optimised second time delay τ₂ (i.e.

In other embodiments the relationship between the first and second repeat periods and the first and second time delays may be different, depending on the exact method that is used to derive the optimised first and second time delays.

The first and second repeat periods maybe approximately, but not exactly, equal to the optimised first and second time delays respectively. Depending on the accuracy that is required in a particular application, the first and second time delays may be taken to be the first and second repeat periods. In other embodiments a more accurate value for the first and/or second repeat period may be determined using one of the above- described methods, once the first and second time delays have been optimised. An example of a general-purpose algorithm is shown in Fig. 17. A method such as the one illustrated in Fig. 17 can be used to determine the second time delay for an unknown type of compound signal. First, in step S1701 a time delay transformation is applied to the plurality of samples of the periodic data to obtain a plurality of transformed signal samples. Each transformed signal sample w is obtained from a difference in values between one of the plurality of samples of the periodic data at a first point in time and one of the plurality of samples of the periodic data at a second point in time t₂, as follows: w - x{t_\ ) - x(t₂ ) where t₂ = ( -Ti-H), H being a shift parameter. The difference between the first and second points in time is therefore set based on the optimised value of the first time delay τ,. In the present embodiment the difference between the first point in time and the second point in time t₂ is not exactly equal to T but is offset by an amount equal to the shift parameter H. If a signal consists of only one period and Τι = Tt exactly, where T_t is the first repeat period, then if the shift parameter H is set to zero then the w coordinate will be exactly zero. On the other hand, if the first time delay Τι differs slightly from the first repeat period ΤΊ then the w coordinate will not lie exactly on the y axis. In practice, when discretised sampling is used, the exact value of the first repeat period T_t is very unlikely to lie on one of the discrete values that the first time delay Τι may take by the discretised sampling. For all types of compound signal, the amount the w coordinate is away from the horizontal is related to the second slower frequency component of the signal. The use of a shift parameter H causes the w coordinate to have greater deviation from the horizontal, enabling easier detection of the slower period. However, in other embodiments the shift parameter could be omitted.

This transformed signal w is sensitive to outliers, and so in the present embodiment the absolute values of the w signal are smoothed using a moving average. Also, to enable reliable detection of AM and FM signals, logs are taken of the smoothed signal. This simplifies the analysis for AM and FM signals, which can be represented as the product of two frequencies A and B, since log(AB) can be decomposed into a simple sum log(A) + log(B). However, in other embodiments the step of taking logarithms may be omitted. Also, in order to enable logs to be taken, the signal w is first squared to ensure that the signal only takes positive values, which results in the period being halved. As such, in embodiments in which the signal w is squared, the repeat period calculated from the optimised second time delay τ₂ must be doubled to find the original slow repeat period. The signal defined by the logs of the smoothed signal squared is hereinafter denoted as the signal M.

Next, in step S1702 the first and second vectors of delay coordinates are recalculated using the transformed signal samples M obtained using the method described above in relation to step S1601. Then, in step S1703 the predefined vectors Vi and v₂, as defined above in relation to step S1602, are recalculated using the first and second vectors of delay coordinates obtained in step S1702.

Finally, in step S1704 the optimised value of the second time delay τ₂ is determined from the recalculated predefined vectors Vi and v₂. In the present embodiment the optimised value of the second time delay τ₂ is determined by minimising the function:

C = v₁ ² (i)+ v₂ ² (i)

However, in other embodiments a different function could be chosen for optimising the second time delay τ₂. For example, in some embodiments the signal C given by the equation shown above may be smoothed, and the second time delay τ₂ can be determined as the position of a minimum in the smoothed signal.

The signal C given by the above equation is periodic. If the signal C has been derived by taking logs of absolute values as described above in step S1701, then C is periodic with a single period of T2/2. Alternatively, if the signal C is derived directly from the transformed signal w, then C may be periodic with a repeat period equal to T₂. Since the signal C is itself periodic, in some embodiments the signal C can be further processed using a time-delay transformation technique to determine the optimised second time delay. For example, if C has a repeat period of T2/2, a transformed signal Q can be obtained as follows: e(_f) = C(f) - C(f - 2r₂ ) This transformation is similar to that used in step S1701 to derive w. The mean of Q(t) can then be minimised to find the second repeat period τ₂. One method by which Q(t) can be minimised is as follows. First, a search is carried out for peaks in the signal Q. For example, all peaks higher than a certain threshold maybe identified, or all peaks with a spacing of at least τ₂/4 may be identified. If more than a certain number of peaks are detected, for example more than nine peaks, the peaks at one or both ends of a plot of the signal Q may be discarded since these are typically the least accurate. The distances between any remaining peaks is then calculated and stored. The distance between two peaks is related to the second time delay τ₂. For example, the second time delay τ₂ may be calculated by multiplying the mean distance between peaks in the signal Q by four, or by multiplying the distance to the first peak by four. The reason for using a multiplication factor of four in the present embodiment is as follows. When taking the absolute values of w, the period of the resulting signal is halved with respect to the original signal. Furthermore, the three point transformation used to derive the signal Q from the signal C has minima (i.e. inverted peaks) with a spacing equal to half the second time delay (τ₂/2). Since the slow repeat period T₂ = ΐτ₂, the peak spacing in Q is therefore approximately equal to VAT₂. Thus, any value for the second repeat period T₂ obtained from the peak spacing in Q must be multiplied by four in this embodiment.

Embodiments such as the ones described above with reference to Figs. 16 and 17 can therefore be used to determine the repeat periods of two different frequency components in a compound signal, whether or not the particular type of compound signal is known. Furthermore, in some embodiments an n=2 or 71=3 method such as those described above with reference to Fig. 4, or a similar method using a different embedding dimension, may be applied to any transformed periodic variable (e.g. v, w, Q to determine an individual repeat period.

Whilst certain embodiments of the invention have been described herein with reference to the drawings, it will be understood that many variations and modifications will be possible without departing from the scope of the invention as defined in the

accompanying claims.

Claims

1. Apparatus for determining a repeat period in periodic data, the apparatus comprising:

a delay coordinate unit configured to obtain a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window;

a coordinate transforming unit configured to obtain a plurality of transformed data points by transforming each of the vectors into a coordinate system comprising a plurality of predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points; and

a repeat period calculating unit configured to determine the repeat period based on the transformed data points.

2. The apparatus of claim l, wherein the delay coordinate unit is configured to obtain each vector of delay coordinates based on a time delay τ, the predefined vectors are selected such that a transformed data point for an exactly periodic signal lies on one of the predefined vectors( v ) and the magnitude of another one of the predefined vectors ( w ) is zero when a relationship between τ and the repeat period is satisfied, and the apparatus further comprises:

a time delay optimising unit configured to determine an optimised value of τ for which the spread of the transformed data points around said one of the predefined vectors ( ) is minimised, by finding the value of τ which minimises a norm of the values of said other one of the predefined vectors ( w ) for the plurality of transformed data points,

wherein the repeat period calculating unit is configured to determine an average repeat period in the approximately periodic data based on the optimised value of τ.

3. The apparatus of claim 2, wherein the periodic data is blood pressure data, and the apparatus is further configured to determine a density in said one of the predefined vectors ( v ) as a measure of time, over a plurality of repeat periods in the periodic data, and determine a measure of cardiac contractility based on the area under a peak in the density, or based on the magnitude of the peak.

4. The apparatus of claim 1, wherein the delay coordinate unit is configured to obtain each vector of delay coordinates based on a first time delay τΐ, and is further configured to obtain a second delay coordinate from one of said one or more predefined vectors based on a second time delay τ2,

wherein the time delay optimising unit is further configured to determine an optimised value of the second time delay τ2 which minimises a function related to the difference between said one of the predefined vectors and the second delay coordinate, and

wherein the repeat period calculating unit is further configured to determine an average repeat period in the approximately periodic data from the optimised value of the second time delay τ2.

5. The apparatus of claim 4, wherein the function related to the difference between said one of the predefined vectors and the second delay coordinate is a scale invariant function.

6. The apparatus of claim 1, wherein the plurality of transformed data points define a projection of an attractor of the periodic data along an axis normal to the plurality of predefined vectors, and

wherein the repeat period calculating unit is configured to determine the repeat period by determining a time taken to complete one loop of the obtained attractor in the projection along said axis.

7. The apparatus of claim 6, wherein the repeat period calculating unit is configured to determine a first time point as being a time at which a trace of the periodic data crosses a section arranged to cut one side of the attractor in the projection along said axis, determine a second time point as being a time at which the trace of the periodic data next crosses the section, and determine the time taken to complete one loop as the difference in time between the first time point and the second time point.

8. The apparatus of claim 7, wherein the repeat period calculating unit is configured to determine a time at which the trace of the periodic data crosses the section by interpolating between two consecutive data points falling on opposite sides of the section.

9. The apparatus of claim 7 or 8, wherein the repeat period calculating unit is configured to arrange the section to be perpendicular to said one side of the attractor.

10. The apparatus of any one of claims 6 to 9, wherein the delay coordinate unit is configured to obtain three-dimensional vectors of delay coordinates [x, y, z], and the axis normal to the predefined vectors is x = y = z.

11. The apparatus of claim 1, wherein the delay coordinate unit is configured to obtain first and second vectors of delay coordinates for each one of the plurality of samples of the periodic data in the time window, the first vector of delay coordinates being obtained using a first time delay and the second vector of delay coordinates being obtained using a second time delay longer than the first time delay,

wherein the apparatus further comprises:

a time delay optimising unit configured to obtain an optimised value of the first time delay while holding a value of the second time delay constant, and to obtain an optimised value of the second time delay using the optimised value of the first time delay, and

wherein the repeat period calculating unit is configured to determine a first repeat period from the optimised value of the first time delay and to determine a second repeat period from the optimised value of the second time delay.

12. The apparatus of claim 11, wherein the time delay optimising unit is configured to obtain the optimised value of the second time delay by applying a time delay transformation to the plurality of samples of the periodic data to obtain a plurality of transformed signal samples, each transformed signal sample being obtained from a difference in values between one of the plurality of samples of the periodic data at a first point in time and one of the plurality of samples of the periodic data at a second point in time, the time delay optimising unit being configured to set the difference between the first and second points in time based on the optimised value of the first time delay;

wherein the delay coordinate unit is configured to recalculate the first and second vectors of delay coordinates for the transformed signal samples,

wherein the coordinate transforming unit is configured to transform each of the recalculated vectors into the coordinate system comprising the plurality of predefined vectors, and

wherein the time delay optimising unit is configured to obtain the optimised value of the second time delay from said plurality of predefined vectors.

13. A method of determining a repeat period in approximately periodic data, the method comprising:

obtaining a vector of delay coordinates for each one of a plurality of samples of the periodic data in a time window;

obtaining a plurality of transformed data points by transforming each of the vectors of delay coordinates into a coordinate system comprising one or more predefined vectors arranged to factor out baseline variation in the periodic data from the transformed data points; and

determining the repeat period based on the transformed data points.

14. The method of claim 13, wherein each vector of delay coordinates is obtained based on a time delay τ, the one or more predefined vectors are arranged such that a transformed data point for an exactly periodic signal lies on one of the predefined vectors ( v ) and the magnitude of another one of the predefined vectors ( w ) is zero when a relationship between τ and the repeat period is satisfied, and the repeat period for the approximately periodic data is determined by:

determining an optimised value of τ for which the spread of the transformed data points around said one of the predefined vectors ( ) is minimised, by finding the value of τ which minimises a norm of the values of said other one of the predefined vectors ( w ) for the plurality of transformed data points; and

determining an average repeat period in the approximately periodic data from the optimised value of τ.

15. The method of claim 14, wherein the periodic data is blood pressure data, and the method further comprises:

determining a density in said one of the predefined vectors ( v ) as a measure of time, over a plurality of repeat periods in the periodic data; and

determining a measure of cardiac contractility based on the area under a peak in the density, or based on the magnitude of the peak.

16. The method of claim 13, wherein each vector of delay coordinates is obtained based on a first time delay τΐ, the method further comprising:

obtaining a second delay coordinate from one of said one or more predefined vectors based on a second time delay τ2; determining an optimised value of the second time delay τ2 which minimises a function related to the difference between said one of the predefined vectors and the second delay coordinate; and

determining an average repeat period in the approximately periodic data from the optimised value of the second time delay τ2.

17. The method of claim 16, wherein the function related to the difference between said one of the predefined vectors and the second delay coordinate is a scale invariant function.

18. The method of claim 13, wherein the plurality of transformed data points define a projection of an attractor of the periodic data along an axis orthogonal to the one or more predefined vectors, where movement along said axis corresponds to baseline variation in the periodic data, and

wherein the repeat period is determined by determining a time taken to complete one loop of the obtained attractor in the projection along said axis.

19. The method of claim 18, wherein determining the time taken to complete one loop of the obtained attractor comprises:

determining a first time point as being a time at which a trace of the periodic data crosses a section arranged to cut one side of the attractor in the projection along said axis;

determining a second time point as being a time at which the trace of the periodic data next crosses the section; and

determining the time taken to complete one loop as the difference in time between the first time point and the second time point.

20. The method of claim 19, wherein a time at which the trace of the periodic data crosses the section is determined by interpolating between two consecutive data points falling on opposite sides of the section.

21. The method of claim 19 or 20, wherein the section is arranged to be

perpendicular to said one side of the attractor.

22. The method of any one of claims 17 to 20, wherein the obtained vectors of delay coordinates are three-dimensional vectors [x, y, z], and the axis normal to the predefined vectors is x = y = z.

23. The method of claim 13, wherein obtaining the vector of delay coordinates comprises obtaining first and second vectors of delay coordinates for each one of the plurality of samples of the periodic data in the time window, the first vector of delay coordinates being obtained using a first time delay and the second vector of delay coordinates being obtained using a second time delay longer than the first time delay, the method further comprising:

obtaining an optimised value of the first time delay while holding a value of the second time delay constant; and

obtaining an optimised value of the second time delay using the optimised value of the first time delay,

wherein a first repeat period is determined from the optimised value of the first time delay and a second repeat period is determined from the optimised value of the second time delay.

24. The method of claim 23, wherein the optimised value of the second time delay is obtained by:

applying a time delay transformation to the plurality of samples of the periodic data to obtain a plurality of transformed signal samples, each transformed signal sample being obtained from a difference in values between one of the plurality of samples of the periodic data at a first point in time and one of the plurality of samples of the periodic data at a second point in time, a difference between the first and second points in time being set based on the optimised value of the first time delay;

recalculating the first and second vectors of delay coordinates for the transformed signal samples;

transforming each of the recalculated vectors into the coordinate system comprising the plurality of predefined vectors; and

obtaining the optimised value of the second time delay from said plurality of predefined vectors.

25. The method of any one claims 13 to 24, wherein the periodic data is

physiological data.

26. A non-transitory computer-readable storage medium arranged to store software instructions which, when executed on one or more processors, perform the method of any one of claims 13 to 26.

27. Apparatus for determining a repeat period in periodic data, the apparatus comprising:

processing means for executing computer program instructions; and

computer-readable memory arranged to store computer program instructions executable by the processing means, wherein the stored computer program instructions include instructions for performing the method of any one of claims 13 to 26.