## Abstract

The continuous wavelet transform (CWT) is specifically efficient in the analysis of transient and non-stationary signals. As such, it has become a powerful candidate for time–frequency analysis of cardiovascular variability. CWT has already been established as a valid tool for the analysis of single cardiovascular signals, providing additional insights into the autonomous nervous system (ANS) activity and its control mechanism. Intercorrelation between cardiovascular signals elucidates the function of ANS central control and the peripheral reflex mechanisms. Wavelet transform coherence (WTC) can provide insight into the transient linear order of the regulatory mechanisms, via the computation of time–frequency maps of the time-variant coherence. This paper presents a framework for applying WTC for quantitative analysis of coherence in cardiovascular variability research. Computer simulations were performed to estimate the accuracy of the WTC estimates and a method for determining the coherence threshold for specific frequency band was developed and evaluated. Finally, we demonstrated, in two representative situations, the dynamic behaviour of respiration sinus arrhythmia through the analysis of the WTC between heart rate and respiration signals. This emphasizes that CWT and its application to WTC is a useful tool for dynamic analysis of cardiovascular variability.

## 1. Introduction

Spectral analysis of heart rate (HR), respiration and blood pressure (BP) is a well-established tool for the non-invasive investigation of cardiovascular control mechanisms (Akselrod *et al.* 1981; Task Force 1996; Malliani 2000). The interrelation and influence between HR, BP and respiratory signals provides information about the functionality of autonomic nervous system (ANS) central control, the peripheral reflex mechanisms, and the origin and characteristics of the driving sources. There are two major inter-modulations with HR: respiratory sinus arrhythmia (RSA; Hainsworth & Mark 1993) and baroreflex, which describes the changes occurring in HR following changes in arterial BP. In addition to these two bivariate basic regulatory mechanisms, there are descriptions of additional multivariate regulatory mechanisms (Triedman & Saul 1994; Parati *et al.* 1995; Baselli *et al.* 2001; Elstad *et al.* 2001). Even though the relationships between the ANS and the responses of the cardiovascular and respiratory system still require further understanding, it is a reasonable assumption to state that some of the regulatory mechanisms have linear characteristics. Therefore, many methods for investigating the regulatory system examine linear relationships and their quantitative properties (linear feedback, correlation, phase, etc.; Saul *et al.* 1991). Other methods focus on nonlinear characteristics (Voss *et al.* 2009). The coherence function is a straightforward method used to assess the existence and strength of linear coupling between two signals in the frequency domain (Kay 1988). As such, it is often applied to the analysis of the ANS functionality via the examination of cardiovascular and respiratory behaviour (Saul *et al.* 1991; Stanley *et al.* 1996). A major limitation of the classical spectral analysis is the steady-state assumption. This stipulation limits the analysis to periods or conditions with evident steady-state external conditions. Therefore, it requires, at some level, disregarding of the dynamic, complex and noisy characteristics of the continuous physiological adjustments required to maintain homeostasis. Both parametric and non-parametric time–frequency methods have been introduced to detect dynamic behaviour such as: the short-time Fourier transforms (STFT); time-dependent autoregressive models; Wigner–Ville distribution (Cerutti *et al.* 2001); and the selective discrete Fourier transform algorithm (Keselbrener & Akselrod 1996). The wavelet transform (Daubechies 1992) also acts as a method for analysing cardiovascular signals in the time–frequency or time–scale domain and has already proven itself as a powerful and robust tool in the analysis of transient phenomena of the ANS (Pichot *et al.* 1999; Davrath *et al.* 2003; Toledo *et al.* 2003). Although not necessarily superior to other methods, the wavelet transform is specifically efficient in the analysis of transient and non-stationary signals by enabling simultaneous optimal interpretation of spectral and temporal information using a window of variable width. This provides us with a more flexible temporal–spectral aspect of the transform than the STFT. Desired resolution can be obtained simultaneously for each signal feature; higher temporal resolution for higher frequency and higher spatial resolution for lower frequency. The use of both discrete and continuous wavelet transforms (CWTs) in cardiovascular signals has been reviewed by Addison (2005).

The use of wavelet transform in bivariate analysis of cardiovascular signals has the potential to offer additional information about the complex functioning of autonomic regulatory mechanisms. The wavelet transform coherence (WTC) is a known, though infrequently used, tool in scientific areas, where random-like, yet deterministic, signals are created by complex and not fully understood mechanisms. Areas such as weather in geophysics (Grinsted *et al.* 2004), turbulence in plasma physics (Van Milligen *et al.* 1997) and brain dynamics in physiology (Klein *et al.* 2006) are using wavelet coherence in order to find transient correlations between signals. Unlike the use of coherence in cardiovascular variability, the wavelet coherence was mainly used for the purpose of finding correlated areas between signals that are uncorrelated most of the time.

The use of WTC in cardiovascular variability research has already been introduced (Keissar *et al.* 2006) and further explored (Ostlund *et al.* 2007; Cnockaert *et al.* 2008), focusing on the dynamic analysis of RSA.

The purpose of the following work is to establish the wavelet coherence as a fundamental method in physiological research, particularly for use in the analysis of ANS control. We have reviewed and adapted the Morlet wavelet transform for the quantitative analysis of wavelet coherence of cardiovascular signals. Computer simulations were performed to asses the bias and standard deviation (s.d.) of the WTC estimates over variable wavelet coefficient and coherence real value. Next, a method for determining the coherence threshold for a specific frequency band was developed and false positive and false negative error rates were used to evaluate this threshold. Finally, we evaluated, in two representative situations, the dynamic behaviour of RSA through the analysis of the WTC between HR and respiration signals.

## 2. Methods

### (a) Mathematical overview

This section describes the method of wavelet analysis based on the Morlet wavelet and gives a mathematical overview for the analysis of the WTC. A more detailed overview was given in previous work (Keissar *et al.* 2006) and is based on the work of Torrence & Compo (1998).

WT decomposes a time-series signal into a time–scale plane. The continuous wavelet function based on the Morlet wavelet function consists of a plane wave modulated by a Gaussian,(2.1)where *t* and *s* are time and scale, respectively. The Morlet coefficient *ω*_{0} shifts the balance between frequency resolution and time resolution. The Morlet coefficient is limited to *ω*_{0}≥6 owing to the admissibility condition (Farge 1992). We found that, for cardiorespiratory signals, the preferable choice of time and frequency resolution combination is 6≤*ω*_{0}≤30.

The CWT is defined as the convolution of a scaled parent wavelet function (2.1) with the analysed function *g*(*t*),(2.2)The discrete form of the CWT of a sequence *x*_{n} is(2.3)(2.4)where *Ψ*^{*} is the complex conjugate of the normalized wavelet function; *n* is the time-series index; and *δt* is the sampling time.

The wavelet power density estimator of *x*_{n} is defined as(2.5)where is the complex conjugate of the wavelet coefficient.

For the Morlet wavelet, the relationship between the scale and the frequency is established through the Fourier wavelength *λ* (Torrence & Compo 1998),(2.6)Since the analysis is usually applied to a discrete signal with a finite length, it should be noted that CWT has edge artefacts because of the wavelets' width in time. The cone of influence (COI) is the region of the spectrum in which the edge effect is significant and the power spectrum should be considered as dubious. The COI is chosen where the wavelet power for the discontinuity at the edge drops by a factor of *e*^{−2}, which is for the Morlet wavelet (derived from the power of the transform in equation (2.1)). The COI also gives a measure of the decorrelation time for a single spike in the time series. The COI is larger for larger *ω*_{0}.

WTC is the logical continuation of applying the WT as an estimate of the time-dependent spectrum. The naive form of coherence is not a valid option, since it is identically equal to 1 (Torrence & Webster 1999). Accordingly, some kind of smoothing integration is required even at the cost of reduced localization. We chose the smoothing operator described below, as suggested by Torrence & Webster (1999) and also applied by Grinsted *et al.* (2004).

The squared wavelet coherence estimator is defined as the squared absolute value of the smoothed cross-wavelet spectrum, normalized by the smoothed wavelet power spectrum of the two signals,(2.7)where and are the wavelet spectral density function; is the cross-wavelet spectrum (2.8); and the angular brackets indicate the smoothing operator. This definition of the wavelet coherence corresponds with the Fourier-based coherence and maintains its value between 0 and 1,(2.8)

The smoothing operator in (2.7) is achieved by a convolution in time and scale (2.9). The time convolution is performed with a Gaussian , which is the absolute value of the wavelet function in each scale. The time convolution will double the edge artefact to . The scale convolution is performed by a rectangular window with a length of *δj*_{0}.*s*, where *δj*_{0}=0.6 is the empirical scale decorrelation length for the Morlet wavelet (Torrence & Compo 1998),(2.9)

In equation (2.9), *c*_{1} and *c*_{2} are the normalization factors and Π is the rectangular function. It should be noted that the estimator characteristics, such as the significance threshold and the bias, strictly depend on the smoothing parameters.

### (b) Assessment of the WTC estimates and threshold levels

This section deals with validation of the WTC estimates using computer simulations and surrogate data, in order to draw any quantitative conclusion about the reliability of the estimates and the actual detection accuracy of the WTC.

#### (i) Theoretical coherence for a linear system

The theoretical coherence estimation is based on the model of a simple input–output of a linear time-invariant (LTI) system with the addition of uncorrelated noise, which accounts for all the nonlinear and uncorrelated effects. In the case of time–frequency analysis, we assume that the system properties are time invariant for each specific instance but can change as time progresses. In the case of our simulations, the transfer function of the system is constant and equals unity. Therefore, the theoretical coherence between two signals *X*(*t*) and *Y*(*t*) is given by (Bendat & Piersol 1971)(2.10)where *G*_{XX}(*f*) and *G*_{NN}(*f*) are the spectral density functions of the input signal *X*(*t*) and the noise that was added to the output signal *Y*(*t*), respectively. In order to apply the model for time–frequency analysis, we would use the wavelet power spectrum,(2.11)Equation (2.11) shows that the coherence is basically determined by the signal-to-noise ratio (SNR) value of the system. Therefore, the best controllable way to obtain different coherence levels is to change the variance of the output noise, keeping other parameters constant (Pinna & Maestri 2001; Porta *et al.* 2002).

#### (ii) Simulations for assessing bias and s.d.

In order to assess the statistical error of the WTC with regard to the theoretical coherence, computer simulations were performed. As our principal simulated signal (*X*(*t*)), we chose a bandpass white noise signal in the 0.15–0.4 Hz band to simulate an arbitrary respiration signal. *X*(*t*) signals were created by applying a low-pass filter, with a cut-off frequency of 0.4 Hz, followed by a high-pass filter (HPF), with 0.15 Hz cut-off over a unit amplitude white noise. In accordance with the LTI model of transfer function equals unity, the output signal *Y*(*t*) signal was obtained by adding uncorrelated white noise to *X*(*t*). The simulations were performed repeatedly with different noise levels (different amplitudes of the added white noise) to give *N*>1000 realizations for each coherence level, in 0.1 steps. Unless mentioned otherwise, the coherence levels were computed by averaging the coherence over the 0.15–0.4 Hz band (band coherence).

The bias and s.d. were obtained as(2.12)(2.13)where *k* is the theoretical coherence levels (*γ*^{2}) from 0.05 to 0.95 in 0.1 steps and *C*^{2}(*k*) is the WTC level associated with the theoretical coherence level equals *k*. The procedure was repeated for different wavelet coefficients (*ω*_{0}): 6, 10, 15, 20 and 30.

#### (iii) WTC threshold (for zero coherence) and error rates

The distribution of the coherence values of two uncoupled signals, using wavelets with a smoothing operator, as applied in this study, cannot be described analytically (Maraun & Kurths 2004). Therefore, in order to determine the existence of coupling, we chose a null hypothesis that the signals have zero coherence at the examined frequency. Using a magnitude threshold, all the time–frequency regions having WTC values, which are unlikely to occur if the two signals were uncoupled, will be considered as regions in which a linear coupling exists between the signals. The numerical approach was to use surrogate data for the estimation of the null hypothesis. The distribution was calculated for 10^{6} values, for each frequency, and the critical level was set at 95 per cent. As previously established (Faes *et al.* 2004), coherence threshold levels and their ability to predict error rate depend on the choice of the surrogate signals used and their similarity to the analysed signals. Therefore, the evaluation of the threshold for the band coherence is performed using uncorrelated signals created in the same way as the simulated signals described above, but with SNR≪1, which gives a theoretical band coherence of virtually zero. Using surrogate data with similar characteristics to the analysed data ensures better evaluation of the false positive rate. False positive rate (*α*) is the probability of accepting coherence existence when the signals are actually not coupled. The threshold of the band coherence is determined as a function of the wavelet coefficient and specific frequency band to be analysed. The threshold level was analysed for the HF band (0.15–0.4) for *ω*_{0}: 6, 10, 15, 20 and 30. The false negative rate (*β*) assesses the probability of not rejecting the null hypothesis (accepting that there is no coupling), when the signals are actually coupled. In order to evaluate the error rates of the threshold level, simulations similar to those performed for valuating bias were used; realizations of the signals were analysed with different SNR values to impose the full scale of coherence levels. For each theoretical coherence level, the WTC band coherence values were compared with the zero threshold level. *β* was calculated by assessing the relative number of coherence estimations, which resulted in a value lower than the threshold level. The simulations were performed repeatedly and gave *N*>1000 realizations for each coherence level, in 0.1 steps. *β* was assessed for different wavelet coefficients: 6, 10, 15, 20 and 30.

### (c) Application to real cardiorespiratory data

#### (i) Patients and protocols

The application of the WTC was investigated in experimental data obtained from a group of ‘normal’ subjects. This study was approved by the Institutional Review Board of Tel Aviv University. All subjects signed a written informed consent form before participation. The group consisted of eight healthy male subjects (C1–C8), age 41±8 years. The protocol included 30 min rest in a supine position, followed by an active change of posture (CP) from supine to standing position (brief transition of 5 s) and 5 min in a standing position. Relevant signals extracted during these experiments were: electrocardiogram (ECG) leads I and II and respiratory signals from the ribs (Resp).

#### (ii) Data acquisition and preprocessing

The ECG and respiration were monitored non-invasively using a BIOPAC A/D system: the ECG was acquired with the BIOPAC's preamplifier and respiration with a Respitrace inductance plethysmograph. Respiration signals were low-pass filtered before analysis. The HR signal was obtained by detecting exact R-wave timing, by correcting any unrecognized arrhythmic beats, and by creating a continuous 10 Hz signal by interpolation.

#### (iii) Data analysis

We analysed two different physiological conditions, which are as follows.

Supine rest for 20 min, which represents a static situation without external stimulation. The first and last 5 min were omitted to ensure a physiological steady state.

CP from a supine to a standing position, which induces a strong response of the autonomic control system, exciting compensatory feedback mechanisms (Ewing

*et al.*1980; Wieling & Wesseling 1993). All the data in the time segments of interest were analysed in the time–frequency domain applying Morlet wavelets and the time–frequency coherence. The WTC procedure was applied to the signal pairs (HR–Resp). The analysis was focused on the HF band but extended between 0.02 and 1 Hz. Unless mentioned otherwise, all analyses were performed using*ω*_{0}=20, which was chosen based on results in §3*a,b*. The results were presented as means ±s.d. Significance calculations between different time regions were calculated using the Student*t*-test for two samples. The results were considered significant for*p*<0.05.

## 3. Results

### (a) Bias and s.d.

The bias and s.d. of the band WTC with regard to the theoretical band coherence levels were analysed for different wavelet coefficients. The bias is positive for low coherence levels, decreases almost linearly with increasing coherence and crosses zero for a small negative bias for all cases apart from *ω*_{0}=6 (figure 1*a*). The s.d. value increases as the coherence level decreases up to a maximum value. After the maximum point, the s.d. declines to a slightly lower value. As in the bias curve, the s.d. curves of lower coefficients have larger slope and higher maximum value (figure 1*b*).

### (b) Coherence threshold and detection error rate

The coherence threshold was 0.79, 0.64, 0.5, 0.45 and 0.3 with respect to *ω*_{0} of 6, 10, 15, 20 and 30. False negative error rates (*β*) are improved for larger *ω*_{0} (figure 1*c*). For *ω*_{0}=30, the error rate remains close to zero for a coherence level of as low as 0.55. Other slopes exhibit the same behaviour with lesser performance. For *ω*_{0}=6, the error rate rises sharply even for higher coherence levels. As expected, the probability extrapolates to 0.05 for zero coherence levels regardless of the value of the wavelet coefficients. This is in accordance with the critical level choice of 95 per cent for the WTC threshold estimation, which sets the false positive rate to be *α*=0.05.

### (c) Application to real cardiorespiratory data

Results from simulated signals have shown that higher coefficient values give better performance regarding bias and detection error rate. The Morlet coefficient *ω*_{0}=20 was subsequently used for the analysis of the human data to keep a high performance with a relatively high time resolution.

#### (i) Supine rest

Coherence between the cardiac and respiratory signals occurred mainly within the normative HF band (0.15–0.4 Hz) throughout the supine rest, as expected (figure 2). Although consistent, the coherence in the HF band is far from stability, with epochs of weak (weak power or thinner frequency bandwidth) or no coherence at all (below the threshold level) (figure 2*c*). The areas of inconsistency are usually due to a sudden change in breathing and/or HR (e.g. approx. 600 and 900 s). The epochs of inconsistency may vary in the range of 10 s up to 100 s. Not all visible changes in Resp or HR time series give a detectable change in the coherence map or in the HF band-averaged coherence. The subjects' group average value of the HF band coherence between HR and Resp, in the HF band, was 0.68±0.09. The group average of the band coherence s.d. (over time for each subject) was 0.13±0.03. One subject showed very low coherence (close to noise level), possibly because of an incorrect functioning of the respiration belt and was omitted from other analyses. Most of the subjects showed a weak but definite trace of another band just above 0.5 Hz (figure 2*a*), which is likely to be a second harmonic of the HF peak. The analysis of this coherence band was excluded from the scope of this study.

#### (ii) Change in posture

The dynamic trend of the coherence in the case of a time-defined physiological stimulus was demonstrated during CP (figure 3). Dynamics were best captured by observing the frequency-averaged WTC HF band. Since each subject performed the CP at slightly different timing, the zero time was determined at the onset of the CP, for each subject.

The band WTC of HR with RESP responds with a drop of coherence, which reaches its minimum value 20 s after CP (figure 4*a*) and recovers to a value significantly lower than that in the supine position (*p*<0.001). All subjects experienced a drop beyond the band threshold, which implies a loss of coherence in the HF band. The drop in coherence starts at approximately 20 s before the onset. Starting from a value of approximately 0.65 for the supine position, the WTC of HR–RESP stabilizes after the drop at approximately 0.55.

Examining the coherence power maps of the CP event, it was evident in all subjects, but one, that the loss of coherence in the HF band was accompanied by elevated coherence levels in the LF band. The coherence was mostly centred in the lower area of the LF band (0.04–0.07 Hz; figure 3*e*) and was consistent when averaged over subjects (figure 4*b*). The significance level for the lower area of the LF band was calculated using the same method as described in methods for the HF band, only using the appropriate frequency band. The coherence threshold was 0.46 for *ω*_{0} of 20.

## 4. Discussion

Time–frequency coherence by WTC was applied on simulated and human signals in order to quantify and validate its abilities to handle the dynamic analysis of RSA on the respiration band. The simulations focused on the HF band, which is more relevant to RSA analysis. The performance of the WTC estimator in other bands is expected to be similar (not shown), as can be seen by the close value between the threshold levels of the HF band and the lower LF band for *ω*_{0}=20.

Applying this approach, we analysed two different typical physiological situations. The WTC was used in three forms of presentation. (i) The full time and frequency-dependent coherence map, which gives a more qualitative yet comprehensive picture of the linear relationship between the two signals. Understanding this relationship can help in choosing optimal parameters for more qualitative tools (figure 2*a*). (ii) The time average WTC, which gives a steady-state representation, is similar to the classical Fourier-based coherence (figure 2*b*). (iii) The frequency average WTC, which gives a quantitative dynamic value of the coupling strength of two signals in a frequency band that usually reflects a specific control mechanism (figure 2*c*).

### (a) Bias and s.d.

We tested the WTC method as a function of the wavelet coefficient *ω*_{0}, which ‘calibrates’ the wavelet towards better frequency resolution for high *ω*_{0} levels and better time resolution for lower *ω*_{0} levels with obvious degradation on the opposite parameter. As expected from frequency analysis (Bendat & Piersol 1971), the bias and s.d. increase for smaller *ω*_{0} as well as for lower coherence levels. The small negative bias exhibited on higher coherence levels was less expected but has the same magnitude of s.d. The preferred wavelet coefficient will be the larger one (above 15), with no major differences between them.

### (b) Coherence threshold and detection error rate

The coherence threshold level acts as a marker for the existence of coupling. Although its absolute value is not relevant, as long as it is used with its current parameters (frequency band and *ω*_{0} in our case), higher threshold values do limit the ‘dynamic range’ of the estimates. Therefore, lower values should be preferred if possible, which means preference towards higher *ω*_{0} values. In the current work, we used a threshold level derived from simulated signals, which were fixed in time. A more dynamic approach, which takes into account the actual spectral characteristics of the analysed signals (e.g. Porta *et al.* 2002), should be considered in the future.

Error rates are higher for smaller *ω*_{0}, again owing to the averaging effect of longer (in time) wavelet function and longer smoothing.

### (c) Application to real cardiorespiratory data

*The first physiological condition* examined was half an hour of supine rest, which represents an ‘external steady-state’ condition. The Resp–HR WTC occurred, as expected, in the HF band representing RSA, and was strong and consistent overall. Nevertheless, almost all of the subjects exhibited areas of inconsistency in the coupling between the two signals even during an external steady state. Those inconsistencies were usually caused by a sudden change in the breathing rhythm and depth, and lasted for up to 100 s. The length of the epochs of the inconsistencies should be taken into account with caution, since WTC can have edge artefacts that can last up to (approx. 36 s in the HF band for *ω*_{0}=20).

For example, a deep sigh causes loss of coherence, which takes approximately 30 s to return to baseline. Taking into account the expected bias and *β* error rate of the coherence baseline levels exhibited during the supine state, it can be assumed that most of those below threshold inconsistencies represent real coherence loss. The demonstrated results suggest that the control of the ANS is more complex and that the assumption of the linearity is locally in time.

*The second physiological condition* examined was the CP, which is a time-defined external stimulus, with a known strong influence on the physiological state and the ANS regulatory systems. For all subjects, the CP event caused a general loss of coherence between HR and respiration. The loss of coherence signifies a full system readjustment, which is evidently a nonlinear process (Hainsworth & Mark 1993). Averaging in the HF band over all the subjects provides a clear insight into the dynamics of the CP event. The early response of the WTC can be explained as a cumulative effect of the window edge artefact and of an anticipatory response (the subjects were notified ahead of time about the coming change in posture).

Although not validated, we have explored and reported the evaluation of the above threshold coherence in the lower part of the LF band during CP. The origin of the elevated coherence exhibited should be cautiously interpreted, since a step followed by a decay was evident in the respiration signal baseline, which could be due to a belt artefact (figure 3*a*).

The coherence value after the transition to standing recovers to a significantly lower value (*p*<0.001) than in the initial supine position. The decrease in coherence power is in agreement with other steady-state methods, which quantified the RSA magnitude in supine and standing positions (Saul *et al.* 1991; Gilad *et al.* 2005).

## 5. Conclusion

The WT has already proven itself as a powerful tool in univariate analysis of dynamic situations in cardiovascular research. The present framework advances WTC as an intuitive and straightforward complementary tool for the bivariate spectral analysis of the ANS regulatory mechanisms. WTC could be applied as a direct tool for the analysis of the dynamic linear coupling between physiological signals, and as a marker for the ANS function. It can also be used for the detection of irregularities in the data when *a priori* assumptions of linear coupling are made, as well as for choosing accurate regions of linear coupling to improve the SNR in analyses such as transfer function or phase calculations.

## Acknowledgments

This work was partially supported by a grant from the Nicholas and Elizabeth Slezak Super Center for Cardiac Research and Biomedical Engineering at Tel Aviv University. The authors wish to thank Ori Gilad for discussion and review of this manuscript.

## Footnotes

↵† Deceased 2008.

One contribution of 15 to a Theme Issue ‘Addressing the complexity of cardiovascular regulation’.

- © 2009 The Royal Society