A method for continuously assessing the autonomic response to music-induced emotions through HRV analysis

Abstract

Interest in therapeutic applications of music has recently increased, as well as the effort to understand the relationship between music features and physiological patterns. In this study, we present a methodology for characterizing music-induced effects on the dynamics of the heart rate modulation. It consists of three steps: (i) the smoothed pseudo Wigner-Ville distribution is performed to obtain a time–frequency representation of HRV; (ii) a parametric decomposition is used to robustly estimate the time-course of spectral parameters; and (iii) statistical population analysis is used to continuously assess whether different acoustic stimuli provoke different dynamic responses. Seventy-five healthy subjects were repetitively exposed to pleasant music, sequences of Shepard tones with the same tempo as the pleasant music and unpleasant sounds overlaid with the same sequences of Shepard tones. Results show that the modification of HRV parameters are characterized by an early fast transient phase (15–20 s), followed by an almost stationary period. All kinds of stimuli provoked significant changes compared to the resting condition, while during listening to pleasant music the heart and respiratory rates were higher (for more than 80% of the duration of the stimuli, p < 10⁻⁵) and the power of the HF modulation was lower (for more than 70% of the duration of the stimuli, p < 0.05) than during listening to unpleasant stimuli.

Appendix

Consider a signal component x(t), whose ACF is the damped complex sinusoid y(τ):

$$ y(\tau) = C e^{-d|\tau| + j2\pi \nu_0\tau} $$

(10)

with d > 0, τ and $C\in{\mathbb{R}}.$ The Fourier transform of y(τ) is equal to:

$$ Y(\nu) = \int_{-\infty}^{\infty} C e^{-d|\tau| + j2\pi (\nu_0-\nu)\tau}d\tau = \frac{2Cd}{4\pi^{2}(\nu-\nu_{0})^{{2}} + {\hbox{d}}^2}.$$

(11)

Note that y(τ) is an hermitian function and Y(ν) is real. The function Y(ν) is the power spectral density of x(t) and has a Lorentian shape with a peak centered on frequency ν₀. The power of the signal component x(t), P _x , is then closely related to coefficient C and it can be analytically obtained as the total area of Y(ν):

$$\begin{aligned} P_x&= \int_{-\infty}^{\infty}Y(\nu){\hbox{d}}\nu = \int_{-\infty}^{\infty}\left[\frac{2Cd}{4\pi^2(\nu-\nu_{0})^{\hbox{ 2}} + d^{\hbox{ 2}}}\right]{\hbox{d}}\nu\\ &= \frac{2Cd}{4\pi^2}\int_{-\infty}^{\infty}\left[\frac{1}{(\nu-\nu_{0})^{\hbox{ 2}} + \left(\frac{d}{2\pi}\right)^2}\right]{\hbox{d}}\nu\\ &= \frac{2Cd}{4\pi^2}\frac{2\pi}{d} \left[ \arctan\left(\frac{(\nu-\nu_{0})2\pi}{{\hbox{d}}} \right) \right]_{-\infty}^{\infty} = C \end{aligned} $$

(12)