# Stochastic Modeling and Optimal Time-Frequency Estimation of Task-Related HRV
## Methods
### Inference Method  
Since a closed-form expression of the process distribution is not known, a maximum likelihood approach for estimating the model parameters is not feasible. In a previous contribution [23], we addressed the problem of fitting a real-valued LSP to sampled data by proposing the inference method HAnkel-Toeplitz Separation (HATS). The method hereby described is based on the same principle of dividing the inference problem into lower-dimensional sub-problems. The estimation procedure is summarized in Algorithm 1 and described in detail below. ```
Input : Zero mean time-series data, both respiratory signal and corresponding HRV. Output : Estimated chirp parameters \(\hat{m},\hat{d}\) and estimated LSP parameters \(\hat{a},\hat{b},\hat{c}\). Chirp estimation
1 Estimate the instantaneous frequency of the chirp from the the respiratory signal, by taking the frequency of maximum power in the chirp range from a TF representation of the signal;
2 Find \(\hat{m},\hat{d}\) and the corresponding covariance matrix \(\hat{F}_{I}\) by fitting the linear chirp to the instantaneous frequency of the respiratory data; LSP estimation
3 Compute the instantaneous power \(\hat{P}\) from the HRV data;
4 Find \(\hat{a},\hat{b}\) by fitting the model function \(q\) to \(\hat{P}\);
5 Find \(\hat{c}\) by fitting \(R\circ\hat{F}_{I}\) to \(\hat{R}_{F}=\hat{C}_{SCM}\odot\hat{Q}\) where \(R\) is the stationary covariance matrix of the model function \(r\), \(C_{SCM}\) is the estimated sample covariance matrix of the HRV data and \(\hat{Q}\) is the corresponding Hankel matrix based on the model \(q\) with parameters \(\hat{a},\hat{b}\); Return : \(\hat{m},\hat{d},\hat{a},\hat{b},\hat{c}\). ```  
**Algorithm 1**Inference for LSCP parameters  
The first step is estimating the instantaneous frequency of the chirp from the respiratory signal, by taking the frequency of maximum power in the chirp range (frequencies above 0.12 Hz) of the spectrogram with Hanning window of 128 samples (32 s). The parameters \(\hat{m}\) and \(\hat{d}\) are estimated through a robust least square fitting of the linear chirp \(mt-md\) to the estimated instantaneous frequency of the respiratory signal. A corresponding chirp covariance matrix \(\hat{F}_{I}\) is constructed according to the model Equation (2), where the covariance matrix values are \(\hat{F}_{I}\left(t_{1},t_{2}\right)=f_{I}\left(t_{1},t_{2}\right)\). The estimation of the chirp from the respiratory signal of each TP allows accounting for the individual differences in the interpretation of the task of breathing accordingly to the metronome. An example of linear chirp fitted to the instantaneous frequency extracted from the spectrogram is presented in Figure 3.  
At this point, a variation of the inference methods HATS can be used to compute the LSP parameters. First the function \(q\) in Equation (4) is fitted to the instantaneous power estimate \(\hat{P}(t)=x^{2}(t)\), where \(x(t)\) is the zero-mean adjusted HRV at time \(t\), as in Equation (6), obtaining the estimated parameters \(\hat{a}\) and \(\hat{b}\), from which the Hankel matrix \(\hat{Q}\), with values \(\hat{Q}\left(t_{1},t_{2}\right)=\hat{q}\left(\left(t_{1}+t_{2}\right)/2\right)\), can be calculated. Figure 2: Function q fitted to the mean instantaneous power of the 92 HRV time-series. The next step is computing the ordinary sample covariance matrix of the HRV data \(\hat{C}_{SCM}\). Thereafter, the Hadamard division of \(\hat{C}_{SCM}\) by \(\hat{Q}\) provides the estimate \(\hat{R}_{F}=\hat{C}_{SCM}\odot\hat{Q}\), where \(R_{F}\) corresponds to a non-stationary covariance matrix with values \(R_{F}\left(t_{1},t_{2}\right)=r\left(t_{1}-t_{2}\right)\cdot f_{1}\left(t_{1},t _{2}\right)\). The parameter \(\hat{c}\) is obtained by fitting the model stationary covariance matrix \(R\) obtained from Equation (5), Hadamard multiplied with the chirp matrix \(\hat{F}_{I}\), to \(\hat{R}_{F}\).