rhosa – Higher-Order Spectral Analysis in R

This package aims to provide functions to analyze and estimate higher-order spectra or polyspectra of multivariate time series, such as bispectrum and bicoherence [1].


You can install the released version of rhosa from CRAN with:


Alternatively, the development version from GitHub with remotes:

# install.packages("remotes")


This is an example, based on the outline at Figure 1 of [2], which shows how to use rhosa’s functions in the simplest way.

First of all, let’s load and attach rhosa.


With four consinusoidal waves having arbitrarily different phases (omega_a, omega_b, omega_c, and omega_d), but sharing a couple of frequencies (f_1 and f_2), we define two functions func_v and func_w to simulate stationary time series:

f_1 <- 0.1
f_2 <- 0.25
omega_a <- runif(1, min = 0, max = 2 * pi)
omega_b <- runif(1, min = 0, max = 2 * pi)
omega_c <- runif(1, min = 0, max = 2 * pi)
omega_d <- runif(1, min = 0, max = 2 * pi)
wave_a <- function(t) cos(2 * pi * f_1 * t + omega_a)
wave_b <- function(t) cos(2 * pi * f_2 * t + omega_b)
wave_c <- function(t) cos(2 * pi * f_1 * t + omega_c)
wave_d <- function(t) cos(2 * pi * f_2 * t + omega_d)
func_v <- function(t) wave_a(t) + wave_b(t) + wave_a(t) * wave_b(t)
func_w <- function(t) wave_c(t) + wave_d(t) + wave_c(t) * wave_b(t)

Both are oscillatory of course:

plot(func_v, 0, 100, n = 10001, col = "green", ylim = c(-3, 3), xlab = "t", ylab = "value")
plot(func_w, 0, 100, n = 10001, col = "orange", add = TRUE)

func_v and func_w.

func_v and func_w.

Then, obtain 10,000 samples from these time series in a fixed rate, with adding noises generated by the standard normal distribution; call them v and w.

raw_v <- sapply(1:10000, func_v)
raw_w <- sapply(1:10000, func_w)
v <- raw_v + rnorm(length(raw_v))
w <- raw_w + rnorm(length(raw_w))

Their first 100 samples looks like as follows:

plot(1:100, v[1:100], type = "l", col = "green", ylim = c(-3, 3), xlab = "t", ylab = "value")
lines(1:100, w[1:100], col = "orange")

v and w.

v and w.

It is noteworthy that these time series’ power spectrum densities are basically identical as shown in their spectral density estimation:

spectrum(v, main = "v", col = "green")
spectrum(w, main = "w", col = "orange")

Spectral density estimation from periodograms.Spectral density estimation from periodograms.

Spectral density estimation from periodograms.

On the other hand, we are going to see that their bispectra are different. Before estimating bispectrum, you may consider dividing samples into epochs or stretches of the same length for smoother results. rhosa’s bispectrum function accepts a matrix whose column represents an epoch, and returns a data frame as its result.

m_v <- matrix(v, ncol = 100)
m_w <- matrix(w, ncol = 100)
bs_v <- bispectrum(m_v)
bs_w <- bispectrum(m_w)

Because bispectrum is complex-valued at each frequency pair, you will find 3D plots are useful for visualizing results. Here we use rgl:

# install.packages("rgl")

… and define a couple of R functions to plot estimated bispectrum in the polar coordinates:

plot_bispectrum <- function(bs) {
    with(bs, {
        mfrow3d(nr = 1, nc = 2, sharedMouse = TRUE)
        plot3d(f1, f2, Mod(value), col = ifelse(Mod(value) > 1, "red", "blue"))
        grid3d(c("x", "y+", "z"))
        plot3d(f1, f2, Arg(value), col = ifelse(Mod(value) > 1, "red", "green"))
        grid3d(c("x", "y+", "z"))

Each line of the following snippets shows an interactive 3D point plot. The f1 and f2 axis represent normalized frequencies in unit cycles/sample of range [0, 1). While the modulus (or the absolute value) of the estimated complex values is plotted on the left side, the argument (or the angle) of them is on the right. A red point in the both sides indicates the higher absoute value than 1.


You must enable Javascript to view this page properly.


You must enable Javascript to view this page properly.

Unlike the power spectrum case, the bispectrum plots differ from each other, particularly the angle of v’s bispectrum at frequency pair (0.1, 0.25) is nearly 0 while w’s is far from 0.

Next, we would like to see if any quadratic phase coupling (QPC) can be detected from samples by identifying significant non-zero values of bicoherence. With rhosa’s bicoherence function, estimating (magnitude-squared) bicoherence is as easy as in the bispectrum case:

bc_v <- bicoherence(m_v, window_function = 'hamming')
bc_w <- bicoherence(m_w, window_function = 'hamming')

Note that in the above code the additional argument to bicoherence is given for tapering with Hamming window function.

plot_bicoherence <- function(bc) {
    with(bc, {
        plot3d(f1, f2, value, col = ifelse(significance, "red", "blue"))
        grid3d(c("x", "y+", "z"))

Frequency pairs of red points in the following plots indicate the existence of QPC adequately.


You must enable Javascript to view this page properly.


You must enable Javascript to view this page properly.

Please read help(bispectrum) or help(bicoherence) for more details.




The author thanks Alessandro E. P. Villa for his generous support to this project.


[1] Brillinger, D. R., & Irizarry, R. A. (1998). An investigation of the second- and higher-order spectra of music. Signal Processing, 65(2), 161–179. https://doi.org/10.1016/S0165-1684(97)00217-X

[2] Villa, A. E. P., & Tetko, I. V. (2010). Cross-frequency coupling in mesiotemporal EEG recordings of epileptic patients. Journal of Physiology-Paris, 104(3), 197–202. https://doi.org/10.1016/j.jphysparis.2009.11.024