Method and device for calculating at least one index of variability of the heart rate from an rr series

Abstract

The calculation method is performed by a processing unit on the basis of an initial RR series or of an RR series derived from the initial RR series. The processing unit is configured to detect local minima (Pji) and / or local maxima (Mji) in a sliding main window (F), to divide the main window (F) into n (n ≥ 3) sub-windows (Fj), to calculate a lower curve (Cmin) and / or an upper curve (Cmax), and to calculate at least one index (IV) of variability of the heart rate from the portion of the lower curve (Cmin) in the central window(s) (F2, ..., Fn-1), excluding the first sub-window (F1) and the last sub-window (Fn), and / or from the portion of the upper curve (Cmax) in the central window(s), excluding the first sub-window (F1) and the last sub-window (Fn).

Description

[0001] METHOD AND DEVICE FOR CALCULATING AT LEAST ONE HEART RATE VARIABILITY INDEX FROM AN RR SERIES

[0002] Technical field

[0003] The present invention relates to the calculation of at least one characteristic index of heart rate variability from samples of an RR series constructed from a cardiac signal. This variability index can in particular be used to quantify the activity of the Autonomic Nervous System (ANS)

[0004] Prior art

[0005] From a physiological point of view, the heart of a living being, isolated from any external influence, contracts automatically in a very regular way like a metronome, under the action of the sinus node which generates an independent nerve impulse, and causes a spontaneous contraction of the cardiac muscle. The heart is not, however, isolated, but is connected to the Autonomic Nervous System (ANS).

[0006] The Autonomic Nervous System (ANS) is a part of the nervous system that regulates certain physiological processes such as, for example, blood pressure, heart rate, body temperature or respiratory rate. It is composed of two branches: a "sympathetic" branch, accelerating, and a "parasympathetic" (or vagal) branch, inhibiting. Its functioning is "autonomous", that is, without conscious effort from the person. The ANS is connected to the organs and sensory receptors of the body from which it receives information about the internal and external environment of the subject. In response to this information, it stimulates physiological processes, through its sympathetic branch, or inhibits them, through its parasympathetic branch. This regulation is necessary to maintain the internal homeostasis of the body.

[0007] One of the important functions of the ANS is cardioregulation, that is, the regulation of heart rate and blood pressure to meet the oxygen needs of tissues, for example, during physical exertion, stress, pain, emotion or other physiological events. This is a constant, continuous regulation, ordered by the two branches of the ANS, sympathetic (stimulation) and vagal (inhibition).

[0008] This phenomenon is particularly visible in the heart rate, which is subject to these influences, acceleration and deceleration, causing variability in the instantaneous heart rate. However, outside of these influences, the heart rate is under the control of the "sinus node", a nerve node located in the tissues of the right atrium of the myocardium, acting as a biological clock at a very regular rhythm. This variability in the heart rate is therefore a direct reflection of the activity of the ANS.

[0009] Furthermore, heart rate is a physiological parameter that is simple to implement and can be measured non-invasively, for example using electrodes for measuring the electrocardiogram (ECG), a plethismographic pulse sensor (PPG) placed on the finger, ear, forehead or wrist depending on the use of the device, or any other sensor such as ultrasound, microphones, piezoelectrics, smartphone cameras, etc.

[0010] In many situations, it is useful to know the activity of the ANS. This is particularly true in medicine, during surgical interventions, in resuscitation or intensive care situations, or even in outpatient monitoring, but also in non-medical fields such as sport, physical activities in general or relaxation. In particular, knowledge of the activity of the ANS can be of great help in developing a diagnosis of many clinical situations. On this subject, we can refer for example to the following publication: Lacroix D, Logier R., Kacet S., Hazard JR, Dagano J. (1992): “Effects of consecutive administration of central and peripheral anticholinergic agents on respiratory sinus arrhythmia in normal subjects, J. of the Autonomie Nervous System”, Vol 39, pages 211-218.To study these fluctuations in heart rate, for example, since 1970, different techniques have been developed for spectral analysis of a signal which represents the evolution over time of the instantaneous heart rate (or frequency), and which is obtained after sampling an analog bioelectric signal, characteristic of the heart rate of a living being, and subsequently called a cardiac signal.

[0011] An example of a heart rate variability analysis method known to date is:

[0012] - to acquire a cardiac signal, by any invasive or non-invasive means [for example, and without limitation, acquisition of an electrocardiographic (ECG) signal by means of an electrocardiograph, or use of a blood pressure sensor connected to a catheter introduced into an artery, or use of an infrared pulse sensor],

[0013] - to construct from this signal a series called RR which is made up of a plurality of samples (RRi) representing the time intervals which separate two successive heartbeats,

[0014] - and to carry out a spectral analysis of the RR series.

[0015] More specifically, the spectral analysis of an RR series from a cardiac signal is usually carried out in two main stages.

[0016] In a first step, the spectral density curve of the RR series is calculated, for example between 0 and 2 Hz, using various known methods. The most commonly used method is to calculate the fast discrete Fourier transform of the RR series, in predefined time windows, weighted by means of a predefined weighting window. Depending on the implementation envisaged, this can be a rectangular weighting window, or for example a Kaiser, Hamming, Hanning or Bartlett weighting window. The calculation time windows can also be predefined and fixed, or it can be a calculation time window, of predetermined size, which is slid over time. For example, the Fourier transform is performed in a sliding time window of 256 seconds, applied to the RR series, and subjected to Kaiser weighting to limit the edge effects due to windowing.

[0017] In a second step, from the spectral density curve obtained at the end of the first step, the spectral powers (areas under the spectral density curve) are automatically calculated between predetermined frequency limits, which can be adjusted by a user. These spectral power calculations make it possible to obtain quantitative information, which is characteristic of the activity of the Autonomic Nervous System (ANS), and thus constitute a means of investigating and analyzing cardiac regulation by the ANS.

[0018] The above-mentioned spectral analysis method, however, has several drawbacks.

[0019] An RR series is non-stationary. Therefore, applying a Fourier transform to this type of non-stationary series gives results that are imprecise, or even erroneous, and that cannot be interpreted with certainty. The calculation of the spectral density curve by fast Fourier transform (or equivalent), is relatively heavy in terms of computing power and / or computing time, which currently makes this spectral analysis method unsuitable and difficult to implement in real time, particularly in embedded systems.

[0020] To achieve acceptable frequency resolution, the fast Fourier transform must also be calculated in relatively large time windows (e.g., 256s), which corresponds to a large number of samples in the RR series. As a result, this spectral analysis method is accompanied by a memory effect that delays the recognition of a change in the RR series.

[0021] Furthermore, a method for calculating a heart rate variability index has already been proposed, in particular in international patent application WO2006 / 032739, which can more particularly be used to quantify the activity of the ANS. This method makes it possible to overcome the aforementioned drawbacks inherent in prior art methods using spectral analysis and is based on calculating a variability index from local minima or local maxima of a resampled, filtered and normalized RR series, which are detected in a sliding main window.More particularly, in this method the main window is divided into several sub-windows, a curve is determined by linear interpolation connecting all the local minima (when the samples of the RR series are respectively equal to the time intervals ôti which separate two successive heartbeats) or connecting all the detected local maxima (when the samples of the RR series are respectively equal to the inverse 1 / ôti of the time intervals which separate two successive heartbeats), and heart rate variability indices are calculated from this curve.

[0022] In this type of process, the sliding of the main window is typically performed with an overlap, the sliding step of the main window having a duration less than the duration of the main window. This overlap avoids the loss of information and promotes the real-time calculation of the variability index.

[0023] By implementing this type of process, however, oscillations or untimely variations in the value of the variability index have been observed at certain times over more or less long periods. These untimely and temporary oscillations or variations do not correspond to physiological reality and constitute artifacts which can harm the reliability of this index.

[0024] Presentation of the invention

[0025] The main objective of the present invention is to propose a method for calculating at least one heart rate variability index from local minima and / or local maxima of an RR series which are detected in a sliding main window with overlap, which method has been improved so as to reduce, and preferably to eliminate, the aforementioned untimely and temporary oscillations or variations of the variability index which do not correspond to physiological reality.

[0026] Summary of the invention

[0027] The invention thus has as its first object a method for calculating at least one heart rate variability index, implemented automatically by a processing unit, from an initial RR series consisting of a plurality of samples (RRi) which are a function respectively of the time intervals (eti) which separate two successive heartbeats or from an RR series derived from the initial RR series, in particular a resampled and / or normalized and / or filtered RR series.

[0028] This method of the invention comprises a step of calculating (a) at least one heart rate variability index in a sliding main window (F) of duration (T), with a sliding step of at least one sample and of a duration (d) less than the duration (T) of the main window; the processing unit is configured to, during said calculating step (a), detect local minima (Pji) and / or local maxima (Mji) in the main window (F), to divide the main window (F) into n (n>3) sub-windows (Fj) consisting of a first sub-window (Fi), one or more central windows (F2, ... ,F n -i) and a last sub-window (F n ), to calculate a lower curve (Cmin) connecting at least the last local minimum of the first sub-window (F1), at least each local minimum (Pji) included in the central sub-window(s) (F2, ... , F n-i), and at least the first local minimum of the last sub-window (F n ) and / or an upper curve (Cmax) connecting at least the last local maximum of the first sub-window (F1), at least each local maximum (Mji) included in the central sub-window(s) (F2, ... ,F n -i), and at least the first local maximum of the last sub-window (F n ), and to calculate said at least one index (IV) of heart rate variability from the portion of the lower curve (Cmin) in the central window(s) (F2, ... ,F n -i), excluding the first sub-window (F1) and the last sub-window (F n ) and / or from the portion of the upper curve (Cmax) in the central window(s), excluding the first sub-window (Fi) and the last sub-window (F n ).

[0029] To achieve the invention, it was demonstrated that the aforementioned untimely and temporary oscillations or variations of the variability index could be linked to an uncertainty in the calculation of the aforementioned curve at the extreme edges of the main window, this uncertainty in the calculation being able to be amplified by the overlap in the sliding of the main window. In the invention, the exclusion of the first and last sub-windows makes it possible to overcome this drawback.

[0030] In this text, the term "cardiac signal" refers to any physical signal characteristic of the instantaneous heart rate of the living being. To implement the invention, various invasive or non-invasive techniques may be used to acquire this cardiac signal. A known invasive technique consists, for example, of using a blood pressure sensor connected to a catheter introduced into an artery. Among the known non-invasive methods (and those which will be preferably chosen), we can cite, for example, the use of an infrared pulse sensor, the use of an ultrasound sensor allowing the detection of cardiac cycles, of the type of sensor implemented in a cardiotocograph, or the acquisition of an electrocardiographic signal (ECG).The acquisition of an electrocardiographic signal (ECG) is in practice the most commonly used method, because in addition to its non-invasive nature, it allows obtaining a more precise cardiac signal than that obtained, for example, by means of an infrared pulse sensor.

[0031] In this text, the term "initial RR series" generally refers to a series of several successive samples (RRi), obtained after sampling an analog cardiac signal characteristic of the heart rate of the living being, each sample (RRi) being generally a function of the time interval (<5ti) between two successive heartbeats. More particularly, each sample (RRi) may be equal or proportional to the time interval (éti) between two successive heartbeats or may be equal or proportional to the inverse (1 / ôti) of this time interval.

[0032] In the preferred embodiment described below with reference to the appended figures, this initial RR series is more particularly constructed from the R waves of an ECG signal. This is not, however, limiting of the invention. In the case of an ECG-type cardiac signal, the so-called “RR” series can be constructed by using the other depolarization waves (P, Q, S or T) of the ECG signal to construct the RR series, the precision being, however, less good than by using the R waves of the ECG signal. Also, when the cardiac signal is not an ECG signal, the samples of the initial RR series are not calculated by determining the time interval separating two successive R waves of the ECG signal, but are more generally determined by detecting in the cardiac signal the time interval between two successive heartbeats.

[0033] The variability indices calculated using the method of the invention can advantageously be used to quantify the activity of the ANS or to characterize or detect any stimulus having an effect on the activity of the ANS, and resulting in a variation in the variability of the heart rate.

[0034] More particularly, the method of the invention may include the following additional and optional features, each additional and optional feature being able to be combined with any of the other technical features of the method:

[0035] - The processing unit is configured to determine, in each central sub-window (F2, ..., F n-i), an area (Aj) which is delimited at least by the lower curve (Cmin) and is above the lower curve (Cmin), and which extends over the entire duration (D) of the sub-window and / or an area (Aj) which is delimited at least by the upper curve (Cmax) and is under the upper curve (Cmax), and which extends over the entire duration (D) of the sub-window and to calculate said at least one index (IV) of variability of the heart rate from said areas (Aj). - The processing unit is configured to determine, in each central sub-window (F2, ... ,F n -i), an area (Aj) which is delimited by the lower curve (Cm in) and by the upper curve (Cmax) and which extends over the entire duration (D) of the sub-window, and to calculate said at least one index (IV) of variability of the heart rate from said areas (Aj)

[0036] - The main window (F) is preferably divided into at least four sub-windows (n>4) and the processing unit is configured to select from among said areas (Aj), the one with the lowest value, and to calculate said at least one index (IV) of heart rate variability from this minimum area (Smin).

[0037] - The processing unit is configured to calculate the ratio (Smin / S) between said minimum area (Smin) and the total area (S) of the main window (F), which is the height (H) of the main window F multiplied by the total duration (T) of the main window F and to calculate said at least one index (IV) of heart rate variability from this ratio (Smin / S).

[0038] - The duration (T) of the main window (F) is preferably between 32s and 256s, and more preferably is equal to 64s.

[0039] - The number (n) of sub-windows is preferably greater than or equal to 4.

[0040] - Sub-windows can have identical durations.

[0041] - The duration of the first sub-window (F1) is preferably equal to the duration of the last sub-window (F n ).

[0042] - the calculation step (a) is carried out on an RR series which is obtained at least by resampling the initial RR series or an RR series derived from the initial RR series, at a sampling frequency (Fe) greater than or equal to 4Hz, and preferably equal to 8Hz.

[0043] - The duration (d) of the sliding step is greater than or equal to the inverse (1 / Fe) of the sampling frequency (Fe), and preferably is equal to 1 s. - the calculation step (a) is carried out on an RR series which is obtained at least by normalizing the samples (RRi) of the initial RR series or the samples (RRi) of an RR series derived from the initial RR series, over the entire duration (T) of the main window (F).

[0044] - the calculation step (a) is carried out on an RR series which is obtained at least by filtering the initial RR series or an RR series derived from the initial RR series in at least one predefined frequency band [f1 - f2].

[0045] - Filtering allows filtering in at least one high frequency band [0.15Hz-0.4Hz],

[0046] - Filtering allows filtering in at least one low frequency band [0.04Hz; 0.15Hz],

[0047] - The processing unit is configured to display the heart rate variability index(es) (VI) in real time on a screen.

[0048] Another subject of the invention is a device for calculating at least one index (IV) of heart rate variability, comprising at least one processing unit configured to automatically implement the above-mentioned calculation method.

[0049] Preferably, the device further comprises means for acquiring a cardiac signal and means for sampling this analog cardiac signal, and the processing unit is configured to also construct said initial RR series from the signal sampled by the sampling means.

[0050] Preferably, the processing unit is configured to automatically implement the above-mentioned calculation method in real time, as the cardiac signal is acquired and the initial RR series is constructed.

[0051] Another subject of the invention is a computer program product comprising program code instructions and making it possible, when executed by at least one processing unit, to carry out the above-mentioned calculation method. Description of the figures

[0052] Other characteristics and advantages of the invention will appear more clearly on reading the detailed description which is given by way of non-limiting and non-exhaustive example, and with reference to the appended drawings in which:

[0053] - figure 1 schematically represents the main elements of an example of a device for acquiring and processing a cardiac signal, which can be configured to calculate indices of variability of the heart rate in accordance with the method of the invention,

[0054] - Figure 2 represents the set of waves (PQRST) characteristic of a heartbeat in an ECG signal,

[0055] - Figure 3 represents an example of a digital ECG signal, obtained after sampling an analog ECG signal,

[0056] - figure 4 represents an example of an initial RR series (also referred to as the initial RR signal) constructed from the signal of figure 3, and

[0057] - Figure 5 represents an example of area calculation on an RR series after normalization and filtering of the signal and after centering and amplitude reduction in a range of unitless values ​​between -0.1 and +0.1.

[0058] - Figure 6 shows another example of area calculation on an RR series after normalization and filtering of the signal and after centering and amplitude reduction in a unitless range of values ​​between -0.1 and +0.1.

[0059] Detailed description

[0060] Device for acquiring and processing a cardiac signal

[0061] Figure 1 shows an example of a device which allows the acquisition and processing of the cardiac signal of a living being (hereinafter referred to as “patient”) in accordance with the invention.

[0062] This system comprises: - standard means for acquiring an ECG signal, comprising several measuring electrodes 1 connected at the input to an electrocardiographic (ECG) monitor 2,

[0063] - means 3 for processing the ECG signal output by the ECG monitor 2.

[0064] The ECG signal processing means 3 comprise an analog / digital converter 4, and a processing unit 5. The input of the converter 4 is connected to the output of the ECG monitor 2, and the output of the converter 4 is connected to an input port of the processing unit 5.

[0065] The processing unit 5 can be implemented with various technologies and is preferably a programmable electronic processing unit, for example of the microcomputer type, programmable controller, electronic card with microprocessor(s) or microcontroller(s) or electronic card with programmable electronic circuit(s) of the FPGA type. The processing unit 5 can also be implemented with specific integrated electronic circuit(s) of the ASIC type.

[0066] In the particular example of embodiment illustrated in Figure 1, and not limiting the invention, the processing unit 5 is constituted by a microcomputer, the converter 4 being connected to a communication port of this microcomputer (for example RS232 serial port or USB port).

[0067] In other embodiments, the analog / digital converter 4 and the processing unit 5 can be integrated into an on-board system, for example of the “holster” type and capable of being worn by a patient.

[0068] In operation, the electrodes 1 are applied to the body of a patient, and the ECG monitor outputs in the usual manner an analog electrical signal, called an ECG signal, which for each heartbeat, for example, has the form of the signal shown in Figure 2.

[0069] Referring to Figure 2, for each heartbeat, this electrocardiographic signal (ECG) consists of a set of electrical waves:

[0070] - the P wave, which corresponds to the depolarization of the atria, and which has a low amplitude and a dome shape;

[0071] - the PQ space which reflects the atrioventricular conduction time;

[0072] - the R wave considered in practice as a marker of ventricular systole, or heartbeat, the QRS complex reflecting ventricular contraction, and the T wave which reflects ventricular repolarization.

[0073] This analog ECG signal is digitized by converter 4, with a predetermined sampling frequency (Fe), for example 256 Hz.

[0074] The sampled signal delivered at the output of the converter 4, (signal shown in FIG. 3) is processed by the processing unit 5, by means of specific calculation software which is described in detail later. This calculation software is stored in the memory of the processing unit 5, and allows, when executed, to automatically calculate, from the digital signal delivered by the analog / digital converter 4, one or more indices (IV) of variability of the heart rate.

[0075] A preferred variant of this calculation software will now be detailed.

[0076] Example of a software algorithm for calculating heart rate variability index(es)

[0077] In a preferred embodiment of the invention, the main successive steps of the algorithm of the software for calculating heart rate variability index(es) are as follows:

[0078] 1. Acquisition of RRi samples from the initial RR series and resampling

[0079] 2. Selection of RRk samples from the resampled RR series included in a main time window F of duration T [T (in seconds) >1 / f] 3. Normalization of the signal.

[0080] 4. Filtering

[0081] 5. Division of the main time window F into n [n >3] sub-windows Fj

[0082] 6. Detection of local minima Pji and local maxima Mji in the main window F

[0083] 7. Calculation of Cmin and Cmax curves

[0084] 8. Calculation of areas Aj

[0085] 9. Calculation of one or more IV indices of heart rate variability

[0086] 10. Sliding (with overlap) of the main window F and repeating the calculation from step 2

[0087] In practice, the system can be programmed to be used in real time or in delayed time.

[0088] When the system is used in deferred time, step 1 is first carried out in real time so as to acquire all the RRi samples over the entire desired analysis period; all of these successive RRi samples are stored in an acquisition file in the processing unit's memory. In a second step, steps 2 to 10 are carried out in a loop, in deferred time, on the RRi interval values ​​stored in the acquisition file.

[0089] When the system operates in real time, the acquisition step

[0090] 1 on the one hand, and the other processing steps 2 to 10 on the other hand are executed by two separate software modules operating in parallel, the first acquisition module (step 1) feeding the second processing and calculation module (steps 2 to 10) by means of a buffer file or equivalent.

[0091] Steps 1 to 10 will now be detailed.

[0092] Acquisition of RRi samples from the initial RR series and resampling

[0093] The acquisition of the RRi samples of the initial RR series is carried out by a first software sub-module which is supplied as input with the successive digital data constituting the digitized ECG signal (signal of figure 3) delivered by the analog-digital converter 4. Each data (or point) of the ECG signal is defined by the instantaneous amplitude ECGi of the ECG signal, and by the sampling instant (t, = ni / Fe, with ni sample number and Fe representing the sampling frequency of the converter 4).

[0094] The first R sample acquisition sub-module is designed to automatically detect each successive peak Ri in the digital signal delivered by the converter 4, and to automatically construct a series RR (figure 4) consisting of a succession of samples RRi. Each sample RRi is defined by the pair of coordinates: ti [a sampling instant (or number)]; time interval ôti (expressed for example as a multiple of the sampling period 1 / Fe) separating a peak Ri from the following peak Ri+i (in another variant it could be the previous peak RM).

[0095] Usually in itself, the R wave being most often the thinnest and widest part of the QRS, it is preferably used to detect the heartbeat with very good precision, the time interval oti corresponding in practice to the time separating two successive heartbeats. However, in another variant, one could consider using other waves (for example Q wave or S wave) of a heartbeat of the ECG signal to detect and construct the RR series.

[0096] The initial RR series (figure 4) delivered by the first sub-module mentioned above is resampled by a second software sub-module at a predefined frequency f. The objective of this resampling is to obtain as output an RR series whose samples RRk are equidistant from a temporal point of view, that is to say in other words an RR series whose sampling instants are equidistant temporally.

[0097] This resampling can be carried out in a manner known per se by interpolation, and for example by linear interpolation, in order to constitute a sequence of samples which are temporally equidistant.

[0098] The resampling frequency f is preferably greater than twice the patient's physiological maximum heart rate.

[0099] In order to take into account all possible physiological situations, the resampling frequency f is preferably set to a value greater than 4Hz, and is for example 8Hz.

[0100] Step 2: Selection of RRi samples included in a main time window of duration T [T (in seconds) >1 / f]

[0101] This step involves isolating a number N of successive RRi samples (N=Tf).

[0102] The duration T of the main window F is preferably greater than or equal to 32s. Below 32s, there is a risk that the heart rate variability index is not significant from a physiological point of view.

[0103] The duration T of the main window F is preferably less than or equal to 256s, in order to reduce calculation times and allow implementation in substantially “real time”.

[0104] As an indication, we choose for example and preferably a main window F having a duration T of 64 seconds, which corresponds to 512 successive RRk samples (N=512), for a resampling frequency f of 8Hz.

[0105] The following steps 3 to 9 are applied to the RRk samples of the resampled RR series included in this main window. Step 3: Signal Normalization

[0106] This step is carried out by means of a software sub-module which first calculates the S norm of the signal from step 2 in accordance with the following usual formula:

[0107] Where p is an integer greater than 1, and Si are the discrete values ​​of the signal. Alternatively, the norm can also be calculated using the following formula:

[0108] Preferably, for the implementation of the invention, p will be chosen equal to 2.

[0109] In a second step, the software sub-module performs a normalization of the signal by dividing each value Si of the signal by the previously calculated norm S.

[0110] This normalization makes it possible to free oneself from the influence of the average value of the heart rate.

[0111] It is essential that this normalization is performed across the entire width of the main window, i.e., taking into account the (N) RRk samples of the main window in the normalization. In contrast, a normalization performed on each sub-window (Fj) would be completely ineffective.

[0112] Preferably, during this step, the software sub-module also performs centering and reduction of the amplitude of the RRk samples in a unitless range, for example between -0.1 and +0.1, such as in the example of Figure 5. This operation makes it possible to overcome the influence of inter-individual amplitude variations.

[0113] Step 4: Filtering

[0114] This step consists of applying a band-pass filter to the RRk samples of the resampled and normalized RR series included in the main window so as to keep only the frequencies included in a predefined frequency band [f1; f2].

[0115] More particularly, the frequency band [f1; f2] is equal to or is included in the frequency band [0.04Hz; 5Hz],

[0116] To carry out this bandpass filtering step, we will use, for example, a high-pass digital filter having a cutoff frequency at frequency f 1 , in series with a low-pass digital filter having a cutoff frequency at frequency f2. It is also possible to use a recursive selective filter with infinite impulse response (IIR filter) centered on the frequency band [f1; f2].

[0117] The high-pass filter (cutoff frequency f1) aims to filter very low frequencies, and incidentally to remove many artifacts in the signal. In practice, the cutoff frequency f1 is greater than or equal to 0.04 Hz.

[0118] The low-pass filter (cut-off frequency f2) aims to filter high frequencies, typically above 5 Hz and preferably above 1 Hz, because in practice they do not contain any interesting information. The implementation of this low-pass filtering, although preferential, is however optional and not essential to the realization of the invention.

[0119] In a preferred embodiment, the frequency band [f1; f2] is equal to [0.15Hz; 0.4Hz] in order to retain only the high frequency variations of the RR series (respiratory variations) relating to the parasympathetic tone.

[0120] In another preferred embodiment, the frequency band [f1; f2] may be equal to [0.04Hz; 0.15Hz]. This low frequency range is known to be influenced by both the sympathetic system and the parasympathetic system.

[0121] Step 5: Splitting the main time window F into n [n >3] sub-windows Fj

[0122] The software sub-module corresponding to step 5 automatically divides the main time window F into n [n >3] sub-windows Fj.

[0123] In an alternative embodiment, the duration of the first sub-window Fi is preferably (but not necessarily) equal to the duration of the last sub-window F n .

[0124] In a simplified embodiment, all the sub-windows Fj can have identical durations.

[0125] As an indication and non-limiting example of the invention, in the particular example of figure 5, the main window has been divided into four (n=4) sub-windows Fi, F2, F3 and F4 of identical durations D (duration D of each sub-window equal to 16s for a main window F of 64s), namely a first sub-window F1, two central sub-windows F2 and F3 and a last sub-window F4.

[0126] In this particular example in Figure 5, the two windows F1 and F4 are the two extreme windows that are excluded and are not considered for the calculation of the heart rate variability indices IV in step 9.

[0127] Step 6: Detection of local minima Pji and local maxima Mji in the main window

[0128] Step 6 is performed by means of a software sub-module which detects local minima Pji and local maxima Mji of this signal.

[0129] This detection is preferably carried out on the principle of measuring the maximum amplitude between two successive zero crossings of the RR series. Thus, a local minimum is Pji or a local maximum Mji is detected as the point having the maximum amplitude relative to zero between two successive zero crossings of the RR series.

[0130] Detection can also be achieved by means of an algorithm detecting a slope reversal (or equivalently a change in the sign of the signal derivative).

[0131] The detection of local minima and local maxima can also be optimized to retain only the Pji minima whose value is lower than a predefined or calculated threshold and to retain only the Mji maxima whose value is higher than a predefined or calculated threshold, in order not to retain local minima or local maxima that are not significant and linked to non-physiological rebounds.

[0132] In the example of Figure 5, the local minima that were detected are as follows: - First sub-window Fi: Pu, P12, P13

[0133] - Central sub-window F2: P21, P22

[0134] - Central sub-window F3: P31, P32, P33

[0135] - Last sub-window F4: P41, P42

[0136] The local maxima that were detected are as follows:

[0137] - First sub-window F1: Mu, M12, M13

[0138] - Central sub-window F2: M21, M22, M23

[0139] - Central sub-window F3: M31, M32

[0140] - Last sub-window F4: M41, M42, M43

[0141] Step 7: Calculation of the Cm in and Cmax curves

[0142] In step 7, processing unit 5 automatically calculates:

[0143] - the lower curve Cmin which connects the local minima Pji detected in the previous step and

[0144] - the upper curve Cmax which connects the local maxima Mji detected in the previous step

[0145] These Cmin and Cmax curves are determined by interpolation, and preferably by linear interpolation, that is to say by calculating the straight line segments connecting the local minima (Cmin curve) and the straight line segments connecting the local maxima.

[0146] These Cmin and Cmax curves can also be determined by any other type of mathematical interpolation, including cosine, polynomial or polynomial by parts interpolations.

[0147] Step 8: Calculation of the areas Aj

[0148] Processing unit 5 automatically calculates for each central sub-window F2 to F n -i [i.e., windows F2 and F3 in the example in Figure 5], excluding the two extreme subwindows [first subwindow F1 and last subwindow F n, that is to say excluding windows F1 and F4 in the particular example of figure 5 with four sub-windows), the area Aj which extends, over the entire duration D of the central sub-window between the lower curve portion Cmin and the upper curve portion Cmax in the central sub-window. Calculation of one or more indices IV of heart rate variability

[0149] The processing unit 5 automatically calculates at least one heart rate variability index IV from the areas Aj calculated in the previous step.

[0150] In a preferred embodiment, the processing unit 5 automatically calculates at least one heart rate variability index IV from the smallest (Smin) of the areas calculated in the previous step, i.e. in the particular example of Figure 5, from the smallest of the areas A2 and A3 [Smin= min (A2; A3)]. In this particular example of Figure 5, the minimum area Smin is equal to A2.

[0151] More particularly and in a non-limiting manner of the invention, the index IV is a function of the minimum area Smin, and for example is calculated by means of an affine function with Smin as a variable, and can more particularly be calculated by means of the following formula:

[0152] IV =100.(a.Smin +b) / S with:

[0153] - S representing the total surface area of ​​the main window F

[0154] - a and b are coefficients that can be determined empirically.

[0155] In the particular example of figure 5 (main window F of duration T equal to 64 seconds and height H equal to 0.2 (between -0.1 and +0.1), the total surface area S (S=HT) of the main window F is equal to 12.8.

[0156] This IV index advantageously allows quantification of the activity of parasympathetic tone.

[0157] Preferably, the value of this index IV can be visualized by being displayed on a screen by the processing unit 5 Sliding (with overlap) of the main window F and repeating the calculation from step 2.

[0158] The processing unit 5 shifts the main time window of duration n seconds, by a time step equal to d seconds less than the duration T of the main window F (d < T), and loops back to step 2.

[0159] The step d for sliding the main calculation window influences the sensitivity of the index.

[0160] As an indication, and in practice, very satisfactory results have been obtained by choosing a step d worth 1 s.

[0161] Other variants

[0162] In another embodiment variant, the aforementioned filtering step (Step 4) could be carried out before the selection step 2, for example continuously as the RRi samples of the initial RR series are acquired. Also, in another embodiment variant, this filtering step could be carried out before the normalization step (abovementioned step 3).

[0163] In another variant, an index IV allowing the quantification of the activity of the parasympathetic tone can also be calculated from only the curve Cmin connecting the local minima Pji, for example by calculating for each central sub-window an area Aj, which is delimited at least by the lower curve Cmin and is above the lower curve Cmin, and which extends over the entire duration D of the central sub-window. This is for example an area Aj which extends between this curve Cmin and a predefined horizontal line, which can more particularly but not necessarily be the abscissa axis, as illustrated for example in Figure 6.

[0164] In another variant, an index IV can also be calculated from only the Cmax curve connecting the local maxima Mji, for example by calculating for each central sub-window an area Aj, which is delimited at least by the upper Cmax curve and is under the upper Cmax curve, and which extends over the entire duration D of the central sub-window. This is for example an area Aj which extends between this Cmax curve and a predefined horizontal line, which can more particularly but not necessarily be the abscissa axis. In another variant, the Cm in curve can be calculated by taking into account only:

[0165] - the last local minimum of the first sub-window (i.e. for figure 5, point P13 of sub-window F1), the other local minima of the first sub-window not being taken into account,

[0166] - the first local minimum of the last sub-window (i.e. for figure 5, point P41 of sub-window F4), the other local minima of the last sub-window not being taken into account,

[0167] - the local minima of the other central sub-windows (i.e. sub-windows F2 and F3 for figure 5).

[0168] In another variant, the Cmax curve can be calculated taking into account only:

[0169] - the last local maximum of the first sub-window (i.e. for figure 5, point M13 of sub-window F1), the other local maxima of the first sub-window not being taken into account,

[0170] - the first local maximum of the last sub-window (i.e. for figure 5, point M41 of sub-window F4), the other local minima of the last sub-window not being taken into account,

[0171] - the local maxima of the other central sub-windows (i.e. sub-windows F2 and F3 for figure 5).

[0172] In another variant, a heart rate variability index IV can be calculated in the same way as the AUCmoy parameter or the AUCmax parameter described in international patent application WO2006 / 032739, but excluding in the calculation the areas calculated in the first sub-window and in the last sub-window of the main window.

[0173] In the exemplary embodiments which have been described with reference to the appended figures, the calculation algorithms are applied to an initial RR series whose samples (RRi) are equal or proportional respectively to the time intervals (éti) which separate two successive heartbeats. In another variant embodiment, the invention can be implemented from an initial RR series whose samples (RRi) are equal or proportional respectively to the inverse (1 / ôti) of the time intervals between two successive heartbeats. In this case, by applying the aforementioned steps 1 to 7, curves are obtained which are the inverse of those of Figure 5.

[0174] The above-mentioned IV variability indices can generally be used to measure the variability of the heart rate in a living being (conscious or unconscious), and therefore to study the effects on the ANS (Autonomic Nervous System) of any stimulus likely to modify the activity of the ANS.

[0175] As non-limiting and non-exhaustive examples of the invention, the IV indices can be used in the medical or surgical field for: the evaluation of the pain felt by a conscious living being that is not mechanically ventilated (i.e. a living being whose respiratory frequency is arbitrary and variable, and is not imposed by a controlled ventilation device in contrast in particular with a patient under general anesthesia), for example a conscious living being under local anesthesia of the epidural type,

[0176] - the assessment of pain felt by an unconscious living being, and in particular a living being under general anesthesia; in the latter case, the assessment of the level of pain makes it possible to indirectly know the level of analgesia during general anesthesia.

[0177] As non-limiting and non-exhaustive examples of the invention, the IV indices can also be used, outside the medical or surgical field, for example for the evaluation of stress felt by a living being.

Claims

CLAIMS 1. Method for calculating at least one index (IV) of heart rate variability, implemented automatically by a processing unit (5), from an initial RR series consisting of a plurality of samples (RRi) which are a function respectively of the time intervals (eti) which separate two successive heartbeats or from an RR series derived from the initial RR series, which method comprises a step of calculating (a) at least one index of heart rate variability in a sliding main window (F) of duration (T), with a sliding step of at least one sample and of a duration (d) less than the duration (T) of the main window, and in which the processing unit (5) is configured to, during said calculation step (a), detect local minima (Pji) and / or local maxima (Mji) in the main window (F),to divide the main window (F) into n (n>3) sub-windows (Fj) consisting of a first sub-window (Fi), one or more central windows (F2, ... ,F, n -i) and a last sub-window (F n ), to calculate a lower curve (Cmin) connecting at least the last local minimum of the first sub-window (F1), at least each local minimum (Pji) included in the central sub-window(s) (F2, ... ,F n -i), and at least the first local minimum of the last sub-window (F n ) and / or an upper curve (Cmax) connecting at least the last local maximum of the first sub-window (F1), at least each local maximum (Mji) included in the central sub-window(s) (F2, ... ,F n -i), and at least the first local maximum of the last sub-window (F n), and to calculate said at least one index (IV) of heart rate variability from the portion of the lower curve (Cmin) in the central window(s) (F2, ..., F n -i ), excluding the first sub-window (F1) and the last sub-window (F n ) and / or from the portion of the upper curve (Cmax) in the central window(s), excluding the first sub-window (Fi) and the last sub-window (Fn).

2. Method according to claim 1, in which the processing unit (5) is configured to determine, in each central sub-window (F2, ..., F n-i), an area (Aj) which is delimited at least by the lower curve (Cmin) and is above the lower curve (Cm in), and which extends over the entire duration (D) of the sub-window and / or an area (Aj) which is delimited at least by the upper curve (Cmax) and is under the upper curve (Cmax), and which extends over the entire duration (D) of the sub-window and to calculate said at least one index (IV) of heart rate variability from said areas (Aj).

3. Method according to claim 1, in which the processing unit (5) is configured to determine, in each central sub-window (F2, ..., F n -i), an area (Aj) which is delimited by the lower curve (Cmin) and by the upper curve (Cmax) and which extends over the entire duration (D) of the sub-window, and to calculate said at least one index (IV) of variability of the heart rate from said areas (Aj).

4. Method according to claim 2 or 3, in which the main window (F) is divided into at least four sub-windows (n>4) and in which the processing unit (5) is configured to select from among said areas (Aj), the one whose value is the lowest, and to calculate said at least one index (IV) of variability of the heart rate from this minimum area (Smin).

5. Method according to claim 4, wherein the processing unit (5) is configured to calculate the ratio (Smin / S) between said minimum area (Smin) and the total surface area (S) of the main window (F) and to calculate said at least one index (IV) of heart rate variability from this ratio (Smin / S).

6. Method according to any one of claims 1 to 5, in which the duration (T) of the main window (F) is between 32s and 256s, and preferably is equal to 64s.

7. Method according to any one of claims 1 to 6, in which the number (n) of sub-windows is greater than or equal to 4.

8. Method according to any one of claims 1 to 7, in which the sub-windows have identical durations.

9. Method according to any one of claims 1 to 8, in which the duration of the first sub-window (Fi) is equal to the duration of the last sub-window (F n ).

10. Method according to any one of claims 1 to 9, in which the calculation step (a) is carried out on an RR series which is obtained at least by resampling the initial RR series or an RR series derived from the initial RR series, at a sampling frequency (Fe) greater than or equal to 4Hz, and preferably equal to 8Hz.

11. Method according to claim 10, in which the duration (d) of the sliding step is greater than or equal to the inverse (1 / Fe) of the sampling frequency (Fe), and preferably is equal to 1 s.

12. Method according to any one of claims 1 to 11, in which the calculation step (a) is carried out on an RR series which is obtained at least by normalizing the samples (RRi) of the initial RR series or the samples (RRi) of an RR series derived from the initial RR series, over the entire duration (T) of the main window (F).

13. Method according to any one of claims 1 to 12, wherein the calculation step (a) is carried out on a series (RR) which is obtained at least by filtering the initial RR series or an RR series derived from the initial RR series, in at least one predefined frequency band [f1 - f2].

14. The method of claim 13, wherein the filtering allows filtering in at least one high frequency band [0.15Hz- 0.4Hz], 15. Method according to claim 13 or 14, in which the filtering allows filtering in at least one low frequency band [0.04Hz; 0.15Hz], 16. Method according to any one of claims 1 to 15, wherein the processing unit (5) is configured to display the heart rate variability index(es) (IV) in real time on a screen.

17. Device for calculating at least one index (IV) of heart rate variability, comprising at least one processing unit (5) configured to automatically implement the calculation method of any one of claims 1 to 16.

18. Device according to claim 17, further comprising means (1, 2) for acquiring an analog cardiac signal and means (4) sampling this analog cardiac signal, and in which the processing unit (5) is configured to also construct said initial RR series from the signal sampled by the sampling means (4).

19. Device according to claim 18, wherein the processing unit (5) is configured to automatically implement the calculation method of any one of claims 1 to 16 in real time, as the cardiac signal is acquired and the initial RR series is constructed.

20. Computer program product comprising program code instructions and allowing, when executed by at least one processing unit, to carry out the calculation method of any one of claims 1 to 16.

Images