Methods, computer program products, computer-readable storage media, and systems

Abstract

A computer-implemented method is provided to estimate a periodic feature, such as a period length or a period location, of a physiological rhythm, such as a cardiac rhythm or an autonomic respiratory rhythm. The computer-implemented method comprises receiving a signal representing the physiological rhythm, estimating a first set of periodic feature estimates based on the signal, generating a prior probability of the periodic feature based on at least a subset of the first set of periodic feature estimates, and estimating a second set of periodic feature estimates based on the first set and the prior probability. By using such a computer-implemented method, the accuracy of the signal is improved, in particular when the signal is obtained from an accelerometer arranged on the chest of a subject.

Description

TECHNICAL FIELD

[0001] The present invention relates to a method of processing a signal representing a circadian rhythm of a subject. The circadian rhythm has a periodic character. The present invention further relates to a computer program product and a computer readable storage medium, each having instructions which, when executed by a computer, cause the computer to perform the method. Furthermore, the present invention relates to a system for processing a signal representing a circadian rhythm of a subject. BACKGROUND

[0002] Accurate knowledge of circadian rhythms, such as heart rate or respiratory rate, is of crucial importance in many healthcare applications, including diagnosis, patient monitoring and treatment. Therefore, estimating the time-varying frequency of these rhythms from conventional clinical measurement modalities, such as electrocardiography (ECG), pulse oximetry or respiratory inductance plethysmography, is a well-studied problem. While for some applications it can be sufficient to estimate only the average frequency over multiple cycles of the underlying rhythm, other applications, such as heart rate variability (HRV) analysis, require resolution between events.

[0003] Due to changes in the population structure and the increasing number of patients with chronic diseases, home monitoring of the health status of patients will play an important role in future healthcare systems. Recent studies have shown that home telehealth approaches can substantially reduce hospitalization and nursing days.

[0004] Since traditional measurement modalities are rarely applicable outside of the clinical environment, there is currently a move towards unobtrusive and unconstrained measurement systems to allow for continuous long-term monitoring of the health status of patients. These systems aim to provide two advantages over current approaches, as they generally do not require user interaction or compliance to perform the measurements, and they do not interfere with the user's daily routine. Unobtrusive measurement systems based on various different principles are known. These include capacitive coupling electrodes for ECG measurements in a chair or car seat, capacitive respiration measurements in a bed, radar, millimeter wave interferometry, and ballistocardiography (BCG) or seismocardiography (SCG) sensors integrated in a bed, chair and weighing scale. Other types of unobtrusive measurement systems include wearable systems. For example, these measurement systems are implemented as a watch, a pendant and a chest belt.

[0005] While these measurement systems provide the above-mentioned advantages over existing methods, they also have their own set of shortcomings. Most notably, due to the fact that the environment in which the measurements are taken is usually not controlled, the signal quality can vary greatly and be unreliable. For certain types of sensors, the morphology of the recorded signal can vary significantly depending on the orientation and position of the user relative to the sensor. This can be observed, for example, in the case of bed-mounted BCG sensors that record the cardiac vibrations caused by the contraction of the heart and the ejection of blood into the aorta. In the case of wearable measurement systems, the signal can vary dramatically with changes in the patient's posture.

[0006] In these cases, a reliable estimate of the instantaneous frequency (or its reciprocal: the local interval length) can pose a significant challenge to algorithms originally developed for clinical applications. These algorithms are usually based on detecting specific features of the signal that are related to the events of interest, such as the QRS complex in an ECG. The instantaneous frequency is then computed by differencing the times of successive occurrences of these events. In the case of an ECG, sophisticated methods have been proposed to improve the robustness of the instantaneous heart rate estimate. However, these methods are usually limited to one type of signal and require extensive prior knowledge of the expected morphology of the waveform. As mentioned above, these requirements do not hold for unobtrusive sensor types. SUMMARY

[0007] It is an object of the present invention to provide an improved method for processing a signal representing a physiological rhythm of a subject. Optionally, it is an object of the present invention to provide a method for obtaining improved information representing a physiological rhythm for different postures of a subject.

[0008] According to a first particular aspect, there is provided a computer-implemented method of processing a signal representing a physiological rhythm of a subject, the physiological rhythm having a periodic characteristic, the method comprising:

[0009] receiving the signal;

[0010] estimating a first set of periodic characteristic estimates based on the signal;

[0011] generating a prior probability for the periodic characteristic based on at least a subset of the first set of periodic characteristic estimates;

[0012] estimating a second set of periodic characteristic estimates based on the first set and the prior probability.

[0013] A first set of periodic characteristic estimates is estimated based on the signal. The first set comprises estimates of the periodic characteristic for a plurality of cycles of the physiological rhythm. In most cases, the first set of periodic characteristic estimates provides good estimates of the periodic characteristic for some cycles, while the first set provides poor estimates of the periodic characteristic for other cycles. For example, the poor estimates are due to movements of the subject or any other type of event that negatively affects the signal accurately representing the physiological rhythm.

[0014] Based on the first set of periodic feature estimates, a prior probability is generated. Based on information from the first set of periodic feature estimates, the prior probability distinguishes between estimates of a possible occurrence of a periodic feature and estimates of an impossible occurrence of a periodic feature. The prior probability is based on the entire first set or on a subset of the first set. For example, the subset comprises less than optimal estimates of periodic features than the first set. For example, the subset set only contains good estimates of periodic features. Physiological rhythms have some predictable characteristics. For example, an upper and lower limit of a heart rate of a subject are known, and the heart rate can only change gradually over time. In another example, an upper and lower limit of an autonomous respiration rhythm of a subject are known. For example, estimates from the first set that exceed these limits are excluded from the subset.

[0015] By estimating a second set of periodic feature estimates based on the first set and the prior probability, the periodic features of the physiological rhythm are more accurately estimated. In particular, for signal portions for which the first set initially provided less than optimal estimates, the accuracy is improved. Thus, even in case of motion or changing posture of the subject, more accurate estimates of the periodic features of the physiological rhythm are provided. This way, accurate estimates of the physiological rhythm are obtained even using an unobtrusive and unconstrained measurement system. Optionally, for signal portions for which the first set initially already provided good estimates, the periodic features of the physiological rhythm are more accurately estimated. In this way, an improved accuracy of the estimates of the periodic characteristics is obtained.

[0016] Physiological rhythms, such as a heartbeat and an autonomous respiration rhythm, have a repetitiveness. The repetitiveness of a physiological rhythm is characterized by a period. A period is defined as the time between two repetitions. For example, a period is the time between two events in a physiological rhythm. For example, a period refers to the time between two heartbeats, or the time between starting two consecutive breaths. A period has at least one periodic feature, such as a period length, a rhythm frequency, or a timing of a certain phase in the period, such as a start or an end of the period. For example, a periodic feature is a physiological event that normally occurs in each period. An example of such a physiological event is the opening of an atrioventricular valve, or the change of airflow from inhalation to exhalation. Due to a disease or a condition of the subject, such a physiological event can not occur in each period.

[0017] A signal comprises information about a physiological rhythm. The signal is generated, for example, by a measurement device that is configured to generate the signal in response to measuring a physiological characteristic of a subject. The measurement device is configured to measure, for example, a heartbeat, a pulse, or a breathing motion or airflow caused by breathing. For example, the signal is an electrical signal. For example, the signal comprises signals from multiple sensors. For example, the multiple sensors measure the physiological rhythm of the subject simultaneously.

[0018] Based on the signal, the method estimates a first set of periodic feature estimates, for example, by determining a local periodicity spectrum.

[0019] Due to the predictable nature of the circadian rhythm, a prior probability is generated based on the first set of periodicity feature estimates. Alternatively, the prior probability is referred to as a prior probability distribution, or simply a prior. The prior probability gives a probability of a value of a periodicity feature estimate in the first set using statistical data.

[0020] For example, the prior probability uses additional information. For example, the additional information includes information about a limitation of the circadian rhythm of the subject, such as a maximum heart rate. For example, the additional information includes information about an activity the subject is performing at the time the signal is generated. For example, the activity is sleeping or sitting or lying or walking, etc.

[0021] For example, the prior probability uses only a subset of the first set of periodicity feature estimates. For various reasons, a portion of the signal can not reflect the circadian rhythm well. For example, during that portion of the signal, the measurement device is placed improperly, or the subject is overly active. The first set of periodicity feature estimates can give very unrealistic estimates based on that portion of the signal. By excluding these unrealistic estimates, the prior probability is improved.

[0022] Based on the first set and the prior probability, a second set of periodicity feature estimates is estimated. By using the prior probability, a probability of a value of a periodicity feature estimate of the first set is assessed over some period. If the prior probability indicates that the probability of the value is high, then the periodicity feature estimate is an accurate estimate. If the prior probability indicates that the probability of the value is low, then the periodicity feature estimate is an inaccurate estimate. In this case, the periodicity feature estimate is re-estimated using the value from the first set and the prior probability. The second set of periodicity feature estimates includes, for example, some values of the periodicity feature estimates from the first set and some re-estimated values of the periodicity feature estimates based on the values of the first set and the prior probability.

[0023] In an example, the second set of periodicity feature estimates is not based on the complete first set, but only on a portion of the first set. For example, the second set is based on a portion of the first set that provides at least a minimum accuracy. For example, the second set is not based on a portion of the first set that has very poor estimates of the periodicity features. The very poor estimates can be due to significant errors in the measurement values that generated the signal. When a signal experiences such severe errors, the signal no longer contains true information representing the circadian rhythm, so that signal is not processed further.

[0024] The method includes determining a spectral likelihood of at least a portion of the signal. The spectral likelihood is a likelihood that the first set of periodicity feature estimates corresponding to the at least a portion of the signal represents a periodicity feature of the circadian rhythm. The method further includes estimating a second set of periodicity feature estimates based on the spectral likelihood of the signal.

[0025] The spectral likelihood converts the signal from the time domain to the frequency domain. The spectral likelihood represents the signal as a series of periodicities. The spectral likelihood compares a periodicity of the signal to a similarity of another periodicity of the signal. For example, the spectral likelihood compares two adjacent periods of the signal. The spectral likelihood represents a likelihood of a value of a period length of a periodicity in the signal. For example, the spectral likelihood represents a most likely value of a period length of a periodicity in the signal. By considering the likelihood of a value of a period length of a periodicity in the signal, an estimate of a periodicity characteristic can be more accurately estimated.

[0026] The method comprises defining a location point in the physiological rhythm. The location point characterizes a phase in a cycle of the physiological rhythm. The method comprises determining a time likelihood of the signal corresponding to the cycle. The time likelihood represents a likelihood that the signal represents the location point. The method comprises estimating the location point based on the spectral likelihood, the time likelihood, and a prior probability.

[0027] To define a timing of a cycle of the physiological rhythm, a location point is defined. The location point characterizes a phase in a cycle of the physiological rhythm. For example, the location point represents a beginning of the cycle or an end of the cycle. For example, the location point indicates a periodic event in the physiological rhythm that causes a significant feature in the signal. For example, the significant feature is a large positive value of the signal or a large negative value of the signal. For example, the significant feature is a maximum gradient in the signal or a zero gradient in the signal. In case the physiological rhythm is a cardiac rhythm, the phase in a cycle of the physiological rhythm is, for example, a heartbeat, or a closing or opening of a heart valve, such as an aortic valve or a pulmonary valve. For example, the phase of the cardiac rhythm causes a large acceleration or a large electrical impulse of the subject, thereby causing the significant feature in the signal. In case the physiological rhythm is an autonomous respiratory rhythm, then the phase in a cycle of the physiological rhythm is, for example, a beginning or end of an inhalation, a beginning of an exhalation, a pause between an inhalation and an exhalation, a pause between an exhalation and an inhalation, a change in direction of airflow, or a pressure caused by airflow from the subject. For example, the phase of the autonomous respiratory rhythm causes a large acceleration or a maximum air pressure or a minimum air pressure of the subject, thereby causing the significant feature in the signal.

[0028] The time likelihood indicates a likelihood that a point in the signal represents the location point in the time domain. By estimating the location point based on the spectral likelihood, the time likelihood, and a prior probability, the location point is accurately estimated. This allows an accurate estimation of a timing of a phase of the physiological rhythm.

[0029] In an embodiment, the location point of the physiological rhythm indicates one of: i) a beginning of the cycle, ii) an end of the cycle, iii) a middle of the cycle.

[0030] In an embodiment, estimating the location point comprises using a maximum a posteriori probability (MAP) estimation. The maximum a posteriori probability (MAP) estimation is based on the spectral likelihood, the time likelihood, and a prior probability.

[0031] The position is estimated by using the MAP estimate, the position is estimated in a robust manner, and the most likely value for the position is provided.

[0032] In an embodiment, the second set of period feature estimates estimates the period feature of the physiological rhythm more accurately than the first set of period feature estimates.

[0033] The period feature is estimated more accurately by the period feature estimates based on the first set and the prior probability.

[0034] In an embodiment, the period feature comprises at least one of a period length and a time stamp.

[0035] The period length indicates the time between two consecutive periods. For example, the period length indicates the time between two heartbeats or the time between the start of two breathing cycles. The time stamp indicates the position of the period in time. The time stamp indicates when in time the period occurs, e.g. when in time a certain feature point of the period occurs. For example, the time stamp provides information about the interbeat interval or the heart rate variability.

[0036] In an embodiment, the physiological rhythm is the autonomous breathing rhythm of the subject. For example, the period feature represents the breathing frequency of the subject.

[0037] In an embodiment, the physiological rhythm is the heart rhythm of the subject. For example, the period feature represents the heart beat frequency of the subject.

[0038] In an embodiment, the method comprises estimating at least one of the interbeat interval and the heart rate variability of the heart rhythm based on the second set of period feature estimates.

[0039] The interbeat interval and the heart rate variability are physiological parameters that provide valuable information about the physical condition of the subject. By accurately estimating one or both of these parameters, accurate information can be obtained to help medical professionals assess the physical condition of the subject.

[0040] In an embodiment, the signal comprises a signal from an accelerometer arranged on the chest of the subject.

[0041] Many people consider an accelerometer placed on the chest of a subject as a comfortable, unobtrusive, and unconstrained way of acquiring information about a physiological rhythm. The accelerometer is configured to generate a signal based on a physiological rhythm, such as a respiratory rhythm or a cardiac rhythm. However, the signal of the accelerometer is also influenced by movements of the subject or the posture of the subject. Furthermore, if the subject moves, the accelerometer can be displaced relative to the subject, resulting in a change of the signal. Due to these influences, the signal from the accelerometer is less accurate than the signal of an ECG. However, by applying the method of the invention to the signal of an accelerometer arranged on the chest of a subject, the accuracy of the signal is significantly improved. This provides an accurate estimate of the periodicity while using a measurement device that is easy to use, i.e. a chest-worn accelerometer.

[0042] In embodiments, the signal represents a physiological rhythm over a period of time, wherein the subject is asleep during the period of time, e.g. asleep in the range of 6-10 hours.

[0043] For example, the period of time is larger than 1 hour, or larger than 2 hours, or larger than 4 hours. During such a long period of time, it is very likely that the subject will move. Such movements can decrease the accuracy of the signal. By applying the method, the signal can be acquired over a long period of time while providing an improved accuracy of the estimate of the periodicity over the long period of time. By acquiring the periodicity information over a long period of time, more information about the physical condition of the subject can be acquired.

[0044] In embodiments, the method comprises determining the occurrence of a sleep disordered breathing (SDB) event based on the second set of estimates of periodicity.

[0045] A sleep disordered breathing (SDB) event is an event in which the subject experiences an interruption in breathing during sleep. SDB includes apneas, such as obstructive sleep apnea (OSA) or central sleep apnea (CSA). During an SDB event, e.g. breathing is paused, obstructed, or significantly reduced. For example, an SDB event is a hypopnea. During a hypopnea event, the airflow of the subject caused by breathing is significantly reduced, such as by at least 30%.

[0046] Due to the SDB event, the oxygen uptake of the subject is reduced and / or the carbon dioxide emission of the subject is reduced. Due to the SDB event, one or more physiological rhythms of the subject are affected. By determining whether an SDB event occurs using the method described in the invention, it can be determined more accurately whether an SDB event occurs. For example, this allows using a measurement device with a lower accuracy that is suitable for use at home while still being able to determine the occurrence of an SDB event with sufficient accuracy.

[0047] In embodiments, the method comprises determining a sleep stage of the subject based on the second set of estimates of periodicity.

[0048] Sleep has four stages: three non-rapid eye movement (NREM) sleep stages and one rapid eye movement (REM) sleep stage. The non-rapid eye movement (NREM) sleep stages are N1, N2 and N3. Sleep stages N1 and N2 are light sleep, while sleep stage N3 is deep sleep. Different sleep stages have certain differences in the subject's circadian rhythm. For example, the heart rate is slower during N3 than during N1 or N2. For example, the heart rhythm during N3 is more regular than during REM, the respiration during N2 and N3 is slower than during N1. For example, during REM, the respiration is more irregular than during N1 -N3. Based on the periodic features estimated according to the method of the application, the sleep stage of the subject is determined. For example, the duration and frequency of each sleep stage are determined.

[0049] In an embodiment, the signal comprises a sum of three-dimensional acceleration data, wherein the sum of three-dimensional acceleration data is band-pass filtered, for example with a band-pass filter having cut-off points of 3 Hz and 25 Hz.

[0050] By using three-dimensional acceleration data, the signal is less affected by the subject's posture or the orientation of the measuring device with respect to the subject. To make the signal more concentrated on the acceleration caused by the circadian rhythm, the three-dimensional acceleration data is band-pass filtered. Acceleration below a certain cut-off point and acceleration above a certain cut-off point are less likely to be caused by the circadian rhythm, and are thus filtered out of the three-dimensional acceleration data. For example, a band-pass filter with cut-off points of 3 and 25 Hz provides a good filtering effect. By using this band-pass filter, the three-dimensional acceleration data mostly or only contains data with a frequency between 3 Hz and 25 Hz.

[0051] In an embodiment, the method comprises receiving information about the subject. The generation of the prior probability is based on the information about the subject.

[0052] By including information about the subject, the prior probability is further improved. This information is provided in addition to the first set of periodic feature estimates. This information includes the subject's gender, age, weight or body mass index, etc. This information includes, for example, information about the subject's physical condition or about previous circadian rhythm measurements or about previous medical examinations.

[0053] In a second aspect of the application, a computer program product is provided, comprising instructions which, when executed by a computer, cause the computer to perform the computer-implemented method of the first aspect of the application.

[0054] In a third aspect of the application, a computer-readable storage medium is provided, the medium comprising instructions which, when executed by a computer, cause the computer to perform the computer-implemented method of the first aspect of the application.

[0055] In a fourth aspect of the application, a system for processing a signal representing a physiological rhythm of a subject is provided, the physiological rhythm having a periodic characteristic, the system comprising a sensor and a computer, wherein the computer is coupled to the sensor, wherein the computer is configured to perform the computer-implemented method of the first aspect of the application. BRIEF DESCRIPTION OF DRAWINGS

[0056] The exemplary embodiments will now be described, by way of example only, with reference to the following drawings in which:

[0057] Figure 1 A first embodiment according to the application is described;

[0058] Figure 2 Acceleration data in a signal according to the first embodiment is described;

[0059] Figure 3 A graph representing a display of a pre-processed signal according to the first embodiment is depicted;

[0060] Figure 4 A spectrogram of periodicity according to the first embodiment is described;

[0061] Figure 5 Spectral likelihood and temporal likelihood according to the first embodiment are described;

[0062] Figure 6 A representation of the computation of the posterior probability according to the first embodiment is described;

[0063] Figure 7 A representation of the likelihood and probability computation is described;

[0064] Figure 8 Various graphs representing steps according to the first embodiment are shown;

[0065] Figure 9 A second embodiment of a method according to the application is described;

[0066] Figure 10 A system for processing a signal representing a physiological rhythm of a subject according to a third embodiment of the application is described. DETAILED DESCRIPTION

[0067] Figure 1 An embodiment of a method according to the application is described. In this method, a signal 100 representing a physiological rhythm of a subject 1001 is processed. The subject 1001 is as described above in relation to the first embodiment of the application. Figure 10The circadian rhythm has periodic characteristics. The method comprises receiving a signal 100; estimating a first set of periodic characteristic estimates based on the signal 100; generating a prior probability for the periodic characteristics based on at least a subset of the first set of periodic characteristic estimates; and estimating a second set of periodic characteristic estimates based on the first set and the prior probability.

[0068] In this embodiment, the circadian rhythm is the cardiac rhythm. The signal 100 is acquired by an accelerometer 1002 attached to the chest of the subject 1001, see Figure 10 In another embodiment, the circadian rhythm is the autonomous respiratory rhythm of the subject 1001.

[0069] The periodic characteristics include the period length and the time stamp.

[0070] The method makes use of a latent variable model 108. A latent variable model 108 is a statistical model that links observed variables to latent (unobserved) variables. The latent variable θ is the periodic characteristic. The observed variable is the acceleration data D in the signal 100. To estimate the latent variable θ from the observed data D, the maximum a posteriori estimate or MAP estimate is used. The purpose of the MAP estimate is to find the most likely value of the heartbeat θ given the acceleration data D:

[0071] θ MAP = argmax θ {P(θ|D)} (1)

[0072] P(D|θ) is the posterior, the probability of θ given D, which is calculated using Bayes’ theorem:

[0073]

[0074] where P(D|θ) is the likelihood, P(θ) is the prior, and P(D) is the evidence. More information on these formulas can be found in Lambert, B. (2018), A Student’s Guide to Bayesian Statistics, Sage Publications. The advantage of using the MAP estimate is that it allows the use of a prior probability during the estimation, so that it is possible to prevent estimates that are not possible to occur, making the method more robust.

[0075] The method as shown in Figure 1 comprises the following: signal 100, pre-processing step 102, spectral conversion step 104, time conversion step 106, latent variable θ, latent variable model 108, prior probability P(θ), computation step 110, MAP search step 112, update step 114 of the latent variable θ. More explanation of the individual parts of the method outlined in Figure 1 follows, with reference to the other figures.

[0076] Figure 2 The acceleration data D in the signal 100 according to the first embodiment is described as a function of time t. The acceleration data D shows a part of the complete signal 100. The part is about 6 seconds in time, from 6 hours 15 minutes 7 seconds to 6 hours 15 minutes 13 seconds. The complete signal 100 represents a circadian rhythm over a period of time in which the sleep time of the subject 1001 is between 6 and 10 hours.

[0077] Upon receiving the signal 100, the signal 100 is pre-processed in a pre-processing step 102. The signal 100 comprises the raw three-dimensional acceleration data x, y and z, which are summed in the pre-processing step 102 to obtain a one-dimensional signal a:

[0078] a = x + y + z. (3)

[0079] Subsequently, the one-dimensional signal a is band-pass filtered to obtain a pre-processed signal s, with cut-off points of 3 Hz and 25 Hz, respectively, and a sampling frequency of 250 Hz for the one-dimensional signal a.

[0080] Figure 3 Two periods 301, 302 of the pre-processed signal s are shown. The periods 301, 302 are estimated by centering the pre-processed signal s at t m The periodicity spectrum is determined as a function of T. T is the period length. A peak is found in the periodicity spectrum at T = 1130 ms, corresponding to the maximum similarity between s[t-T] and s[t], as shown by superimposing the signal 303 in period 302 of the pre-processed signal s on period 301, i.e. signal 304. In this way, the period length T is estimated for period 301. So in period 301, the estimated period length in the first group is 1130 ms. This procedure is repeated for other periods in the pre-processed signal s. This results in a first group of period feature estimates based on the signal 100.

[0081] In a spectral transformation step 104, a local periodicity spectrum is computed using the pre-processed signal s:

[0082] S m (T) = p(s[t-T], s[t]) (3)

[0083] where T is the estimated period length in Figure 3 s is the pre-processed signal and p is a periodicity function. The periodicity function p represents the similarity between s[t-T] and s[t] for t m < t < t m + T, where t mThe center time of the spectrum. The periodicity function is based on a modified version of the autocorrelation, the average magnitude difference function (AMDF) and the maximum magnitude pair. This calculation can be found in Brüser, C., Winter, S. and Leonhardt, S. (2013) "Robust interbeat interval estimation in cardiac vibration signals", Physiological Measurement, 34(2), pages 123-138.

[0084] The periodicity spectrum is calculated every 200 ms, resulting in a periodicity spectrogram 400 as shown in Figure 4 The spectrogram 400 describes a function of the estimated period length over time. The spectrogram 400 is relatively stable over interbeat intervals and interbeat intervals of the same length. When the interbeat interval changes, a change occurs, i.e. a heartbeat occurs. Around these changes, the periodicity gradually changes, which is visible in the spectrogram 400 as artifacts 401 and 403 at t = 06:15:08 and t = 06:15:13, respectively. In contrast, at 402 at about t = 06:15:12, the spectrogram 400 is relatively stable over 2.5 seconds, while the interbeat interval is 1250 ms. In the middle of this interval, a heartbeat occurs, and the interbeat interval before and after this heartbeat is 1250 ms long.

[0085] The spectrum conversion step 104 improves the accuracy of the determination of the interbeat interval, which is the length T of the period between consecutive heartbeats. To further improve the accuracy of the heartbeat timestamps, i.e. the heartbeat localizations, a time conversion step 106 is performed. The purpose of the time conversion step 106 is to localize the heartbeats at a specific phase in the cardiac cycle. By using a certain phase for each period, the time difference between the timestamps corresponds to the time difference between the R-peaks in the ECG signal. In the present embodiment, the maximum magnitude point of the acceleration data D in the pre-processed signal s is used.

[0086] The method comprises defining a localization point in the physiological rhythm, wherein the localization point represents a phase in a period of the physiological rhythm. The method comprises determining a time likelihood of the signal 100 corresponding to the period. The time likelihood represents a likelihood of the signal 100 representing the localization point. The method comprises estimating the localization point based on the spectrum likelihood, the time likelihood and a prior probability.

[0087] For example, the localization point of the physiological rhythm represents a start of a period of the physiological rhythm, an end of a period of the physiological rhythm or a middle of a period of the physiological rhythm.

[0088] The latent variable model 108 computes the degree of fit of the data to the latent variables θ in the form of a likelihood P(D|θ). The latent variable model 108 takes as input the values of the latent variables (heartbeat timestamps and inter-beat intervals) and computes the degree of fit of the spectral and time-converted data to these timestamps and intervals through the following relation:

[0089]

[0090] where f is the likelihood function, D T is the time-converted data, D S is the spectral-converted data, f T is the time-likelihood function, and f S is the spectral-likelihood function. The first equation states that, without loss of generality, the likelihood is a function of D and θ, while the second equation, which factorizes this function, is our modeling choice.

[0091] The total spectral likelihood f Stot is then computed as follows:

[0092] where the latter equation decomposes the likelihood into a product of functions (f S ), which is not unreasonable since the periodic spectrogram (D S ) depends only locally on the occurrence of heartbeats. The likelihood of a heartbeat timestamp t i and associated inter-beat interval T i is computed as follows:

[0093] f S (D S , t i , T i ) = (Π m S m (T i )) 1 / p (6)

[0094] for all m with t i ≤ t m ≤ t i + T i and p is the number of spectra used. The computation includes the multiplication of all spectra with center time t i in the inter-beat interval (t i to t i + T m ), so that (6) resembles a joint probability. The exponent 1 / p represents a geometric mean, which prevents short interval times (p small) from being favored over long interval times (p large).

[0095] The method comprises determining a spectral likelihood of at least a portion of the signal 100. The spectral likelihood represents a likelihood that a first set of period feature estimates corresponding to the at least a portion of the signal 100 represent a period feature of a physiological rhythm. The method comprises estimating a second set of period feature estimates 100 based on the spectral likelihood of the signal. The second set of period feature estimates more accurately estimate the period feature of the physiological rhythm than the first set of period feature estimates.

[0096] Figure 5 In the upper part of the figure, two inter-beat intervals 501 and 502 are shown as horizontal lines. The first interval 501 has a duration of 880 ms. The first interval 501 covers 4 vertical lines, i.e. 4 spectral lines (p = 4), all having a peak at 880 ms. Hence, the spectral likelihood (6) is maximal at T = 880 ms. The second interval 502 has a duration of 1225 ms. The second interval 502 covers 6 vertical lines, i.e. 6 spectral lines, all having a peak at 1225 ms. Hence, the spectral likelihood is also maximal at T = 1225 ms.

[0097] For the total time likelihood f Ttot , the following is a decomposition method similar to the one used for the spectral likelihood:

[0098]

[0099] where the latter equation decomposes the likelihood into a product of functions f T equal to the value of the time likelihood at t i :

[0100] f T (D T , t) = llh(t) (8)

[0101] In Figure 5 , the following figure is the pre-processed signal s, split into two parts 503 and 504. Part 503 shows the positive values of the pre-processed signal s in Figure 3 . Part 504 shows the negative values of the pre-processed signal s in Figure 3 , but represented as positive values. Figure 5 Two consecutive heartbeat timestamps, vertical dashed lines 505 and 506, are shown. The probabilities llh(t i ) and llh(t i + T i ) are located at the maximum values of the pre-processed signal s in part 503 at lines 505 and 506, respectively. Since the spectral likelihood and the time likelihood of the two heartbeats are at maximum values in Figure 5 , the total likelihood (6) is also at maximum value.

[0102] The prior gives information about the probability of a value of the latent variable θ and helps to prevent silly estimates, e.g. a heart beat interval of 100 ms (corresponding to 600 bpm). The prior uses knowledge of the distribution of the latent variable θ over the measurement population. The prior is explicit, e.g. by using the knowledge that the subject 1001 is a human being, e.g. by using the knowledge that the subject 1001 is lying down, so that shorter heart beat intervals (due to a higher heart rate) are discouraged. Alternatively, the prior considers the probability of a sequence, e.g. by using a Markov model or information from a neural network model. For example, the prior is a personalized prior.

[0103] For example, the prior is a Gaussian prior where μ = μ(t). The value of σ is empirically set to be twice the standard deviation of the difference between the input and output of the median filter length in the first group.

[0104] Based on the prior likelihood, the spectral likelihood and the temporal likelihood, the posterior probability is computed and maximized over the latent variable θ. To reduce the computational complexity of the problem, the approach is to decompose the problem into smaller sub-problems and to solve each sub-problem iteratively. Assume that there is a partial solution θ n :

[0105] θ n = (t n , T n-1 ) = (t1,..., t n , T1,..., T n-1 ) (9)

[0106] θ1= t1, see Figure 6 . Extend θ n to one heart beat interval and know the posterior probability of θ n = P(θ n | D), then the posterior probability of the extended solution θ n+1 can be computed as follows:

[0107] P(θ n+1 | D) = P(θ n , T n | D) = P(T n | D, θ n ) P(θ n | D) (10)

[0108] The first equation uses (9) and t n+1 = t n + T nAnd the second equation uses P(A, B) = P(A|B)P(B). Formula (10) is the basis of the MAP search engine because it is an expression for estimating the optimization posterior probability using MAP. The right side of (10) contains the posterior P(T) n |D,θ n In calculation step 110, Bayes' theorem is used to calculate the following:

[0109]

[0110] We discard P(D) because it is constant when maximizing the posterior, leaving only the numerator. The numerator includes the prior P(T) calculated based on the known distribution. n |θ n The likelihood P(D|T) is calculated using the spectral likelihood function and time likelihood function derived from equations (6) and (7). n θ n ), to obtain:

[0111] P(D|T n θ n )=f S (D S , t n T n )f T (D T , t n+1 (12)

[0112] The calculation of the molecule in (11) is as follows: Figure 7 As shown, Figure 7 The likelihood and probability in the logarithmic field are shown, such that the multiplication in (11) and (12) becomes addition, because if A = BC, then log(A) = log(B) + log(C)). Figure 7 The top figure 700 shows the time log-likelihood f. T Figure 702 in the middle shows the spectral log-likelihood f. S With logarithmic prior P(T) n |θ n The sum of ) and the bottom figure 703 is the posterior P(T) n |D,θ n The logarithm of ).

[0113] In the expansion part, the solution θ n When T is not considered n All values, but only consider T in the posterior probability (13) where there is a peak. n This value allows for the discretization of the search space. Figure 7 This is outlined in the text. Figure 7 Arrows 704, 705, and 706 represent peak values.

[0114] The forward recursive evaluation with constraints starts from the start time t1 and forms the first partial solution θ1. The posterior probability of θ1 is computed as:

[0115] P(θ1|D) = f T (D T ,t1) (13)

[0116] which is the initial value in (13). Starting from the first partial solution, a tree search is performed, where branches are peaks in the posterior (as shown in Figure 7 ) and leaves are partial solutions θ d where d is the forward depth. For example, d = 8 for a good balance between performance and computational load. During the recursive tree exploration, the posterior probability P(θ d |D) in this process is built up by the loop application of (11). Then, the path with the highest posterior probability is selected and θ1 is expanded with the first beat interval of the selected path, so that θ2 is obtained. This process is repeated until the partial solution covers all data and becomes the full solution.

[0117] Equation (11) defines the posterior probability P(T n |D, θ n ) for a single beat interval T n . The posterior probability of the numerator of (11) is approximated as:

[0118] Q(T n ) = P(D|T n , θ n ) P(T n | θ n ) = f S (D S , t n , T n ) f T (D T , t n+1 ) P(T n | θ n ) (14)

[0119] The method includes estimating the time stamps using maximum a posteriori (MAP) estimation. The maximum a posteriori (MAP) estimation is based on the spectral likelihood, the temporal likelihood, and the prior probability. The MAP estimation is computed in a MAP search step 112.

[0120] The result of the MAP search step 112 is used to update the latent variable θ in a computation step 114.

[0121] Figure 8 Various figures are shown that represent steps in the method according to the invention. The figures are functions of time t.

[0122] The first graph 801 starts from the raw acceleration data D, x, y and z. These lines are plotted as visualized adjusted offsets. The peak-to-peak acceleration is about 50 mg.

[0123] The second graph 802 shows the pre-processed one-dimensional signal. This is the pre-processed signal s.

[0124] The third graph 803 shows the spectral transform of the pre-processed signal s.

[0125] The fourth graph 804 shows the time transform of the pre-processed signal s.

[0126] The fifth graph 805 shows the results of the MAP estimation. The MAP estimation provides the most likely values of the latent variable θ, i.e. the inter-beat-intervals (IBI) in [ms] as horizontal lines 807 and the beat time stamps as dashed vertical lines 808. For the sake of clarity, only two lines 807 and only two vertical lines 808 are numbered. The method allows estimating the inter-beat-intervals and the heart rate variability of the cardiac rhythm based on the second set of periodicity feature estimates.

[0127] The sixth graph 806 shows the ECG of the same cardiac rhythm represented by the signal. The comparison of the fifth graph 805 and the sixth graph 806 shows that the inter-beat-intervals and the beats estimated according to the method of the present application match the ECG accurately.

[0128] In an embodiment, the method comprises determining the occurrence of a sleep disordered breathing (SDB) event based on the second set of periodicity feature estimates. The SDB event can be identified by a typical change in the cardiac rhythm.

[0129] In an embodiment, the method comprises determining the sleep stage of the subject 1001 based on the second set of periodicity feature estimates.

[0130] Figure 9 A method according to a second embodiment of the present application is described. The method is for processing a signal representing a physiological rhythm of a subject 1001. The physiological rhythm has periodicity features. The method comprises receiving the signal 100 in a receiving step 901. The method comprises estimating a first set of periodicity feature estimates based on the signal 100 in a step 902. The method comprises generating a prior probability for a periodicity feature based on at least a subset of the first set of periodicity feature estimates in a step 903. The method comprises estimating a second set of periodicity feature estimates based on the first set of periodicity feature estimates and the prior probability in a step 904.

[0131] In an optional step 905, the method comprises determining a spectral likelihood of at least a portion of the signal 100. The spectral likelihood represents a likelihood that the first set of periodicity feature estimates corresponding to the at least a portion of the signal 100 represent a periodicity feature of the physiological rhythm. The method comprises estimating the second set of periodicity feature estimates 100 based on the spectral likelihood of the signal.

[0132] In optional step 906, the method comprises defining a landmark in the circadian rhythm. The landmark is a phase in the period of the circadian rhythm. The method comprises determining a time likelihood of the signal 100 corresponding to the period. The time likelihood represents a likelihood of the signal 100 representing the landmark. The method comprises estimating the landmark based on the spectral likelihood, the time likelihood and the prior probability.

[0133] Figure 10 A system 1000 for processing a signal 100 representing a circadian rhythm of a subject 1001 1001 is described. The circadian rhythm has a period characteristic. The system 1000 comprises a sensor 1002 and a computer 1010. The computer 1010 is coupled to the sensor 1002 by a signal line 1008. The computer 1010 is configured to perform the method of any of the above described embodiments. The sensor 1002 is an accelerometer attached to the chest of the subject 1001 1001. The sensor 1002 is attached to the chest with a belt 1002. The subject 1001 1001 is lying in bed and is asleep. The computer 1010 has an input device 1012 connected to the sensor 1002 and receiving the signal 100 from the sensor 1002. The input device 1012 forwards the signal 100 to a processor 1014. The processor 1014 processes the signal 100 according to the method of any of the above described embodiments. The processor 1014 is connected to an output device 1016. The output device 1016 outputs an output signal about the period characteristic estimated by the method. For example, the output device 1016 produces a video signal or a Bluetooth signal or any other type of signal suitable for use in conveying information about the estimated period characteristic to a person.

[0134] The computer 1010 executes a computer program product comprising instructions which, when executed by the computer 1010, cause the computer 1010 to perform the method according to the present application. The instructions are provided to the computer 1010, for example, by a computer-readable storage medium containing the instructions. For example, the computer 1010 has a memory storing the computer program product. The processor 1014 is configured to receive the instructions of the computer program product from the memory.

[0135] Variations of the disclosed embodiments can be understood and implemented by those skilled in the art by a study of the drawings, the present disclosure, and the appended claims. In the claims, the word "comprising" does not exclude other elements or steps, and the indefinite articles "a" or "an" do not exclude a plurality. A single processor or other unit can fulfil the functions of several items recited in the claims. The mere fact that certain measures are recited in mutually different dependent claims does not indicate that a combination of these measures cannot be used to advantage. A computer program product can be stored / distributed on a suitable medium, such as an optical storage medium or a solid-state storage medium supplied together with or as part of other hardware, but can also be distributed in other forms, such as via the Internet or other wired or wireless telecommunication systems. Any reference signs in the claims should not be construed as limiting the scope.

Claims

1. A computer-implemented method for processing signals (100) representing the physiological rhythms of a subject (1001), said physiological rhythms having periodic characteristics (T, t). i The method includes: Receive the signal (100); Based on the signal (100), estimate the first set of periodic features (T, t). i )estimate; Based on the first set of periodic features (T, t) i At least a subset of the estimated values ​​are used to generate the periodic features (T, t) i The prior probability (P(θ)) of ). Based on the first set of periodic features (T, t) i Estimate the prior probability (P(θ)) and estimate the second set of periodic features (T, t). i )estimate; Determine at least a portion of the spectral likelihood (f) of the signal (100). S ); The spectral likelihood (f) mentioned above S ) represents the first set of periodic features (T, t) i Estimate the periodic characteristics (T, t) representing the physiological rhythm. i The likelihood of the first set of periodic features (T, t) i The estimate corresponds to at least a portion of the signal (100); The second set of periodic characteristics (T, t) is estimated. i The estimation is based on the spectral likelihood (f) of the signal (100). S ); Define a location point in the physiological rhythm, wherein the location point represents a phase in the cycle of the physiological rhythm; Determine the time likelihood (f) of the signal (100) corresponding to the period. T ); The time likelihood (f) mentioned above T The signal (100) represents the likelihood of the location point; Based on the aforementioned spectral likelihood (f) S The time likelihood (f) T The location point is estimated using the prior probability (P(θ)) and the prior probability (P(θ)).

2. The computer implementation method according to claim 1, comprising: The estimation of the localization points includes estimation using maximum a posteriori probability (MAP). The maximum a posteriori probability (MAP) estimation mentioned above is based on the spectral likelihood (f S The time likelihood (f) T ) and the prior probability (P(θ)).

3. The computer implementation method according to claim 1 or 2, wherein, Compared to the first set of periodic features (T, t) i Estimate the second set of periodic characteristics (T, t) i To estimate more accurately the periodic characteristics (T, t) of the physiological rhythm. i ).

4. The computer implementation method according to claim 1 or 2, wherein the periodic feature (T, t) i This includes the period length (T) and the timestamp (t). i At least one of them.

5. The computer implementation method according to claim 1 or 2, wherein the physiological rhythm is the spontaneous breathing rhythm of the subject (1001).

6. The computer implementation method according to claim 1 or 2, wherein the physiological rhythm is the heart rhythm of the subject (1001).

7. The computer implementation method according to claim 6, comprising: Based on the second set of periodic features (T, t) i Estimate, at least one of the heart rhythm's beat interval and heart rate variability.

8. The computer implementation method according to claim 5, wherein the signal (100) includes a signal from an accelerometer disposed on the chest of the subject (1001).

9. The computer-implemented method according to claim 1 or 2, comprising: Based on the second set of periodic features (T, t) i Estimate and determine the occurrence of sleep-disordered breathing (SDB) events.

10. The computer-implemented method according to claim 1 or 2, comprising: Based on the second set of periodic features (T, t) i The sleep stage of the subject (1001) was estimated and determined.

11. A computer program product comprising instructions that, when executed by a computer (1010), cause the computer (1010) to perform the computer implementation method according to any one of claims 1 to 10.

12. A computer-readable storage medium comprising instructions that, when executed by a computer (1010), cause the computer (1010) to perform the computer-implemented method according to any one of claims 1 to 10.

13. A system (1000) for processing signals (100) representing the physiological rhythms of a subject (1001), the physiological rhythms having periodic characteristics (T, t). i The system (1000) includes a sensor (1002) and a computer (1010), wherein the computer (1010) is coupled to the sensor (1002), and wherein the computer (1010) is configured to perform a computer implementation method according to any one of claims 1 to 10.

Images