JP2013111467A - Method and device for measuring rsa component from heart rate data - Google Patents

Abstract

PROBLEM TO BE SOLVED: To provide a method and a system that estimate RSA component strength from human subject's heart rate data.SOLUTION: The method includes: sampling a heart rate by a transducer device such as a photoplethysmograph (PPG) in digital processing to obtain a time series of inter beat interval time; performing autoregressive (AR) analysis to windowed segments of data; and identifying a filter pole for correlating the strength with the present RSA component in each data window. The method specifies the calculation of a robust measurement value of the RSA component by formula in the figure, where, R is a component output measurement value, Eis the total window energy, k is a proportionality constant, and Pis the filter pole with the strongest correlation with the RSA component.

Description

  The present invention relates to a method for measuring the respiratory sinus arrhythmia (RSA) component of heart rate variability (HRV) and a system utilizing the method. The present invention further relates to biofeedback or anxiety management applications based on methods and systems for measuring the RSA component of heart rate variability.

  Biofeedback based on RSA has been shown to be useful in controlling anxiety and stress, and clinical trials have shown its effectiveness in reducing hypertension. Instant RSA can be used as a biofeedback parameter to modulate the audiovisual display to the subject or user to indicate how a certain pace of breathing exercise affects the person's heart pattern Can do. Simple examples of this include modulating the screen color displayed to the user and modulating the frequency of the dial tone emitted to the subject.

A typical measurement, RSA intensity S ′ (t), is the following equation (1):
Can be defined as the maximum value of the cross-correlation function of the functions of heart rate signal H (t) and respiration rate B (t) calculated over a synchronized time window, where * is Represents the cross-correlation function applied to the two signals over the same time interval, and max represents the maximum value function.

  Since the respiration rate is less convenient to measure than the heart rate, the measure S ′ (t) is defined a priori as a sinusoidal function of B (t) with an amplitude of 1.0 so that the user breathes. In many cases, the target respiration frequency to be instructed is estimated after being defined as a range of 0.05 to 0.2 hertz. This allows S ′ (t) to be estimated without directly measuring expiratory airflow or lung volume.

  Several sensor technologies can be used to convert cardiac activity into electrical signals that can be used for processing and RSA feature analysis. These sensors include a microphone (voice heart signal), a pressure sensor (pulse pressure), an electrocardiograph (ECG), a photoelectric pulse wave meter (PPG), and a non-contact sensor utilizing RF technology.

  Current methods for determining RSA are based on cross-correlation or power spectrum analysis, which are often too computationally intensive to run in real time on many mobile devices.

  From the above, in the market, there is a need for a method that can be used in a wide range of computing devices, including mobile devices with limited resources, in the method of extracting the RSA component of the heart rate signal.

  The present invention discloses an RSA biofeedback application that guides the user about techniques for breathing that are paced so that a state of cardiac coherence is achieved. The present invention further discloses a scoring system that rewards the subject by assigning points when the goal is met.

  The present invention provides a method for efficiently estimating the metric of the RSA component of the heart rate of a human subject from a raw beat interval data stream.

The RSA intensity estimation method divides the digital heart rate signal into overlapping windowed data segments and calculates the best-fit all-pole linear filter (autoregressive model) for each window of the beat interval. A single pole in the target sector of the complex plane is identified as the RSA pole, and the RSA intensity estimate is the following equation:
Where | P _RSA | is the absolute value of the identified RSA pole, E _w is the total window energy, and k is a proportionality constant.

  The main advantage of this method is by scaling the frequency of the filter model, which allows the method to scale its processing efficiency. This is particularly useful when trying to offload processing to a computing device such as a mobile device that has insufficient computing power to perform complex Fast Fourier Transform (FFT) processing in real time. . The autoregressive methodology further allows for higher resolution RSA identification using shorter data segments. A shorter data window may be preferred for certain biofeedback applications because it can speed up the update of the subject's physiological state.

  These embodiments and other aspects of the invention will be readily apparent from the following detailed description and the accompanying drawings, which are intended to illustrate the invention by way of example. It is not intended to limit the invention.

1 is a system diagram of an illustrative embodiment of the invention. FIG. FIG. 4 is a flow diagram of one preferred method for extracting RSA metrics from raw heart rate sensor data. An Argan diagram showing the pole positions of a well-matched all-pole linear filter is drawn, with the identified RSA pole A and the phase angle B of this pole displayed.

  The invention will be more fully understood when the following detailed description is read in conjunction with the accompanying drawings. Although detailed embodiments of the present invention are disclosed herein, the disclosed embodiments are merely illustrative of the present invention, and the present invention may be embodied in various forms. I want you to understand that Accordingly, specific functional details disclosed herein are not to be construed as limiting, but merely as a basis for the claims and various aspects of the invention. It should be construed as a representative basis for teaching one skilled in the art in an embodiment that is substantially and adequately detailed to employ.

  FIG. 1 is a system diagram of an exemplary embodiment 100 of the present invention. The sensor 104 measures the heart rate of the human subject 102. Most preferably, the sensor 104 communicates raw heart rate sensor data to the computing device 108 via the communication link 106, which computes the heart rate data and extracts RSA metrics. In certain preferred embodiments, the communication link 106 is a short-range communication link, for example, through the use of Bluetooth communication. However, the communication link 106 is not limited to near field communication, for example, a wide area network such as the Internet, a wireless communication, or a local area communication network. In general, communication link 106 may be provided by any type of communication link or network that supports data communication. In a preferred embodiment, the RSA metric is used by a biometric application that is running on the computing device 108 and provides biofeedback to the subject 102 based on the processed sensor data. Although sensor 104 is depicted as a commercially available oximeter, the invention is not so limited and other heart rate sensors may be used. In addition, although the computing device 108 is depicted as a smartphone, it includes any mobile phone, laptop computer, pad computer, game console, set top box such as used in cable TV systems, and personal computers. A type of computing device may be used. The computing device 108 has at least a processor, a data storage for storing program codes and data, a communication performance for communicating with the sensor 104, and a display for displaying a graphic expression result by the biometric application. I have it.

Method for Computing RSA Intensity from Heart Rate Signals FIG. 2 is a flow diagram of one preferred method 200 for extracting RSA metrics from raw heart rate sensor data. In the method 200, a raw sensor data feed 202 comes from a heart rate measuring device, or sensor embodiment 104. Examples of such heart rate measurement measuring devices include photoelectric pulse wave (PPG) devices, which are optically acquired pulse wave diagrams, organ volume measuring devices, and electrocardiogram (ECG) devices. In general, the method 200 operates on some data signal from a sensor or transducer that measures a human beat or provides an electrical signal relating to aspects of the human beat. In one successful prototype, the sensor 104 is implemented as a NONIN 9560 oximeter connected to a Nokia N97 mobile device using a Bluetooth short range communication link. The NONIN Bluetooth oximeter data stream is accessed using the JSR82 JAVA Bluetooth software module available on the NOKIA N97 side, and PPG heart rate sensor data expressed as 16-bit unsigned data is sent to the NOKIA N97 at a rate of 75 samples per second. Sent by.

  The raw sampled beat sensor data stream is converted to a time series of time intervals between beats, referred to herein as IBI data, by a first low pass filtering to reduce noise in step 204, and then in step 206. A peak detection algorithm determines the local maximum peak for each beat and outputs the time interval between peaks. Incoming raw data is filtered using a low pass filter and the trend is removed by applying a moving average filter. In one embodiment, an order 4 Chebyshev filter is used for this purpose with a cutoff frequency of 4.0 Hz.

  This allows the zero average signal with reduced noise to pass through the threshold based on the peak detection algorithm in step 206. If the first derivative is equal to zero, it may indicate a peak, but by comparing all points where such first derivative is equal to zero, the peak is detected and each beat The maximum value within the positive range of the cycle is found. This is a known method in the field of digital signal processing (DSP) for robust peak detection on physiological signals.

At step 208, the IBI data in seconds is then analyzed to identify and replace outliers, where the outliers are out of the current normal fluctuations and trends in the IBI data stream. It is a value. Outliers can also occur due to measurement noise, errors in the peak detection algorithm, and even physiological causes such as ectopic beats that generate either very short or very long beats. Ectopic beats also occur in normal healthy subjects, and heart disease is suspected if the ratio of these beats is very high. Here, the abnormal value removal is realized by calculating the expected value of each beat. The first pass for handling outliers uses a simple moving average filter of the IBI data to calculate the average expected value, and then:
As shown in FIG. 4, the value outside the range centered on the expected value is defined as an abnormal value. The abnormal value is then replaced with the expected value in the IBI stream, and a signal with no error is created for later processing.

  In step 212, the expected IBI value is calculated using the AR filter model determined in step 218. Step 218 uses a method known as linear predictive coding (LPC) in speech signal processing, which is a robust measure of expected values in stable and metastable time series. give. The AR model calculation is returned from step 218. Applying the AR model to the previous nth order IBI value, where n is the model order, calculate an estimate of the next IBI value. AR model filter orders between 4 and 15 are used. Step 218 will result in the value, AR model value, corresponding to the previous IBI data window since step 218 occurs after step 212. This is valid based on the assumption that the IBI signal is quasi-stationary and that the new data window is taken at a relatively short interval. This assumption is valid when the window overlap is 5 seconds or less.

  Once the outliers have been removed, at step 210, the IBI data stream is resampled at a uniform rate. This is necessary because the raw data is sampled at each peak in the raw data stream and these peaks are not evenly spaced in time. In order to take advantage of linear system analysis, the data must be provided at a uniform rate, and thus the data must be processed to achieve this before subsequent processing can be applied. Uniform sampling is achieved by taking complementary points in the IBI data at evenly spaced points. Cubic spline interpolation is applied to the non-uniform IBI data, and re-sampling points are established at 0.25 s (4 Hz) intervals.

  In step 214, the uniform IBI data stream is windowed into segments of data of length n, each window overlapping m samples by the previous window. A preferred value for n is 128 and a preferred value for m is 32. Windowing refers to applying a function whose value is zero outside a certain interval. In order to reduce the distortion of the AR parameter caused by the windowing process, a Hamming window function is further applied to the window data.

  At step 216, the windowed data is then averaged to remove the DC offset of the data and detrended to remove any inherent trends from the data over the length of the data window. Is called. Mean removal and detrending of data in the window calculates the best-fit line through each data point in the window using a least mean square algorithm and calculates this centerline from each data point in the window. Done by subtracting.

At step 218, the data window is then processed by an autoregressive (AR) method to calculate the best-fit all-pole filter model of the data. The AR method is in the time domain and has the following equation:
Can be formulated as a linear prediction method used to find the best estimate of the next value based on the linear weighted sum of the previous values in the time series. here
Is the expected value,
Is a prediction coefficient. This is the following equation:
As can be seen, the input signal can be viewed as a transfer function in z-space that relates the output signal.

Equation 5 describes an all-pole filter, and each filter pole is represented by equation (6) below:
Can be described in polar coordinates as having the magnitude and phase angle shown in FIG.

The larger the size r of the pole P _i, the stronger the effect it has on the overall dynamics of the system. For each pole, the phase angle θ is related to the frequency at which this pole exerts its effect. Each pole can be mapped to a component of the data stream being modeled, and the intensity and frequency of this component can be tracked by the magnitude and phase shift of the pole with which it is matched.

There are several ways to calculate the parameters a ₁ to a _{p of the} all-pole model from the window of input data. In a preferred embodiment, the Burg method, also known as the maximum entropy method or MEM, is used. The Burg method is well known in the art and estimates the parameters of an all-pole filter model of order p based on a window of input data. The computational efficiency of the Burg method is advantageous over other methods, especially for low order systems. The Burg method is also particularly advantageous for the analysis of short data segments (eg, 30 seconds) that are more difficult to accurately estimate poles. In one embodiment, the filter order is adapted to the computing power of the computing device that performs the method. In one embodiment used with the aforementioned N97 instrument, a filter order of 9 has been used successfully.

The output of the Burg algorithm is a parameter set a ₁ , a ₂ ,... A _p that describes a polynomial that matches the transfer function, known as Yule-Walker polynomial. At step 220, the poles of the transfer function are calculated by solving the roots of this polynomial. In this step, any standard root solving algorithm may be used. In one embodiment, the Jenkins-Traub method, an efficient algorithm that gives robust results, is used.

Once the transfer function poles are calculated, in step 222 the poles most closely related to the RSA component of the beat data are identified. This pole is in a normal respiratory rate constrained sector of 4 to 30 breaths per second. This range is further constrained by the fact that if the breathing exceeds 20 breaths, the sympathetic nervous system reduces the effects of RSA, resulting in a decrease in the pole size. This is related to the frequency range from 0.07 Hz to 0.33 Hz. (A slow controlled respiration rate on the order of 5-7 breaths per minute has been shown to cause a peak in RSA component intensity.) Given a sampling rate f _s , this is the following equation (7 ), That is,
Corresponds to the sector angle range given by. The poles in this sector are compared in terms of size. The metric of the influence of each pole P _i can be defined as the pole energy E _i , where
It is.

The pole with the largest E _i in the target sector is identified. If there are other poles in the lower phase angle sector with pole energies greater than 0.7 × E _i max, the lowest of these pole angles is identified as the RSA pole.

The total window energy E _w of the averaged and detrended IBI data window is the following equation:
Calculated by

In step 224, the total intensity of the RSA component is calculated by the following equation (10):
Calculated using Where E _RSA is the pole energy associated with the selected RSA pole P _RSA , E _w is the total window energy, and k is a proportionality constant used to normalize the output to the target range. In one embodiment, the value of k is defined as the reciprocal of the number n of samples in the IBI data window, where n = 128 means a value of k = 0.008. Combining the above equation (9) and equation (10), the following equation:
Is given. The intensity S _RSA is a robust measurement of the current RSA component of the heart rate data. In the current embodiment, this value of S _RSA measurement exceeded 1.0 when RSA was strong. An S _RSA below 0.2 was an indication that there was little or no RSA component in the IBI data window.

In step 226, S _RSA , the total strength of the RSA component, is transmitted to computing device 308. In other embodiments, the S _RSA value is stored and sent to the computing device 108 on demand.

FIG. 3 depicts an Argan diagram showing the pole position of the best-fit all-pole filter, with the identified RSA pole A and the phase angle B of this pole displayed. The embodiment shown in FIG. 3 shows the pole positions calculated using a 9th-order all-pole filter model at a sampling rate of 2 Hz for one windowed segment of heart rate data. Within the sector of interest, one pole (labeled A) has the largest magnitude, is above the threshold intensity value, and is identified as the RSA pole. The intensity value S _RSA calculated from this pole position is the output RSA measurement value transmitted as step 426 in FIG.

100 System 102 Human Subject 104 Sensor 106 Communication Link 108 Computing Device

Claims (7)

In a computer-implemented method for calculating RSA from beat sensor data,
(A) sampling the subject's beat during a window of time with a sensor to obtain beat sensor data;
(B) transmitting, by the sensor, pulsation sensor data corresponding to each sample during the window period to a computing device;
(C) calculating an interval between beats from the beat sensor data by the computing device;
(D) obtaining an all-pole filter model that best fits the inter-beat interval data by the computing device;
(E) calculating the position of the filter model pole by the computing device;
(F) identifying, by the computing device, the pole having the largest absolute value in a predefined sector of the unit circle;
(G) identifying an RSA value by the computing device based on the identified poles. The step of identifying includes determining the component intensity represented as S _RSA by the equation:
The method of claim 1 wherein: P _RSA is the identified pole, E _w is the total window energy, and k is a proportionality constant.   The method of claim 1, wherein the sensor is a transducer device selected from the group consisting of an electrocardiograph, a photoelectric pulse meter, a pressure sensor, and a microphone.   The method of claim 1, wherein the calculated inter-beat interval is resampled at a uniform rate.   The method according to claim 1, further comprising removing the average value and detrending the sample in the window by the computing device.   The method of claim 1, wherein a Burg method is used to calculate the best-fit all-pole model filter coefficients. The pre-defined sector of the unit circle is
The method of claim 10, wherein f _s is a sampling frequency.