In this vignette we show how to detect quadratic phase coupling (QPC)
of one-dimensional or multi-dimensional real-valued time series by
`bicoherence`

or `cross_bicoherence`

of rhosa,
respectively. The first section gives an example applying
`bicoherence`

to the data from a simple model exhibiting QPC
at each frequency pair. In the second section we describe a generative
model of three channels with another kind of QPC revealed by
`cross_bicoherence`

. The third section summarizes when and
why `bicoherence`

or `cross_bicoherence`

fails to
recognize a certain type of QPC.

We begin with importing rhosa:

`library(rhosa)`

Our first mathematical model, adapted from [1], is a superposition of six cosine curves of unit amplitude with different frequencies, named \(x_1\):

\[x_1(t) = \sum_{i=1}^{6} \cos(\lambda_i t + \varphi_i)\]

where, for each \(i = 1,2,...,6\), \(\lambda_i\) is given and fixed, namely

\[ \lambda_1 = 0.55; \lambda_2 = 0.75; \lambda_3 = \lambda_1 + \lambda_2; \] \[ \lambda_4 = 0.6; \lambda_5 = 0.8; \lambda_6 = \lambda_4 + \lambda_5. \]

On the other hand, we choose \(\varphi_i\) (\(i = 1, ..., 5\)) independently from the uniform variable of range \([0, 2\pi)\), and define

\[ \varphi_6 = \varphi_4 + \varphi_5. \]

Note that the trigonometric identities implies

\[\cos(\lambda_6 t + \varphi_6) = \cos((\lambda_4 + \lambda_5) t + (\varphi_4 + \varphi_5)) = \cos(\lambda_4 t + \varphi_4) \cos(\lambda_5 t + \varphi_5) - \sin(\lambda_4 t + \varphi_4) \sin(\lambda_5 t + \varphi_5),\]

so \(\cos(\lambda_6 t + \varphi_6)\) is positively correlated with the product of \(\cos(\lambda_4 t + \varphi_4)\) and \(\cos(\lambda_5 t + \varphi_5)\). But \(\cos(\lambda_3 t + \varphi_3)\) is not correlated with the product of \(\cos(\lambda_1 t + \varphi_1)\) and \(\cos(\lambda_2 t + \varphi_2)\) in general as the phase \(\varphi_3\) is randomly assigned.

Once \(\varphi_i\)s are chosen,
\(x_1(t)\) is a periodic function of
\(t\). So it turns out that \(x_1\) is a (strictly) stationary stochastic
process. The wave length of any frequency component of \(x_1\) is shorter than \(4\pi\). Now consider sampling a realization
of \(x_1\) repeatedly during a fixed
length of time, say 256. The sampling rate is 1. That is, the interval
of consecutive samples is 1. The following R code effectively simulates
\(x_1\) as `x1`

.

```
<- function(a, b) c(a, b, a + b)
triple_lambda <- c(triple_lambda(0.55, 0.75), triple_lambda(0.6, 0.8))
lambda <- function(k) {
x1 set.seed(k)
<- runif(5, min = 0, max = 2*pi)
init_phi <- c(init_phi, init_phi[4] + init_phi[5])
phi function(t) do.call(sum, Map(function(l, p) cos(l * t + p), lambda, phi))
}<- function(f) {
observe sapply(seq_len(256), f)
}<- 100
N1 <- do.call(cbind, Map(observe, Map(x1, seq_len(N1)))) m1
```

Each column of matrix `m1`

in the last line represents a
realization of `x1`

. That is, we have taken 100 realizations
of \(x_1\) as `m1`

. Letâ€™s
plot them:

```
<- function(i) {
ith_sample data.frame(i = i, t = seq_len(256), v = m1[,i])
}<- do.call(rbind, Map(ith_sample, seq_len(100)))
r1
library(ggplot2)
ggplot(r1) +
geom_line(aes(t, v)) +
facet_wrap(vars(i)) +
theme(axis.title.y = element_blank(),
axis.text.y = element_blank(),
strip.background = element_blank(),
strip.text.x = element_blank())
```

A realization of \(x_1\) looks like:

`plot(m1[,1], type = "l", ylim = c(-5, 5), xlab = "t", ylab = "")`

Its spectrum estimation shows peaks at six frequencies as expected:

`spectrum(m1)`

Now we approximate \(x_1\)â€™s
bicoherence by `bicoherence`

:

`<- bicoherence(m1) bc1 `

â€¦ and define an R function to plot the heat map of the estimated bicoherence:

```
<- function(bc) {
heatmap_bicoherence ggplot(bc) +
geom_raster(aes(f1, f2, fill = value)) +
coord_fixed() +
scale_alpha(guide = "none")
}
heatmap_bicoherence(bc1)
```

Note that the highest peak is at the bifrequency \((f_1, f_2) = (\frac{\lambda_5}{2 \pi}, \frac{\lambda_4}{2 \pi}) \approx (0.127, 0.095)\).

Another example of QPC consists of three channels, which accept series of periodic input signals and suffer from Gaussian noises. From a couple of the channels, say \(C_1\) and \(C_2\), we observe the summation of input and noises as their output. On the other hand the last channel, called \(C_3\), adds \(C_1\)â€™s output multiplied by \(C_2\)â€™s to its own input and noise.

The following block diagram shows the skeleton of our three-channel model.Here, assume that \(C_1\)â€™s input is a triangle wave of a fixed frequency with varying initial phases. A rectangle wave of another frequency for \(C_2\)â€™s input. A cosinusoidal curve of yet another frequency for \(C_3\)â€™s input. Running the following code simulates the model.

```
<- 1.2
Fcoef1 <- 0.7
Fcoef2 <- 0.8
Fcoef3
<- function(x, p) {2 * asin(sin(Fcoef1 * x + p))}
i1 <- function(x, p) {ifelse(cos(Fcoef2 * x + p) >= 0, -1, 1)}
i2 <- function(x, p) {cos(Fcoef3 * x + p)}
i3
<- 0.3
Qcoef
<- function(k) {
tc set.seed(k)
<- runif(3, min = 0, max = 2*pi)
ps function(x) {
<- i1(x, ps[1]) + rnorm(length(x), mean = 0, sd = 1)
c1 <- i2(x, ps[2]) + rnorm(length(x), mean = 0, sd = 1)
c2 <- Qcoef * c1 * c2 +
c3 i3(x, ps[3]) + rnorm(length(x), mean = 0, sd = 1)
data.frame(c1, c2, c3)
}
}
<- 100
N2
<- function() {
sample_tc Map(function(f) {f(seq_len(256))}, Map(tc, seq_len(N2)))
}
<- function(y) {
c1_data_frame do.call(cbind, Map(function(k) {y[[k]]$c1}, seq_len(N2)))
}
<- function(y) {
c2_data_frame do.call(cbind, Map(function(k) {y[[k]]$c2}, seq_len(N2)))
}
<- function(y) {
c3_data_frame do.call(cbind, Map(function(k) {y[[k]]$c3}, seq_len(N2)))
}
<- sample_tc()
y1 <- c1_data_frame(y1)
d1 <- c2_data_frame(y1)
d2 <- c3_data_frame(y1) d3
```

That is, we obtain 100 series of data with 256 points. To be specific, \(C_1\)â€™s triangle wave has cycle \(\frac{2 \pi}{1.2}\), the cycle of \(C_2\)â€™s rectangle wave is \(\frac{2 \pi}{0.7}\), and the one of \(C_3\)â€™s cosinusoidal wave is \(\frac{2 \pi}{0.8}\).

\(C_1\)â€™s output looks like:

`plot(d1[,1], type = "l", ylim = c(-5, 5), xlab = "t", ylab = "")`

And \(C_2\)â€™s output:

`plot(d2[,1], type = "l", ylim = c(-5, 5), xlab = "t", ylab = "")`

\(C_3\)â€™s output:

`plot(d3[,1], type = "l", ylim = c(-5, 5), xlab = "t", ylab = "")`

Their spectrum estimation show significant peaks at expected frequencies:

`spectrum(d1)`

`spectrum(d2)`

`spectrum(d3)`