## Abstract

In the last decades, one of the main challenges in the study of heart rate variability (HRV) signals has been the quantification of the low-frequency (LF) and high-frequency (HF) components of the HRV spectrum during non-stationary events. At this regard, different time–frequency and time-varying approaches have been proposed with the aim to track the modification of the HRV spectra during ischaemic attacks, provocative stress testing, sleep or daily-life activities. The quantitative evaluation of power (and frequencies) of the LF and HF components has been approached in various ways depending on the selected time–frequency method. This paper is an excursus through the most common time–frequency/time-varying representation of the HRV signal with a special emphasis on the algorithms employed for the reliable quantification of the LF and HF parameters and their tracking.

## 1. Introduction

The heart rate variability (HRV) is the beat-to-beat oscillation of RR intervals around its mean value. It is the result of complex regulation mechanisms through which the autonomic nervous system (ANS) rules heart rate (HR) and keeps cardiovascular parameters within physiological ranges. Power spectral analysis of RR series is a well-established tool for the analysis of this kind of data, and it provides significant non-invasive information on ANS activity (Task Force 1996).

A typical spectrum of RR series is characterized by two main components (figure 1): a high-frequency (HF) component (synchronous with respiration, generally in the band 0.2–0.45 Hz) and a low-frequency (LF) one (centred approx. 0.1 Hz, band 0.03–0.15 Hz). While the HF peak is mainly modulated by the parasympathetic (vagal) branch of the ANS, producing the respiratory sinus arrhythmia (RSA), the LF component is mainly the result of sympathetic modulation, even if the influence of vagal tone may be observed. Consequently, the ratio of LF to HF power is often considered a measure of sympathovagal balance.

Traditional power spectral analysis of HRV data (employing either parametric or non-parametric approaches) implies the signal to be stationary. To overcome this limitation and to describe the changes in HRV spectra during transients of various nature, in the last two decades, tools for simultaneous time–frequency analysis of HRV data have been proposed. They include, among others, the short-time Fourier transform, time-varying autoregressive (AR) analysis (Bianchi *et al*. 1993), time–frequency representations (TFR; Novak & Novak 1993) and wavelet decomposition (Toledo *et al*. 2003). While in early 2000 (and before), the emphasis was focused on the selection of the most appropriate tool for time–frequency analysis of HRV data, more recently attention has been moved towards the development of algorithms for the reliable quantification of the LF and HF parameters. This paper aims to review the most popular tools for the TFR of the HRV signal and the algorithms employed for the tracking of LF and HF parameters.

## 2. Time-varying AR analysis

Time-varying analysis of HRV supposes the sequence of RR intervals, *x*(.), to be generated by a time-varying AR model of *p*-order(2.1)where is the discrete-time index; is the *parameter vector* (i.e. the model coefficients); is the *observation vector*; and *d*(*t*) is a scalar, zero-mean, white noise disturbance term of variance , . In equation (2.1), the parameter vector is made time dependent to take into account that system dynamics may change with time. Without loss of generality, we may describe these changes by(2.2)where *δθ*(*t*) is known as the *drift term*. If no *a priori* information is available, *δθ*(*t*) is assumed to be a zero-mean variable with and independent of .

The parameter vector is usually unknown and must be estimated from the data. Various algorithms are available for this task and those pertinent to the analysis of HRV series will be reviewed in the next sections. Here, it is sufficient to say that, when the current estimate of is available, time–frequency information is obtained by computing a time-varying spectrum(2.3)where , and *f*_{s} being the sampling frequency, and where , with *z*_{i}(*t*)'s being the model poles.

### (a) Recursive least-squares estimation

The problem of estimating can be solved by using the exponentially weighted least-squares (EWLS) algorithm, which is obtained by minimizing the loss function(2.4)where is the model prediction error. The coefficient *λ* is the *forgetting factor* that progressively discounts old data in favour of the new information and controls the algorithm responsiveness to model changes: the lower the forgetting factor is, the higher the algorithm responsiveness becomes. The minimum of (2.4) is achieved when(2.5)where(2.6)and where *Q*(*t*) is the so-called covariance matrix of the algorithm.

A crucial point is the selection of the appropriate *λ*, which is not a trivial task, and it may affect the performance of the estimation. Campi (1994) showed that the parameter estimation error is bounded(2.7)where *α*_{1} and *α*_{2} are the two coefficients depending on process characteristics and is the memory length of the algorithm indicating approximately how many samples are taken into account to produce the estimate.

The two terms on the right-hand side of (2.7) describe how the parameters' drift and the noise affect the estimation error. The first term—known as tracking error—increases linearly with *L*; thus, small *L* values should be selected to keep the error bounded in the presence of large (i.e. large variations in the parameter vector). The second term—known as steady-state error—is inversely proportional to *L*; thus, high *L* values are required to smooth the effect of noise. This theoretical framework supports the common knowledge that when *λ* is close to 1 (higher *L* value), the EWLS algorithm is insensitive to noise but is not able to track fast changes. Conversely, when *λ* is decreased (lower *L*), the algorithm is more reactive to parameter changes but becomes more sensitive to noise. Equation (2.7) shows that an optimum *L*, which minimizes the parameter estimation error, does exist, but its value depends on the parameters *α*_{1} and *α*_{2}, usually unknown in practice. In the analysis of HRV data, *λ*=0.98 is commonly employed for a memory length of approximately 50 samples (Mainardi *et al*. 2002).

#### (i) Recursive formulation

Even if equations (2.5) and (2.6) provide a way to estimate the model parameters, this solution is not efficient as it requires, at every step, the inversion of matrix *Q*(*t*). A more effective implementation makes use of a recursive updating of *P*(*t*)=*Q*(*t*)^{−1}, leading to the recursive least-squares algorithm (EW-RLS)(2.8)

For each new sample *x*(*t*), the parameter vector is updated by a quantity proportional to the model prediction error. The larger the error, the bigger the variation in model parameters and, in a picturesque way, we can say that the algorithm is learning from its own mistakes. A key role is played by *K*(*t*), the *gain* of the algorithm (updated step by step), which defines the property of the algorithm and modulates the extent to which prediction errors are reflected in the update of model parameters.

#### (ii) Whale forgetting

In order to overcome the limitation imposed by exponential forgetting, a more generalized recursion rule can be considered (Bianchi *et al*. 1997)(2.9)where *c*_{1} and *c*_{2} are the two coefficients such that the two roots of lie in (0, 1). The advantage of this generalized recursion is that the compromise between responsiveness and noise insensitivity can be made less dramatic than that in (2.7). The resulting weighting function (figure 2) is comparable to the EW-RLS ones in discharging past samples (i.e. the algorithm is equally reactive to the changes in the parameters), but the average of most recent data allows to make the algorithm less sensitive to noise. Owing to the particular shape of the resulting weighting function, the algorithm is known as the *whale-forgetting* least-squares method.

### (b) Smoothed Kalman filter

Equations (2.1) and (2.2) form the state-space model for the process *x*(.), where the AR model parameters represent the states. Thus, the minimum mean square estimator of , given the observations *φ*(*t*), can be obtained using the Kalman filter(2.10)where is the observation noise covariance and is the state noise covariance. These two terms have to be estimated. In the analysis of HRV data, Tarvainen *et al*. (2006) suggested to update using the following recursion with *γ*=0.95, and to estimate noise covariance as , being a running estimate of the variance of the HRV signal. The latter normalization makes the updating independent of the signal amplitude. The role of *μ* is analogous to the role of the forgetting factor in EW-RLS. In fact, by increasing *μ*, the adaptation speed is increased, but at the expense of an augmented variance of the state estimates. Thus, *μ* should be specified in such a way that a desired balance between the filter adaptation and the estimate variance is obtained. The analogy between *μ* and *λ* becomes evident if one observes the similarity between algorithms (2.8) and (2.10). For detailed description of the relationships between RLS and Kalman filter, see Sayed & Kailath (1994) and Lippuner & Moschytz (2004).

If a small lag in the processing is allowed and future observations can be used, a more accurate estimate of the state may be obtained using the so-called Kalman *smoother*. Using a fixed time interval *T* and for 1≤*t*≤*T*, the smoothing recursions for the state estimate is(2.11)where . The philosophy of the smoother is to increase the robustness of the estimates through a weighting combination of previous/past estimates.

### (c) HRV parameter estimation

One of the advantages of using a time-varying AR model is the possibility to extract the LF and HF parameters using a spectral decomposition method. The method was originally proposed by Zetterberg (1969) for EEG analysis and successively applied to the decomposition of the HRV spectrum (Baselli *et al*. 1987).

Using Cauchy's residue theorem, the signal power, i.e. the area under the AR spectrum (2.3), can be divided into a sum of *p* components(2.12)where the quantity(2.13)is the pole residue. Consequently, the spectrum can be decomposed into bell-shaped curves, the spectral components (figure 1), resulting in(2.14)

The centre frequency *f* and the power *P* of each spectral component can be computed as(2.15)1where ∠(.) is the phase expressed in radians, and *μ*=2 for complex pole pairs and *μ*=1 for real ones. When dealing with the HRV spectra, the calculation of the LF and HF powers may be performed by summing up the power values of those spectral components whose centre frequency lies within the chosen LF and HF bands, respectively.

Evaluation of equation (2.15) requires the computation of AR poles. Traditional factorization methods for pole estimation are time-consuming procedures, not appropriate for continuous pole detection: they need to be restarted each time a new coefficient set is available, and they do not take into account the previous estimated pole values. Therefore, pole-tracking algorithms were proposed for the analysis of HRV data (Mainardi *et al*. 1995; Meste *et al*. 2005).

The simplest way to evaluate the effect of parameter changes in the pole positions is to linearize the relationship(2.16)obtaining(2.17)

Using (2.17), for each new parameter vector, the pole positions can be updated as (Mainardi *et al*. 1995)(2.18)

The above relation is based on a linearization procedure that applies only to small coefficient changes. As the estimation procedure continues, error accumulation can corrupt the pole estimations. For this reason, once the poles are updated, the parameter vector is recalculated on the basis of the new *z*_{i}'s. The procedure may be iterated, if needed. The provides a measure of the sensitivity of the *j*th pole position to the changes in the *k*th model coefficient. This sensitivity depends on the reciprocal position of the poles and increases when the pole distance is reduced. Thus, in the case of closely grouped poles, even small changes in the parameters may overshoot pole changes (Oppenheim & Schafer 1975). Finally, it can be demonstrated that the algorithm is unable to generate a complex pole pair from two real poles; the desired number of complex poles must be provided during initialization.

To overcome these limitations, a second algorithm was presented, which is based on the iteration procedure of the well-known Bairstow algorithm, commonly employed to compute the roots of a polynomial function. Details and comparisons between the two methods applied to HRV data can be found in Mainardi *et al*. (1995).

In the work of Meste *et al*. (2005), in order to avoid problems related to variations in pole configurations (a complex pole pair becoming two real poles or vice versa), the AR model of even order is considered and the following factorization is used:(2.19)

Tracking is then performed on coefficients *b*_{i} instead of poles *z*_{i}, and reflected on pole movements using classical root calculation formula for the second-order polynomials. A formulation analogue to (2.18) was obtained and the efficient way to compute partial derivative is presented using a matrix formed by multiple convolutions. Finally, the authors propose an evolutive AR modelling (i.e. a model in which the *a*_{i} coefficients are constrained to be a linear combination of some basis functions) to provide a time-continuous description of the model coefficients and obtain a robust track.

## 3. Quadratic TFR

The time–frequency distribution of the signal energy may be obtained by quadratic TFR, among which the Wigner–Ville (WV) distribution plays a central role. It is defined as the Fourier transform of the instantaneous autocorrelation function (ACF), sometimes referred to as the central covariance function2(3.1)

The WV satisfies several desirable mathematical properties for a quadratic TFR (Hlawatsch & Boudreaux-Bartels 1992): it is real valued and preserves time and frequency shift of the signal; it satisfies the marginal properties, and thus the frequency (or time) integral of WV corresponds to the signal's instantaneous power (or power spectrum); and also the instantaneous frequency can be estimated from the first moment of the WV. The properties of (3.1) have been extensively discussed in Classen & Mecklenbrauker (1980).

Unfortunately, WV is affected by the interference terms (or cross-terms) that prohibit the interpretation of the WV as a pointwise time–frequency energy density and widely limit the readability of the distribution. The cross-terms oscillate in the *t*–*f* plane and can be attenuated by using kernels *ψ*, obtaining the *Cohen's class* (Cohen 1989)(3.2)where it is apparent that every member of the class may be interpreted as a two-dimensional, filtered version of the WV and is unequivocally defined by the kernel *ψ*. Among the various transformations belonging to Cohen's class, the smoothed pseudo-Wigner–Ville (SPWV) has been largely employed for the analysis of HRV data.

### (a) Smoothed pseudo-Wigner–Ville

In its discrete-time and discrete-frequency formulations, the SPWV is defined as(3.3)where *f*_{m}=*m*/*K*; *h*(*k*) is the frequency smoothing symmetric normed window of length 2*K*−1; *g*(*p*) is the time smoothing symmetric normed window of length 2*M*−1; and .

The properties of the SPWV depend on the characteristics of *g*(.) and *h*(.), which control time–frequency resolutions of the distribution. In respect to other representations of the Cohen's class, the SPWV is characterized by a separable kernel, in which time and frequency smoothing can be tuned independently, resulting in a great flexibility in the selection of smoothing criteria, application and efficient computation. This aspect is also relevant for the analysis of HRV series where the requirements for time–frequency resolutions are different. While frequency resolution is not problematic, as the spectrum of HRV is characterized by a few spectral components (mainly the LF and HF ones), time resolution should be increased in order to track spectral changes when they occur. In HRV data, therefore, time resolution is stressed, while frequency resolution is not.

The selection of the most appropriate smoothing windows for the analysis of cardiovascular variability series has been originally investigated by Novak & Novak (1993), using rest-to-tilt data and controlled reparation experiments. They suggested using a rectangular time-domain smoothing (window length between 5 and 25 samples) and a Gaussian window for frequency smoothing ( with *α*=2.5 and *N*=128). Comparing the performance of the SPWV with other TFRs, Pola *et al*. (1996) employed a Gaussian (23 samples length) and a hamming (128 samples) window for time- and frequency-domain filtering, respectively. In both works, the length of the time window is short, i.e. time resolution is stressed and limited to a few tens of seconds.

### (b) HRV parameter extraction

The parameters of the LF and HF components are usually extracted from the SPWV distribution using the marginal and the first momentum properties (even if they are not strictly valid), thus obtaining, for the LF components (analogous formulation can be used for HF one),(3.4)and(3.5)where the LF band is delimited by . Alternatively, the instantaneous frequency is obtained as the maximum in the band(3.6)

More recently, the work of Mainardi *et al*. (2004) has proposed a SPWV *decomposition method* for the quantification of the LF and HF contributions. The method is based on the decomposition of the instantaneous ACF and it assumes that the instantaneous frequencies of the LF and HF components are, locally, linearly varying with time. In particular, we assume that the analytic HRV signal is described as(3.7)where *A*_{LF} and *A*_{HF} are the instantaneous amplitudes of the LF and HF components, whose frequencies are and , respectively. If we compute *r*_{x}(*n*, *k*) for the signal in (3.7), we get(3.8)where the last term is the interference term, which is located in between the two components and oscillates in the time–frequency plane as expected (Hlawatsch & Boudreaux-Bartels 1992). Equation (3.8) is therefore composed by the sum of complex sinusoids including both the signal contributions and the interference terms. Note that the signal components contain information on the instantaneous frequencies and amplitudes of the LF and HF components. Equation (3.8) has been derived for two components only, but the extension to an *N*-component signal is straightforward using the quadratic superimposition principle (Flandrin 1984).

Using an exponential window for frequency smoothing and a rectangular window for time smoothing, it can be shown that the quantity(3.9)of equation (3.3) can be still approximated with a sum of complex sinusoids, where the influences of cross-terms are vanished by the application of smoothing windows and where the remaining components should be mainly related to signal contributions. For more details on how the assumption holds for a rectangular time-smoothing window, see Bailón *et al*. (2006*a*).

The problem of estimating the HRV spectral parameters is now reduced to the problem of estimating the *A*_{i}(*t*), *f*_{i}(*t*), *b*_{i} of (3.9). This can be solved using, for example, the method proposed by Kumaresan & Tufts (1982), which provides accurate detection and estimation of exponentially damped sinusoidal signals. The method combines linear backward prediction and singular value decomposition and it is particularly robust against noise. An example of decomposition of the SPWV is shown in figure 3.

Finally, it is noteworthy that *a priori* information on the HRV spectral components can be easily inserted as constraints in the estimation of parameters of (3.9). The inclusion of knowledge about respiratory frequency is discussed in Bailón *et al*. (2006*a*).

## 4. Wavelet analysis

Wavelet analysis has become a very popular tool for time–frequency signal characterization. The continuous wavelet transform (CWT) of a signal is defined as (Rioul & Vetterli 1991)(4.1)and it measures the similarity between the signal and a family of functions generated by scaling (contracting or dilating) and translating a function prototype *ψ*(.), called *mother wavelet*. The parameter *τ* gives the position of the wavelet in time, while *α* is the scaling factor that governs its frequency.

Using an appropriate selection of *ψ*(.), the inverse CWT exists and a signal can be represented by its decomposition into wavelets(4.2)whereis a constant that depends on *ψ*(*t*) only, and *Ψ*(*ω*) being the Fourier transform of the wavelets. The existence of (4.2) is guaranteed when *C*_{ψ}<+∞ (the *admissibility condition*). This condition imposes *Ψ*(0)=0, i.e. the wavelets have zero mean and are oscillating around the *t*-axis. In general, as |*Ψ*(*ω*)|→0 when *ω*→+∞, wavelets have a bandpass-like spectrum.

Using the Parseval equality, we rewrite relation (4.1) in terms of the Fourier transform of *x*(.) and *ψ*(.)(4.3)which provides a clear interpretation of the CWT in terms of a constant-*Q*-filterbank: a rescaling by *α* in the time domain becomes a rescaling by 1/*α* in the frequency domain, which has the effect of moving the centre frequency of the wavelet bandpass filter from *f*_{ψ} to *f*_{ψ}/*α*, with a similar rescaling of the 3 dB bandwidth, which is reduced when the scale increases. Therefore, the very appealing feature of the wavelet is that *α* makes the time resolution dependent on frequency (scale; and vice versa): at larger scales, a dilated wavelet is used to explore the signal reducing the time resolution but increasing the frequency localization of the components; at lower scales, the contraction of the wavelet provides timely location of HF details, at the expense of frequency resolution.

For the purpose of time–frequency analysis of HRV, we make use of the fundamental observation that CWT is isometric (Grossmann & Morlet 1984) and it preserves the signal energy *E*_{x}. We may write(4.4)which evidences that the squared modulus of wavelet coefficients can be interpreted in terms of energy distribution, where represents the ‘energy’ of *x* at lag *τ* and scale *α*. Such a distribution is known as *scalogram*. A detailed discussion of the CWT as energy distribution can be found in Rioul & Flandrin (1992).

Since the representation of the CWT is redundant, wavelet analysis is performed on a discrete set of time-scale parameters. The most common choice set is *τ*=*kα* and *α*=2^{−j}, where *j*, *k* are the integers (i.e. the CWT is sampled on a ‘dyadic’ grid) so that .

If the basis functions *ψ*_{j,k} form an orthonormal basis, a signal can be represented exactly as the weighted sum of these basis functions (Rioul & Vetterli 1991) , where the *d*_{j}(*k*)'s are the wavelet coefficients that can be computed using the discrete wavelet transform (DWT). In this case, the signal energy is obtained as .

### (a) HRV parameter extraction

On the basis of relation (4.4), the quantification of the LF and HF components is generally obtained by integrating the squared modulus of the wavelet coefficients over the desired frequency band . To do that, we have to move from scales to frequency and interpret the time-scale map in terms of a time–frequency distribution. We observe that the point in the time-scale plane is mapped in the point (*t*, *f*), and then using equation (4.4), we may write(4.5)which can be used for the computation of instantaneous power in a given frequency band. An example of tracking of HRV parameters during reperfusion after myocardial infarction (MI) is shown in figure 4.

When the DWT is considered, wavelet expansion divides the frequency axis into dyadic blocks and the frequency contents of each details *d*_{j} is . Instantaneous power of the LF and HF components is measured by(4.6)where the summation is extended to those scales that cover the LF and HF bands. These scales depend on the sampling frequency of the signal. For example, given *f*_{s}=2 Hz, the LF band will be covered by scales *j*=4,5, resulting in a band interval , while the HF band will be covered by scales *j*=2, 3, resulting in the band interval . Using a larger band for HF is generally not problematic, because only one spectral component (and its harmonic) exists at frequencies greater than 0.18 Hz (Toledo *et al*. 2003). However, this is only appropriate when the respiratory frequency lies within this band, but counter-examples exist such as during stress testing (Bailón *et al*. 2007).

## 5. Other approaches

A few other strategies have been pursued for the analysis of the HRV signals. Two of them are worth mentioning: the empirical mode decomposition (EMD) and complex demodulation (CD) techniques. In the EMD, the signal is decomposed as , where *m*_{N} represents the residual trend and *r*_{n}(*t*) are the *intrinsic mode functions* (IMFs), which are extracted iteratively as follows: (i) identify the extrema (local minima and maxima) of *x*(*t*), (ii) interpolate between minima (maxima), obtaining some ‘envelope’ *e*_{min}(*t*) (*e*_{max}(*t*)), (iii) compute the average , (iv) extract the detail , and (v) iterate on the residual *r*(*t*).

In practice, the above procedure has to be refined by a *sifting* process (Huang *et al*. 1998). According to the way they are extracted, the IMFs are constrained to be zero-mean, amplitude modulation, frequency modulation waveforms. Even though they are not necessarily narrow-band signals—as they can be non-stationary—the IMF should be, locally, monocomponent signals. Therefore, if one computes the Hilbert transform of the IMF, , the tracking of instantaneous power and frequency of the IMF are obtained as |*a*(*t*)|^{2} and , respectively.

In the CD approach, the signal *x*(.) is supposed to have a component of interest in the band and the signal is demodulated (i.e. multiplied by exp(*j*−2*πf*_{0})), thus shifting the spectrum and centring the band of interest on the *f*=0 axis. By applying a low-pass filtering with a cut-off frequency *f*_{w}, the desired component, let us name it , is retrieved. If *f*_{w} is selected sufficiently small (i.e. we have band-limited signal), we may write and instantaneous powers and amplitude estimated as above. The procedure must be repeated for all the frequencies of interest. It is noteworthy that this is equivalent to the bandpass (band: ) signal and apply the Hilbert transform to it. Analogy with EMD appears: the CD provides a signal decomposition in which the IMFs are simply bandpass versions of *x*(.). The main limitation of the CD is related to the selection of *f*_{0}'s and the definition of *f*_{w}. This turns out to be a non-trivial task, especially for the HF components. Therefore, *a priori* information about the frequency location of the HRV components must be included. The information may be taken from TFR of the HRV signal (Monti *et al*. 2002) or from the processing of other biological signals (i.e. respiration).

## 6. Applications

Early applications of the above-mentioned approaches dealt with the assessment of ANS responses during standard autonomic tests, in which a selective stimulus elicits sympathetic (rest to tilting, cold pressure test, etc.) or parasympathetic (deep breathing, valsalva manoeuvre) responses. The pronounced, acute changes after the stimulus are important for grading the autonomic performances and for evaluating neuroregulatory disorders.

Early responses to tilt have been first studied by the SPWV (Jasson *et al*. 1997; Mainardi *et al*. 2002), evidencing both vagal and sympathetic changes: the almost instantaneous disappearance of the HF component was accompanied by a slower, progressive increase in the LF components (Jasson *et al*. 1997). The response was dependent on the amplitude of the stimulus: a faster (and larger) sympathetic response was observed for more pronounced tilt angles (Mainardi *et al*. 2002). When the tilt was prolonged, oscillations and non-stationarity characterized the trend of the LF components, being more evident in patients with syncope (Furlan *et al*. 1998). Autonomic unbalance before syncope has also been assessed using the time-varying method (Mainardi *et al*. 1997; Furlan *et al*. 1998) documenting a sympathetic withdrawal preceding syncope, which may be either sudden or progressive. This withdrawal has also been confirmed by the analysis of MSNA variability using the CD (Kamiya *et al*. 2005). More recently, the DWT has been applied to a variety of autonomic tests (Ducla-Soares *et al*. 2007) and documented to be effective for the quantification of the first stage of the autonomic performance. In general, dynamic response to stimulus may help to better differentiate the class of patients: chronic fatigue syndrome with or without orthostatic postural tachycardia (Yoshiuchi *et al*. 2004), patients with predisposition to syncope (Lepicovska *et al*. 1992) or normal versus mild hypertensive patients (Novak *et al*. 1994).

Another field of application is the study of relationships between ANS and myocardial ischaemia, either spontaneous or induced. An early result, obtained by time-varying methods (Bianchi *et al*. 1993), evidenced an LF increase within the 2 min preceding the ischaemic attack. When different ischaemic events are considered, different trends were observed (Bianchi *et al*. 1997; Vila *et al*. 1997), all documenting an important involvement of ANS in either generating or sustaining the attack. An acute variation in autonomic tone preceding the attack has been recently confirmed using wavelet analysis and associated with spontaneous coronary spasm in patients with variant angina (Tan *et al*. 2003). During dipyridamole-induced ischaemia, the trend of spectral parameter was obtained through time-varying models and evidenced, in positive tests, an increase of the LF component in coincidence with the echocardiographic marker of ischaemia, but preceding the ST displacement on the ECG (Petrucci *et al*. 1996). Using wavelets analysis, this LF increase was related to the degree of pathology (being more evident in patients with triple vessel disease) and directly correlated with the wall motion score index (Petretta *et al*. 1999). Finally, the CWT was used for explaining the patterns of cardiac rate control during reperfusion induced by thrombolysis in MI patients (Toledo *et al*. 2003).

The dynamic evolution of RSA and its links to ventilation were investigated during exercise using an evolutive AR model (Meste *et al*. 2005), evidencing a progressive reduction in RSA modulation during the early part or the exercise (low/medium workloads) followed by a successive increase up to exercise peak (high workloads). The RSA persistence at peak exercise is probably a non-autonomic mechanism related to mechanical stretch of the sinus node, in response to ventilation increase. Dynamic behaviour of the HF component and the respiratory frequency in relation to the anaerobic threshold (AT) was investigated using the Hilbert transform (Anosov *et al*. 2000) demonstrating significant changes in the AT region. Interestingly, ventilatory thresholds could be assessed from the HRV data during stress by time-varying analysis of RSA (Blain *et al*. 2005) and during incremental exhaustive running test by the SPWV (Cottin *et al*. 2007).

An emerging field of application is related to the evaluation of HRV modifications during either normal or pathological sleep for the detection of apnoea events (Roche *et al*. 2003), for the sleep staging or for the assessment of autonomic changes during sleep (Spicuzza *et al*. 2003; Mendez *et al*. 2008). Events of apnoea (hypopnoea) are accompanied by concomitant cyclic variations in HR (a bradycardia/tachycardia pattern), which can be used to detect it (Penzel 2003). WT was proposed as a potential diagnostic screening tool in patients with sleep apnoea–hypopnoea syndrome (Hilton *et al*. 1999), and was proved to be very well suited to recognize the sleep apnoea-specific cyclic variability of HR (Roche *et al*. 2003). A comparison of different methods for sleep apnoea detection using ECG, also including time–frequency methods, can be found in Penzel *et al*. (2002). Recently, indices derived from WT of HRV have been proposed as the indicator of the degree of sleep fragmentation (Sforza *et al*. 2007). Markedly different autonomic alterations in the modulation of HR during (and after) obstructive and central apnoeas in patients with sleep-disordered breathing have been observed using the WV transform, suggesting a surge in sympathetic modulation after the obstructive episodes (Spicuzza *et al*. 2003). Arousal-induced changes in the HRV have been studied using TV methods evidencing that the induced sympathoexcitatory cardiovascular effects are relatively long lasting, compared with MSNA recordings, and may accumulate if repetitive arousals occur in close succession (Blasi *et al*. 2003).

Other applications deal with the evaluation of sympathovagal balance before non-sustained ventricular tachycardia (Chen 2002) or the assessment of ANS during daily activity (Chan *et al*. 2006).

## 7. Discussion and conclusion

This paper presents an overview of the most common time–frequency methods for HRV analysis during non-stationary conditions. The attention was focused on the algorithms employed for the continuous quantification of the LF and HF indices, which ranged from the simple integration of power in defined frequency bands (scales) up to ad hoc decomposition methods designed for the HRV signal. A few examples were presented to emphasize the usefulness of these methods in extracting HRV indices during non-stationary situations and to stress advantages/shortcoming of the methods.

At this regard, the parametric time-varying analysis fits the pseudo-stochastic nature of the HRV signals, and allows great flexibility because time resolution (governed by the forgetting factor) and frequency resolution (depending on model order) can be tuned separately. The AR models provide smoother and narrow frequency peaks, permitting precise frequency identification and the recursive identification is useful for real-time applications. However, tracking of power is more problematic especially when the model poles are close to the unit circle because a small change in their position may result in a large spectral modification.

Computation efficiency and excellent time–frequency resolutions are among the advantages of the SPWVD methods. When these properties are combined with a decomposition method, an efficient tracking of the LF and HF parameters is obtained. The main limitation of this approach is the simplicity of the HRV model (3.7).

The appealing feature of the CWT is that time resolution depends on the frequency of interest and this allows a fast tracking of HF spectral power. However, power computation may be problematic when scales and HRV bands do not properly match or when more than one spectral component is present at a given scale. Despite these limitations, the number of HRV studies employing wavelet transform is considerably increasing in the last few years.

All the above-mentioned approaches, when correctly tuned, provide useful information about ANS dynamic, but further work is needed. A dynamic change in ANS control is usually accompanied by changes (sometimes very pronounced) in the mean HR resulting in the trends that need to be removed to emphasize the beat-to-beat variations. While this is automatically obtained using signal decomposition methods (WT and EMD), pre-processing is requested before the application of time–frequency and time-varying methods (Mainardi *et al*. 2002; Tarvainen *et al*. 2002). Another problem, which is often overlooked, is the selection of the most appropriate frequency band for the computation of the LF and HF parameters. In traditional spectral analysis, these bands are fixed, but this choice is obviously suboptimal when dynamical tracking has to be performed. In these situations, the definition of *dynamic* boundaries is important and different solutions have recently been proposed (Goren *et al*. 2006; Bailón *et al*. 2007).

In conclusion, the methods described in this paper were designed for HRV analysis, but they may be of potential interest in other applicative fields where the tracking of frequency and power of the spectral components is needed.

## Footnotes

One contribution of 13 to a Theme Issue ‘Signal processing in vital rhythms and signs’.

↵This approximation is valid when |

*z*_{i}| is sufficiently close to 1. More correctly, the maximum of the spectral component is located at*ω*satisfying .↵From this point on,

*t*represents the continuous time.- © 2008 The Royal Society