Modeling of Heart Rate Variability Using Time-Frequency Representations

: The heart rate variability signal is highly correlated with the respiration even at high workload exercise. It is also known that this phenomenon still exists during increasing exercise. In the current study, we managed to model this correlation during increasing exercise using the time varying integral pulse frequency modulation (TVIPFM) model that relates the mechanical modulation (MM) to the respiration and the cardiac rhythm. This modulation of the autonomic nervous system (ANS) is able to simultaneously decrease sympathetic and increase parasympathetic activity. The TVIPFM model takes into consideration the effect of the increasing exercise test, where the effect of a time-varying threshold on the heart period is studied. Our motivation is to analyze the heart rate variability(HRV) acquired by timevarying integral pulse frequency modulation using time frequency representations. The estimated autonomic nervous system (ANS) modulating signal is filtered throughout the respiration using a time varying filtering, during exercise stress testing. And after summing power of the filtered signal, we compare the power of the filtered modulation of the ANS obtained with different time frequency representations: smoothed pseudo Wigner–Ville representation, spectrogram and their reassignments. After that, we used a student t-test p < 0.01 to compare the power of heart rate variability in the frequency band of respiration and elsewhere.


Introduction
The Heart Rate Variability is the time variation of heartbeats. It reflects the regulations of the Autonomic nervous System [1,2]. The fluctuations of heart rate are modeled by mathematical chaos [3]. The variability leads to the flexibility to rapidly contend with a changing environment. Modeling biological systems reveals spatial and temporal complexity, so any disorder changes this complexity [4]. Higher magnitude of the HRV is not always a sign of healthy biological conditions. For instance, any conduction problem in the heart, makes HRV amplitudes higher. To study the HRV, we need to start with studying electrocardiogram (ECG) waves. A good interpretation of these waves shows whether high HRV values are the consequence of cardiac problems like atrial fibrillation [5] or other side effects that could impact heart rhythm [6]. A normal magnitude of HRV could be interpreted as a sign of adaptability or resilience. Higher magnitudes of HRV are due to the effect of executive functions by the prefrontal cortex [7,8].
The Integral pulse frequency modulation (IPFM) model explains the influence of the heart rate by the autonomic nervous system [9] in several physiological conditions: at rest or during exercise [10,11]. However, the IPFM model considers the heart period as a constant [12,13] which means that it considers a constant mean heart period [14]. That's why we use the time varying threshold IPFM (TVIPFM) model to analyze the heart rate variability during exercise. Applying this model leads to the approximation of a modulating signal taking into account the time varying threshold. Many studies highlighted the close relationship between respiration and HRV [15][16][17]. That's why in this paper, we analyze non-stationary (HRV) signals using time varying filtering based on the frequency band of respiration.
Our study is in stress testing exercise, permitting us to look the respiratory information to estimate parasympathetic activity. Heart rate variability was analyzed in time frequency domain. The Wigner-Ville [18] and the Smoothed Pseudo Wigner-Ville Distributions (SPWVD) have already been used to analyze biological signals [19], including HRV and respiration oscillations. The Smoothed Pseudo Wigner-Ville Distribution is a good technique to study large band nonstationary signals. This technique allows us to determine time-intervals, where the heart rate variability is excited and where it is not. The knowledge of these intervals in time could also be used to find a direct correlation, graphically and analytically, between the heart rate variability and the respiration. But unfortunately, the SPWVD has a cross term which causes interference, that's why we need a reduction of cross terms, which is obtained by use of analytical signal with no frequencies instead of real signal. This procedure is based on the Hilbert Transform, which is a linear transform that extends a real signal to a complex one, to satisfy the equations of Cauchy Riemann.
This method could be applied during exercise with a high stress test, but it could be applied during other conditions.

Methods and Materials
In our study, one of the time frequency representations used in non-stationary conditions, is the spectrogram. The spectrogram is as below [20,21] h is a smoothing window. The spectrogram is the squared magnitude of the short Fourier Transform.
Another time-frequency representation studied is the Pseudo Wigner-Ville Distribution (PWVD) since it uses a short term window, whose effect is to operate on the Wigner-Ville distribution a smoothing frequency, which amounts to limit the interference terms only components that are simultaneously present in the window. The PWVD is defined by [22]: with p is a smoothing frequency window, if this window is factorable, we have: And the result obtained is the pseudo Wigner-Ville, inspired from the Wigner-Ville. This analysis is a slippery short term analysis, which therefore resembles the spectrogram. To find a relation between spectrogram and SPWVD, we introduce the weighted delayed signal: permitting as to define the PWVD: Spectrogram can be written as: Using the spectrogram as a smoothed version of Wigner-Ville distribution presents two major weaknesses: the first one is that passing to spectrogram implies the loss of theoretical advantageous properties of Wigner-Ville distribution, and the second one is that the smoothing introduced obeys the constraint of Heisenberg-Gabor between time and frequency. Since the smoothing time frequency associated to spectrogram has only one degree of freedom.
Considering the two dimensions obtained by the time and the frequency domain, improvement is possible if we move to a smoothing with two "degrees of freedom" respectively to the time and frequency. Both distributions, spectrogram and pseudo Wigner-Ville, have the same principle, a signal segment taken with a short-term window. After, both undergo a Fourier Transform, followed by a quadratic operation. In this paper we use the SPWVD and the balance between time and frequency filtering can be done independently. One of the drawbacks of the SPWVD representation is the presence of interferences.
After we will use the indicator function named GAB resp in the time frequency domain showed as below

The Time Varying Threshold Integral Pulse Frequency Modulation Model
The TVIPFM model considers m(t) as the ANS control. The first beat happens at time t 0 = 0, so t k is solution of [2,23]: the instantaneous heart rate is defined as 1+m(t) Supposing that the oscillating term m(t) T(t) is with higher frequency than the term 1 and the heart rate variability is represented by After dividing the HRV d HRV (t) by d HRM (t) we obtain m(t). The variations of the timevarying threshold T AC (t) are smaller than the mean value T DC . Rewriting the instantaneous heart rate Assuming m(t) 1 we have: We assume This mean heart rate obtained by the TVIPFM model, is more adequate to stress tests and exercise, taking into consideration a mean heart period varying function of time, changing between two successive heartbeats. So First, we calculate k (t) from the pairs (t k , k), then we low-pass filter d HR (t) to obtain the d HRV (t) term. Therefore, we get m (t). Finding this modulating signal using our TVIPFM model combined with different time frequency representation is aiming to calculate the effect of the autonomic nervous system on the sinoatrial node. This control is synchronized with the respiratory sinus arrhythmia. This modulation m (t) plays a key role to understand different pathologies and dysfunctions of the autonomic nervous system to monitor the cardiac rhythm.

Statistical Analysis
The correlation between the respiration and the ANS modulation is investigated by comparing the instantaneous power throughout the respiration using the gabarit function, and the instantaneous power at the rest of the frequencies, other than the frequency band of respiration, using a student t-test. Statistical results have been shown at the time frequency domain, during exercise stress testing.
We calculate the instantaneous power for each time, around the respiration at the first time, and elsewhere at the second time. The student t-test significance was fixed for p < 0.01.

Simulation Study
In Fig. 1, we generate a signal containing three parts to represent the modulating signal. The first part is a sinusoidal signal with constant frequency around 100 Hz and a constant magnitude equal to one. In the second part, we generate a chirp with instantaneous frequency varying linearly from 100 to 350 Hz, and with a magnitude, increasing linearly from one to five, and the third part is like the first one. We define the gabarit function as a function which is equal to one around the frequency band of respiration and zero elsewhere, mathematically, it can be considered as the Indicator function which is equal to one around the respiratory frequency and zero elsewhere. The gabarit plays a key role in filtering, the representation of the obtained signal will be multiplied by the gabarit function in the time frequency domain for spectrogram, SPWVD and their reassignments. After the filtered signal obtained at the frequency band of respiration, will be studied to compare which representation is better, and contains more power. In this simulated signal, we consider that the respiratory signal is concentrated around the frequencies of the signal in each part, the aim of this simulation is to show how the gabarit is well adapted to our simulated modulation [24].
In Fig. 2, after representing in the time frequency domain with spectrogram, SPWVD and both their reassignments [25] for a simulated signal, the result obtained will be multiplied by the GAB resp and after that we integrate over the frequencies weighting by the number of the points used in the gabarit function. The obtained result, is quadratic in the time domain, to compare it with the original signal, we take the root square of this quantity, and we compare with the envelope of the magnitude of the starting signal as shown in Fig. 3.

Real Signals During Exercise Stress Testing
In our study, we need the instants t k , from these beat occurrences, first we obtain the instantaneous heart rate, then we get the time varying mean heart rate from the instantaneous heart rate, using a high-pass filtering, then we obtain the heart rate variability by subtracting the time varying mean heart rate from the instantaneous heart rate. Finally, the obtained HRV will be corrected by the time varying mean heart rate, and we obtain the modulating signal using the TVIPFM model.
The technique proposed in this paper is the time varying signal filtering using the gabarit function at the time frequency domain around the respiration.
The SPWVD signal X (t, υ) of the original signal x(t) is after integrated respect to the frequency υ, for spectrogram and the SPWVD and their reassignments, we expect to have the instantaneous power |x (t) | 2 , and  So, we obtain the magnitude of the signal x (t) after the time varying filtering, that's why we apply the root square on this quantity to retrieve the magnitude of the real signal.
In Fig. 5, we represent the envelope of the magnitude of the modulating signal with different time frequency representations like SPWVD, spectrogram and their reassignments. Results obtained prove the importance of our time varying filtering technique to retrieve the modulation from the heart rate variability due to the respiratory sinus arrhythmia. The use of a time varying model with a time varying filtering around the respiration, gives us the best envelope using the spectrogram method.  There is a high correlation between the respiration and the heart rate variability, with significance p < 0.01. In this figure, for each time frequency representation, the instantaneous power at the frequency band of respiration and elsewhere are computed. We obtain two sets of data. We use a t-test with a window of 200 samples. The curves spwv, r-spwv, spe and r-spe plotted in Fig. 6 represent the result of the t-test, the horizontal straight line represents the threshold of the significance level. The presence of the instantaneous power in the respiratory frequency band is always shown p < 0.01 with each time frequency representation. Therefore, we can use our technique of filtering to retrieve the envelopes of the magnitude of the modulating signal corresponding to the ANS modulation. The spectrogram representation, shows the most significant correlation between the modulating signal and the respiration.
The correlation modeled between the respiration and the HRV plays a key role to understand the control of the ANS on the heart. A recent work [26] studied the link between respiration and HRV during yoga breathing practice. After six months of yoga, HRV showed improvement towards the parasympathetic domain. This study was performed in the time and frequency domains. Modeling the HRV coupled with the respiration has several applications. One current study [27] was applications was about tracking the cardiorespiratory load of firefighters reflecting the respiratory metabolism state to ensure their safety. This study needs to take into account a respiration during the intense physical exercise. The reader may refer to [28], where the HRV is positively linked with compassion. The study HRV associated with soothing emotions.

Conclusion
In this paper, we used the TVIPFM to obtain the modulating signal of the heart rate variability filtered with respect to the respiration frequency. This filtering is assumed in time frequency domain, using different representations like spectrogram, the smoothed pseudo Wigner-Ville distribution and their reassignments. Differences between magnitudes of instantaneous power in all of the representations are due to the reassignment principle moving power dispersed, to reassign it to the nearest point. After that, we used a statistical test to explain our choice to filtrate around respiration, so we involve the student t-test to compare the instantaneous power in the respiratory frequency band and elsewhere using the gabarit function. The significance of the t-test is set at p < 0.01 showing that all time frequency representations used prove the existence of the power in the respiratory frequency band.
HRV has been also studied with Kalman smoothed method to estimate time varying characteristics of HRV. The proposed method, supposes that HRV is obtained by auto regressive model. Using a Kalman smoothed algorithm, parameters are estimated. In this method, the time varying spectrum is a continuous function of frequency, so we can evaluate it at any desired frequency up to the Nyquist frequency f s 2 . The advantage of the Kalman smoother algorithm method is that the spectrum is separable.
However, one of the limits of this method is that the characteristics of the Kalman smoother spectrum depends strongly on the order of the auto regressive model. For example, if the used model is with a higher order, it can lead to some interference in the spectrum. Studying HRV during exercise imposes the use of a non-stationary mean time varying period. So, we studied the HRV using our TVIPFM model, also we used the method of time varying filtering. Finally, with this study, we satisfy the increasing exercise conditions. Funding Statement: This work was supported by College of Engineering and Technology, the American University of the Middle East, Kuwait. Homepage: https://www.aum.edu.kw.

Conflicts of Interest:
The authors declare that they have no conflicts of interest to report regarding the present study.