1. Two dyads

Ten minutes of each partner, sampled at 2 Hz — a resolution an interpolated physiological series comfortably supports.

dyads |>
  filter(t <= 180) |>
  pivot_longer(c(A, B), names_to = "person", values_to = "value") |>
  ggplot(aes(t, value, colour = person)) +
  geom_line(linewidth = 0.4) +
  scale_colour_manual(values = PERSON) +
  facet_wrap(~ dyad, ncol = 1) +
  labs(title = "1. The two dyads",
       subtitle = "First three minutes, so that the faster rhythm is visible at all.",
       x = "time (s)", y = "value (SD units)", colour = NULL)

Both panels show two partners wandering in loose company. Nothing distinguishes them.

2. What a correlation reports

dyads |> summarise(`r over the whole recording` = round(cor(A, B), 3), .by = dyad)

## # A tibble: 2 × 2
##   dyad   `r over the whole recording`
##   <chr>                         <dbl>
## 1 Dyad 1                        0.443
## 2 Dyad 2                        0.328

Two moderate positive correlations. They differ a little, but neither says anything about where the coupling is, and on this evidence there is no reason to think the dyads differ in kind.

3. Where one person’s variance lives

Before asking what the partners share, it is worth asking what each of them has. A power spectrum decomposes one series’ variance by frequency: \(S_{AA}(f)\) is the amount of A’s variance attributable to oscillations at rate \(f\).

spectrum_of <- function(x) {
  sp <- spec.pgram(x, spans = c(9, 9), taper = 0.1, detrend = TRUE, plot = FALSE)
  tibble(freq = sp$freq * FS, power = sp$spec)
}

dyads |>
  group_split(dyad) |>
  map(\(d) bind_rows(spectrum_of(d$A) |> mutate(person = "A"),
                     spectrum_of(d$B) |> mutate(person = "B")) |>
        mutate(dyad = d$dyad[1])) |>
  list_rbind() |>
  filter(freq <= 0.6) |>
  ggplot(aes(freq, power, colour = person)) +
  geom_vline(xintercept = c(F_SLOW, F_FAST), colour = "grey85", linewidth = 0.3) +
  geom_line(linewidth = 0.5) +
  scale_colour_manual(values = PERSON) +
  scale_y_log10() +
  facet_wrap(~ dyad) +
  labs(title = "3. Each person's power spectrum",
       subtitle = paste0("All four people have power at both ", F_SLOW,
                         " Hz and ", F_FAST, " Hz."),
       x = "frequency (Hz)", y = "power (log scale)", colour = NULL)

All four look alike: two peaks apiece. A spectrum describes one person, so it cannot speak to the pair. For that the two series must be brought together.

4. The cross-spectrum, and coherence

The cross-spectrum \(S_{AB}(f)\) is the frequency-domain counterpart of a covariance. At each frequency it is a complex number whose modulus records how strongly the two series co-vary there, and whose argument records the phase offset between them. Dividing by each partner’s own power removes the dependence on amplitude and gives squared coherence:

\[C(f) \;=\; \frac{\lvert S_{AB}(f) \rvert^{2}}{S_{AA}(f)\, S_{BB}(f)} \qquad 0 \le C(f) \le 1\]

The expression is a squared correlation computed separately at every frequency, and it behaves like one: it has no sign, and it is bounded above by 1.

Taking the modulus discards the argument. One consequence is immediate and important: a constant delay between the partners changes the phase of \(S_{AB}\) but not its modulus, so coherence is unaffected by lag. Section 6 returns to the information this throws away.

coherence_of <- function(d, spans = c(9, 9)) {
  sp <- spec.pgram(cbind(d$A, d$B), spans = spans, taper = 0.1,
                   detrend = TRUE, plot = FALSE)
  tibble(freq = sp$freq * FS, coh = as.numeric(sp$coh), phase = sp$phase)
}

coh <- dyads |>
  group_split(dyad) |>
  map(\(d) coherence_of(d) |> mutate(dyad = d$dyad[1])) |>
  list_rbind()

coh |>
  filter(freq <= 0.6) |>
  ggplot(aes(freq, coh)) +
  geom_vline(xintercept = c(F_SLOW, F_FAST), colour = "grey85", linewidth = 0.3) +
  geom_line(linewidth = 0.6, colour = "grey20") +
  scale_y_continuous(limits = c(0, 1)) +
  facet_wrap(~ dyad) +
  labs(title = "4. Squared coherence",
       subtitle = "The two dyads that Section 2 could not tell apart.",
       x = "frequency (Hz)", y = "squared coherence")

in_band <- function(f, lo, hi) f >= lo & f <= hi

coh |>
  mutate(band = case_when(in_band(freq, 0.02, 0.06) ~ "slow (0.02-0.06 Hz)",
                          in_band(freq, 0.20, 0.30) ~ "fast (0.20-0.30 Hz)",
                          .default = NA_character_)) |>
  drop_na(band) |>
  summarise(`mean squared coherence` = round(mean(coh), 3), .by = c(dyad, band)) |>
  arrange(dyad, desc(band))

## # A tibble: 4 × 3
##   dyad   band                `mean squared coherence`
##   <chr>  <chr>                                  <dbl>
## 1 Dyad 1 slow (0.02-0.06 Hz)                    0.949
## 2 Dyad 1 fast (0.20-0.30 Hz)                    0.041
## 3 Dyad 2 slow (0.02-0.06 Hz)                    0.101
## 4 Dyad 2 fast (0.20-0.30 Hz)                    0.765

The dyads are opposites. Dyad 1 is coupled almost perfectly in the slow band and not at all in the fast one; dyad 2 is the reverse.

This is not a sharper version of the Section 2 result. It is information a correlation structurally cannot produce: a correlation sums over all frequencies before reporting, so slow agreement and fast disagreement are added together and the sum is middling. No choice of window or lag recovers the separation, because a window is a unit of time and the thing being separated is a frequency.

5. Smoothing is not optional

Coherence estimated from a single raw periodogram is exactly 1 at every frequency, for every pair of series, however unrelated. At one frequency the two complex numbers always lie on a line through the origin, so the normalisation in \(C(f)\) has nothing to shrink.

The estimate becomes meaningful only when \(S_{AB}\), \(S_{AA}\) and \(S_{BB}\) are each averaged over neighbouring frequencies before the ratio is taken. The spans argument sets the width of that average.

d1 <- filter(dyads, dyad == "Dyad 1")

list(NULL, c(3, 3), c(9, 9), c(25, 25)) |>
  map(\(sp) {
    z <- spec.pgram(cbind(d1$A, d1$B), spans = sp, taper = 0.1,
                    detrend = TRUE, plot = FALSE)
    f <- z$freq * FS
    tibble(spans = if (is.null(sp)) "none" else paste(sp, collapse = ", "),
           `slow band` = round(mean(z$coh[in_band(f, 0.02, 0.06)]), 3),
           `fast band` = round(mean(z$coh[in_band(f, 0.20, 0.30)]), 3),
           `off-band, 0.8-1.0 Hz` = round(mean(z$coh[in_band(f, 0.8, 1.0)]), 3))
  }) |>
  list_rbind()

## # A tibble: 4 × 4
##   spans  `slow band` `fast band` `off-band, 0.8-1.0 Hz`
##   <chr>        <dbl>       <dbl>                  <dbl>
## 1 none         1           1                      1    
## 2 3, 3         0.724       0.169                  0.227
## 3 9, 9         0.949       0.041                  0.073
## 4 25, 25       0.968       0.013                  0.015

Unsmoothed, all three columns read 1 — a perfect and perfectly uninformative answer. Widening the smoother sharpens the contrast, and does so by trading frequency resolution for stability: too little smoothing and everything is coherent, too much and genuinely distinct bands are blurred together.

Two practical consequences. A reported coherence near 1 with no smoothing parameter given should be read as a diagnostic of the analysis rather than as a result. And because the estimate is a function of the smoother, the smoothing must be reported alongside the coefficient.

6. Phase, and the timing that coherence discards

The argument of the cross-spectrum, \(\phi(f) = \arg S_{AB}(f)\), is the phase difference between the two series at frequency \(f\). A pure delay of \(\tau\) seconds shifts phase in proportion to frequency, so the delay can be read back:

\[\tau(f) \;=\; \frac{\phi(f)}{2\pi f}\]

To see this work, take dyad 1 again and let B’s copy of the shared slow rhythm arrive five seconds late.

read_delay <- function(d) {
  z <- coherence_of(d) |> filter(in_band(freq, 0.02, 0.06)) |> slice_max(coh, n = 1)
  tibble(`peak coherence` = round(z$coh, 3),
         `at frequency (Hz)` = round(z$freq, 3),
         `phase (rad)` = round(z$phase, 3),
         `implied delay (s)` = round(z$phase / (2 * pi * z$freq), 2))
}

bind_rows(`no delay` = read_delay(delayed_pair(0)),
          `5 s delay` = read_delay(delayed_pair(5)),
          .id = "simulation")

## # A tibble: 2 × 5
##   simulation `peak coherence` `at frequency (Hz)` `phase (rad)`[,1]
##   <chr>                 <dbl>               <dbl>             <dbl>
## 1 no delay              0.989               0.038            -0.034
## 2 5 s delay             0.977               0.047             1.40 
## # ℹ 1 more variable: `implied delay (s)` <dbl[,1]>

The delay is recovered, and the coherence barely moves. That is the claim of Section 4 made concrete: coherence measures whether a relationship exists at a frequency, and phase measures when. Reporting coherence alone discards the timing entirely.

Two cautions. The relation \(\tau = \phi / 2\pi f\) assumes a single constant delay; if the delay differs across frequencies, each frequency gives its own answer and there is no single \(\tau\). And phase is defined only modulo \(2\pi\), so a delay longer than one period of the frequency being examined is indistinguishable from a shorter one.

7. What coherence cannot see

Coherence integrates over the whole recording. Where a correlation reduces the time axis to one number, coherence reduces it to a curve — more informative, and still silent about when anything happened.

A third dyad makes the cost visible. Nothing about it is stated in advance.

coherence_of(dyad3) |>
  filter(freq <= 0.6) |>
  ggplot(aes(freq, coh)) +
  geom_vline(xintercept = c(F_SLOW, F_FAST), colour = "grey85", linewidth = 0.3) +
  geom_line(linewidth = 0.6, colour = "grey20") +
  scale_y_continuous(limits = c(0, 1)) +
  labs(title = "7. Dyad 3, coherence over the whole recording",
       subtitle = "Two modest bumps, against the near-ceiling single peaks of Figure 4.",
       x = "frequency (Hz)", y = "squared coherence")

coherence_of(dyad3) |>
  summarise(`slow band` = round(mean(coh[in_band(freq, 0.02, 0.06)]), 3),
            `fast band` = round(mean(coh[in_band(freq, 0.20, 0.30)]), 3))

## # A tibble: 1 × 2
##   `slow band` `fast band`
##         <dbl>       <dbl>
## 1       0.391       0.182

Two mediocre bumps, neither near the ceiling. The natural reading is a dyad weakly coupled in both bands. It is wrong, and no amount of re-smoothing will correct it, because the quantity being estimated is an average over a period in which the answer changed.

8. Wavelet coherence

The repair is to estimate coherence locally, using basis functions that are limited in time as well as in frequency. A Morlet wavelet is a wave packet — an oscillation inside a decaying envelope — and correlating both series against a family of these, at many scales and positions, yields coherence as a surface over time and frequency.

The smoothing of Section 5 is still required, and is now done in two directions at once: over neighbouring scales, and over neighbouring times.

wc <- analyze.coherency(
  as.data.frame(select(dyad3, A, B)),
  my.pair = c("A", "B"),
  dt = 1 / FS, dj = 1 / 12,
  lowerPeriod = 2, upperPeriod = 64,
  window.size.t = 100, window.size.s = 1 / 4,   # the smoothing
  make.pval = FALSE, verbose = FALSE
)

keep <- seq(1, NP, by = 4)   # thinned for plotting only

# as.numeric() flattens a matrix column-major, so `period` must vary fastest.
wtc <- expand_grid(idx = keep, period = wc$Period) |>
  mutate(freq = 1 / period,
         t = idx / FS,
         coh = as.numeric(wc$Coherence[, keep]),
         # Near either end a slow wavelet extends past the data; those
         # estimates are contaminated by the boundary.
         edge = t < 1.4 * period | t > SECS - 1.4 * period)

ggplot(wtc, aes(t, freq)) +
  geom_raster(aes(fill = coh), interpolate = TRUE) +
  geom_raster(data = \(d) filter(d, edge), fill = "white", alpha = 0.62) +
  geom_hline(yintercept = c(F_SLOW, F_FAST), colour = "grey30",
             linewidth = 0.3, linetype = "dotted") +
  geom_vline(xintercept = SECS / 2, colour = "grey10", linewidth = 0.4) +
  scale_fill_gradient(low = COH_LOW, high = COH_HIGH, limits = c(0, 1)) +
  scale_y_log10(breaks = c(0.02, 0.04, 0.1, 0.25, 0.5)) +
  labs(title = "8. Wavelet coherence for dyad 3",
       subtitle = paste0("Dotted lines mark ", F_SLOW, " and ", F_FAST,
                         " Hz; the vertical rule is at 300 s. ",
                         "Washed-out margins are the edge zone."),
       x = "time (s)", y = "frequency (Hz, log scale)", fill = "squared\ncoherence")

fr <- 1 / wc$Period
h1 <- 1:(NP / 2); h2 <- (NP / 2 + 1):NP
sl <- which(in_band(fr, 0.03, 0.055))
fa <- which(in_band(fr, 0.20, 0.30))

tibble(half = c("first 300 s", "second 300 s"),
       `slow band` = round(c(mean(wc$Coherence[sl, h1]), mean(wc$Coherence[sl, h2])), 3),
       `fast band` = round(c(mean(wc$Coherence[fa, h1]), mean(wc$Coherence[fa, h2])), 3))

## # A tibble: 2 × 3
##   half         `slow band` `fast band`
##   <chr>              <dbl>       <dbl>
## 1 first 300 s        0.92        0.167
## 2 second 300 s       0.626       0.839

Dyad 3 was coupled in the slow band for the first five minutes and in the fast band for the second. The whole-recording coherence of Section 7 averaged these into two unremarkable bumps; the surface separates them. In the fast band the change is clean, from 0.17 to 0.84.

9. The time–frequency trade

The slow band does not switch cleanly. It reads 0.92 in the first half and 0.63 in the second — the right direction, far from the separation the fast band shows.

This is not a defect of the implementation but a property of the problem. Estimating coherence at a frequency requires observing several cycles of it, and one cycle at 0.04 Hz lasts 25 seconds. Several of those occupy a large fraction of a 300-second half, so an estimate at that frequency cannot be localised to much better than a minute or two, and a change at the midpoint is necessarily smeared across it.

Time resolution and frequency resolution trade against each other, and the exchange rate is set by the frequency itself. A claim that coupling in a slow band began at a particular moment is, for this reason, rarely supportable from the data — whatever the surface appears to show. The washed-out margins of Figure 8 are the same effect at the boundaries: near the ends of the recording, a slow wavelet extends past the available data.