A toy model of QPC
We begin with importing 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.
triple_lambda <- function(a, b) c(a, b, a + b)
lambda <- c(triple_lambda(0.55, 0.75), triple_lambda(0.6, 0.8))
x1 <- function(k) {
set.seed(k)
init_phi <- runif(5, min = 0, max = 2*pi)
phi <- c(init_phi, init_phi[4] + init_phi[5])
function(t) do.call(sum, Map(function(l, p) cos(l * t + p), lambda, phi))
}
observe <- function(f) {
sapply(seq_len(256), f)
}
N1 <- 100
m1 <- do.call(cbind, Map(observe, Map(x1, seq_len(N1))))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:
ith_sample <- function(i) {
data.frame(i = i, t = seq_len(256), v = m1[,i])
}
r1 <- do.call(rbind, Map(ith_sample, seq_len(100)))
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())
x1.A realization of \(x_1\) looks like:
x1.Its spectrum estimation shows peaks at six frequencies as expected:
x1.Now we approximate \(x_1\)’s
bicoherence by bicoherence:
… and define an R function to plot the heat map of the estimated bicoherence:
heatmap_bicoherence <- function(bc) {
ggplot(bc) +
geom_raster(aes(f1, f2, fill = value)) +
coord_fixed() +
scale_alpha(guide = "none")
}
heatmap_bicoherence(bc1)
x1’s estimated bicoherence.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)\).
