# coherence.jl

*FourierAnalysis* estimates coherence both in the *frequency domain*, in this unit, and in the *time-frequency domain*, in unit timefrequencybi.jl.

### Notation

The following notation will be used:

symbol | meaning | Julia function |
---|---|---|

superscript $^*$ | complex conjugate | conj() |

superscript $^H$ | complex conjugate transpose | ' |

$\mid \cdot \mid$ | modulus | abs() |

$\Bbb R(\cdot)$ | real part of the argument | real() |

$\Bbb C(\cdot)$ | imaginary part of the argument | imag() |

$\left<\cdot\right>$ | an average across argument's elements | mean() |

## coherence estimates

**frequency domain (FD)**

In the FD, the input is a continuous multivariate data matrix $X=[x_1 ... x_n]$, where each column is a series (e.g., a time-series). Coherence estimates take the form of a $n$x$n$ coherence matrix $C(f)$ for each discrete Fourier frequency $f$. Those are normalized versions of the $n$x$n$ *cross-spectral matrices* $S(f)$ at corresponding frequencies.

Cross-spectral matrices are estimated via a Welch procedure, that is, they are averaged across windows sliding over $X$. The windows may be overlapping or not. Let $Y_l(f)$ be the discrete Fourier transform at discrete Fourier frequency $f$ of the $l^{th}$ window of $X$. The cross-spectral matrix estimation at frequency $f$ is then given by

$S(f)=\left<Y(f)Y(f)^H\right>$,

where the average is taken across a number of sliding windows.

Therefore, coherence matrices in the FD is a measure of the *synchronization between all $n$ series of $X$ taken pair-wise*, resulting in a symmetric matrix for each discrete Fourier frequency.

**time-frequency domain (TFD)**

There exist many time-frequency representation, *e.g.*, those obtained using wavelets, short-time Fourier transform, etc. In *FourierAnalysis* time-frequency representations are obtained passing a univariate signal through a filter bank and computing the analytic signal of the filter output via the Hilbert transform. This results in a matrix with the band-pass regions of the filer bank in rows and the $t$ time samples in columns (see `filterbank`

and `analyticsignal`

for details).

For estimating coherence in the TFD, the input is composed of two sets of $k$ univariate data vectors $π±_a=\{x_{a1}, x_{a2},...,x_{ak}\}$ and $π±_b=\{x_{b1}, x_{b2},...,x_{bk}\}$, where each element of the set is to be understood as a single window holding $t$ samples and where the elements in corresponding positions forms pairs. Let $Y_r(a)$ and $Y_r(b)$ be the time-frequency analytic signal of pairs $x_{ar}$, $x_{br}$, for all $r=1...k$. Coherence estimates also have a time-frequency representation. They are a normalized version of crossed analytic signals, which are analogous to cross-spectra and in a time-frequency plane are given by

$Z=\left<Y(a) \odot Y(b)^*\right>$,

where the average is taken across the $k$ pairs (realizations).

Therefore, the coherence matrix in the TFD is a measure of the *synchronization between pairs $π±_a$, $π±_b$* for each point in the time-frequency plane.

### kinds of coherence

*FourierAnalysis* estimates several kinds of *linear* and *non-linear* squared coherences, simply referred to as *coherence*. Let $i$ and $j$ indicate two time series, $\left<s_{ij}\right>$ be an average cross-spectrum between $i$ and $j$, and $\left<s_{i}\right>$, $\left<s_{j}\right>$ be the auto-spectrum of $i$ , $j$. In the frequency domain those are function of frequency, whereas in the time-frequency domain they are functions of both time and frequency.

Finally, for time-frequency data let $w$ denote non-negative weights normalized so that their average is 1. Those are weights for the pairs on which the averages are computed. Setting all weights equal to 1, gives the unweighted version of all measures, which is the only supported option in the FD, since therein the average is taken across windows and since windows segmentation is arbitrary, weighting is usually meaningless.

All kinds of coherences estimated in *FourierAnalysis* are summarized in the following table:

kind | linear | non-linear |
---|---|---|

total | $\left<w\mid s_{ij} \mid^2\right>\big /\left<ws_is_j\right>$ | $\left<w\mid s_{ij} \mid^2\right>$ |

real | $\left<w\mid \Bbb R(s_{ij}) \mid^2\right>\big/\left<ws_is_j\right>$ | $\left<w\mid \Bbb R(s_{ij})\mid^2\right>$ |

imaginary | $\left<w\mid \Bbb C(s_{ij}) \mid^2\right>\big/\left<ws_is_j\right>$ | $\left<w\mid \Bbb C(s_{ij}) \mid^2\right>$ |

instantaneous | $\left<w\mid \Bbb R(s_{ij}) \mid^2\right>\big/\left<w(s_is_j-\mid \Bbb C(s_{ij}) \mid^2)\right>$ | $\left<w\mid \Bbb R(s_{ij}) \mid^2\right>\big/\left<w(1-\mid \Bbb C(s_{ij}) \mid^2)\right>$ |

lagged | $\left<w\mid \Bbb C(s_{ij}) \mid^2\right>\big/\left<w(s_is_j-\mid \Bbb R(s_{ij}) \mid^2)\right>$ | $\left<w\mid \Bbb C(s_{ij}) \mid^2\right>\big/\left<w(1-\mid \Bbb R(s_{ij}) \mid^2)\right>$ |

The *linear* **total** coherence is the classical *squared coherence* measure and its *non-linear* counterpart is known as *phase-locking value* or *phase coherence*.

The **real** coherence (Pascual-Marqui, 2007) and the **imaginary** coherence (Nolte *et al.*, 2004), sum up to the total coherence.

The **lagged** coherence has been proposed by Pascual-Marqui (2007). As a completion of the table, we here also define the **instantaneous** coherence in an analogous way.

Corresponding *linear* and *non-linear* measures are the same, but are computed on linear and non-linear cross-spectra in the FD and crossed analytic-signal in the TFD, respectively. In fact, the non-linear expressions are obtained setting $s_i=s_j=1$, which is the case for non-linear cross-spectra and analytic signals. See CrossSpectra and TFAnalyticSignal. For further explanations on the linearity of those objects, see Congedo (2018).

### interpreting coherence

The real part of the cross-spectrum, named the *co-spectrum*, describes the *synchronous* synchronization, that is, in-phase synchronization or out-of-phase synchronization. The imaginary part, named the *quadrature spectrum*, describes the *asynchronous* synchronization, that is, with a quarter of a cycle lead or lag (see Congedo, 2018). Consequently,

- the
*total*coherence is sensitive to all kinds of synchronization, - the
*real*coherence is sensitive only to synchronous synchronization, - the
*instantaneous*coherence is sensitive to all but asynchronous synchronization, - the
*imaginary*coherence is sensitive only to asynchronous synchronization, - the
*lagged*coherence is sensitive to all but synchronous synchronization.

## Coherence

**Categories**: data objects, FDobjects

In the time-frequency domain coherence estimates are given as Julia $Matrix$ objects. In the frequency-domain they are encapsulated in the following structure:

```
struct Coherence
<same fields of the CrossSpectra structure>
end
```

This object has the same structure of the CrossSpectra object, with the difference that its `y`

data field holds real matrices, whereas for cross-spectra holds complex matrices.

By convention the diagonal elements of all FD coherence matrices are filled with 1.0. Note that the FFT coefficients for the DC level and the Nyquist frequency are real, therefore for those frequencies the identity matrix is returned for the *imaginary* and *lagged* coherence. Also, note that at these frequencies the *instantaneous* coherence is equal to the *real* coherence. In the time-frequency domain, we do not have discrete frequencies but a filter bank, thus no special behavior appears at the edges.

A vector of *Coherence* objects is of type CoherenceVector.

**Methods for Coherence and CoherenceVector objects**

method | Coherence | CoherenceVector |
---|---|---|

`bands` | β | β |

`extract` | β | β |

`mean` | β | β |

`smooth` | β | β |

`sameParams` | β |

**Generic Constructors**:

In order to construct a *Coherence* object in the frequency domain from multivariate data, *FourierAnalysis* provides two `coherence`

constuctors from raw data and two constuctors building coherence matrices from CrossSpectra and CrossSpectraVector objects. Those are what you will use in practice most of the time.

Manual constructors are also possible, for which you have to provide appropriate arguments. Generic constructors follow exactly the same syntax as the generic constructors of CrossSpectra objects.

In order to estimate coherence in the time-frequency domain from sets of univariate data, *FourierAnalysis* provides two more `coherence`

functions.

**Constructors from data**:

`FourierAnalysis.coherence`

β Function```
(1)
function coherence(π :: CrossSpectra;
allkinds :: Bool = false)
(2)
function coherence(π’ :: CrossSpectraVector;
allkinds :: Bool = false,
check :: Bool = true)
(3)
function coherence(X :: Matrix{T},
sr :: Int,
wl :: Int;
tapering :: Union{Taper, TaperKind} = harris4,
planner :: Planner = getplanner,
slide :: Int = 0,
DC :: Bool = false,
nonlinear :: Bool = false,
smoothing :: Smoother = noSmoother,
tril :: Bool = false,
allkinds :: Bool = false,
β© :: Bool = true) where T<:Real
(4)
function coherence(π :: Vector{Matrix{T}},
< same argument sr, ..., β© of method (3) > where T<:Real
```

**alias**: coh

(1)

Construct a Coherence object from a CrossSpectra object, allowing coherence estimates using the Welch method. All non-data fields are copied from the cross-spectra object, i.e., all fields but `y`

, which holds the coherence matrices that are computed by this function.

If `π.tril`

is true, only the lower triangular part of the coherence matrices is computed and the matrices in `y`

are LowerTriangular matrices, otherwise the full matrices are computed and `y`

will hold Hermitian matrices (see LinearAlgebra).

If optional keyword argument `allkinds`

is true, all five kinds of coherence are returned. In this case the output is a 5-tuple of Coherence objects, in the order:

*total*coherence,*real*coherence,*instantaneous*coherence*imaginary*coherence,*lagged*coherence.

If `allkinds`

is false (default), only the *total* (classical) coherence is returned as a single `Coherence`

object.

(2)

Construct a CoherenceVector object from the input CrossSpectraVector object `π’`

, allowing coherence estimates using the Welch method. This method calls method (1) for all objects in `π’`

.

The `allkinds`

optional keyword parameter has the same meaning as in method (1). In this method if the argument is passed as true, the output is a 5-tuple of `CoherenceVector`

objects.

If optional keyword argument `check`

is true (default), it is checked that all `CrossSpectra`

objects in `π’`

have the same value in the ``.sr`

, `.wl`

, `.DC`

, `.taper`

, `.nonlinear`

and `.smoothing`

field. If this is not the case, an error message is printed pointing to the first field that is not identical in all objects and `Nothing`

is returned.

(3)

Given a multivariate time series data matrix `X`

of dimension $t$ x $n$, where $t$ is the number of samples (rows) and $n$ the number of time-series (columns), sampling rate `sr`

and epoch length `wl`

, compute the squared coherence matrices of `X`

, that is, the coherence matrices at all Fourier discrete frequencies obtained from the Welch (sliding window) average cross-spectra.

**optioal keyword arguments**:

`sr`

, `wl`

, `tapering`

, `planner`

and `slide`

have the same meaning as for the `spectra`

function.

`DC`

, `nonlinear`

, `smoothing`

, `tril`

and `β©`

have the same meaning as for the `crossSpectra`

function, to which they apply for estimating the cross-spectra.

The `allkinds`

optional keyword parameter has the same meaning as in method (1).

(4)

Construct a CoherenceVector object from a vector of real multivariate data matrices. Compute the coherence matrices from the cross-spectral matrices estimated using the Welch method as per method (3) for all $k$ data matrices in `π`

.

The $k$ matrices in `π`

must have the same number of columns (*i.e.*, the same number of time-series), but may have any number of (at least `wl`

) rows (samples). All other arguments have the same meaning as in method (3), with the following difference:

`β©`

: if true (default), the method is run in multi-threaded mode across the $k$ data matrices if $k$ is at least twice the number of threads Julia is instructed to use, otherwise this method attempts to run each coherence estimation in multi-threaded mode across series as per method (3). See Threads.

If a `Planner`

is not explicitly passed as an argument, the FFTW plan is computed once and applied for all coherence estimations.

**note for methods (1) and (3)**:

If `tril`

is false (default), the output coherence object(s) is of type `Array{Hermitian,1}`

, which is the `βVector`

type used in package PosDefManifold. Since coherence estimates are symmetric positive definite, they can be straightaway used as argument to PosDefManifold's functions, e.g., for computing matrix moves on geodesics, matrix distances, etc. and the the whole vector output to compute matrix means, spectral embedding and more.

If `tril`

is true, the output is of type `Array{LowerTriangular,1}`

, which is the `πVector`

type used in PosDefManifold, that is, only the lower triangle of the coherencee is computed in order to save time and memory.

**note for methods (2) and (4)**:

If `tril`

is false, the output is of type `Array{Array{Hermitian,1},1}`

, which is the `βVectorβ`

type used in PosDefManifold.

If `tril`

is true, the output is of type `Array{Array{LowerTriangular,1},1}`

, which is the `πVectorβ`

type used in PosDefManifold.

**See**: crossspectra.jl, Spectra, Coherence.

**Examples**:

```
## common code for methods (1)-(4)
using FourierAnalysis, LinearAlgebra
function generateSomeData(sr::Int, t::Int; noise::Real=1.)
# four sinusoids of length t samples and sr sampling rate
# peak amplitude: 0.7, 0.6, 0.5, 0.4
# frequency: 5, 7, 13, 27
# phase: 0, Ο/4, Ο/2, Ο
v1=sinusoidal(0.7, 5, sr, t, 0)
v2=sinusoidal(0.6, 7, sr, t, Ο/4)
v3=sinusoidal(0.5, 13, sr, t, Ο/2)
v4=sinusoidal(0.4, 27, sr, t, Ο)
return hcat(v1, v2, v3, v4) + (randn(t, 4)*noise)
end
sr, wl, t = 128, 512, 8192
# (1)
X=generateSomeData(sr, t) # multivariate data matrix 8192x4
# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)
# Only classical coherence
C=coherence(S)
# All 5 kinds of coherence
Ctot, C2real, C3inst, C4imag, C5lag=coherence(S, allkinds=true);
Ctot
# (2)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# cross-spectra using default harris4 tapering window
S=crossSpectra(X, sr, wl)
# Only classical coherence
C=coherence(S)
# All 5 kinds of coherence
Ctot, C2real, C3inst, C4imag, C5lag=coherence(S, allkinds=true);
Ctot
# (3)
X=generateSomeData(sr, t) # multivariate data matrix 8192x4
# coherence using default harris4 tapering window
C=coherence(X, sr, wl)
# check the coherence matrix at frequency 5Hz
C.y[f2b(5, sr, wl)]
# coherence using hann tapering window
C=coherence(X, sr, wl; tapering=hann)
# using Slepian's multi-tapering
C=coherence(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear coherence (phase-locking value)
C=coherence(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of coherence matrices
C=coherence(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# compute all five kinds of coherence
Ctot, Creal, Cinst, Cimag, Clag=coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, allkinds=true);
Ctot
# smooth a-posteriori the coherence
C2=smooth(blackmanSmoother, C)
# or compute coherence already smoothed
C=coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, smoothing=blackmanSmoother)
# mean coherence matrix in 8Hz-12Hz range
M=mean(C, (8, 12)) # or also M=mean(C, (8.0, 12.0))
# extract all coherence matrices in 8Hz-12Hz range
E=extract(C, (8, 12))
# coherence matrices averaged in 2Hz band-pass regions
B=bands(C, 2)
# (4)
# generate 3 multivariate data matrices 8192x4
X=[generateSomeData(sr, t) for i=1:3]
# The examples here below use exactly the same syntax as method (3).
# However, since the input here is a vector of data matrices
# and not a single data matrix, the examples here below create a vector
# of the object created by the examples of method (3).
# For example:
# coherence using the default harris4 tapering window
# this creates a CoherenceVector object
C=coherence(X, sr, wl)
# check the first Coherence object
C[1]
# check the coherence matrix at fr. 5Hz for the first Coherence object
C[1].y[f2b(5, sr, wl)]
# coherence using Hann tapering window
C=coherence(X, sr, wl; tapering=hann)
# using Slepian's multi-tapering
C=coherence(X, sr, wl; tapering=slepians(sr, wl))
# compute non-linear coherence
C=coherence(X, sr, wl; tapering=slepians(sr, wl), nonlinear=true)
# compute only the lower triangle of coherence matrices
C=coherence(X, sr, wl; tapering=slepians(sr, wl), tril=true)
# compute all five kinds of coherence
Ctot, Creal, Cinst, Cimag, Clag=coherence(X, sr, wl;
tapering=slepians(sr, wl), tril=true, allkinds=true);
Ctot
# smooth a-posteriori all coherence objects in S
C2=smooth(blackmanSmoother, C)
# or compute them all already smoothed
C=coherence(X, sr, wl; tapering=parzen, smoothing=blackmanSmoother)
# mean coherence matrix in 8Hz-12Hz range for all coherence objects
M=mean(C, (8, 12)) # or also M=mean(C, (8.0, 12.0))
# grand-average mean of the above across all Coherence objects
meanM=mean(mean(C, (8, 12)))
# extract all coherence matrices in 8Hz-12Hz range for all coherence objects
E=extract(C, (8, 12))
# grand average of coherence matrices in 8Hz-12Hz range for all coh. objects
meanE=mean(extract(C, (8, 12)))
# coherence matrices averaged in 2Hz band-pass regions for all coh. objects
B=bands(C, 2)
# Pre-compute a FFT planner and pass it as argument
# (this interesting if the function is to be called repeatedly).
plan=Planner(plan_exhaustive, 10.0, wl, eltype(X[1])) # wait 10s
C=coherence(X, sr, wl; planner=plan)
# how faster is this?
using BenchmarkTools
@benchmark(coherence(X, sr, wl))
@benchmark(coherence(X, sr, wl; planner=plan))
...
...
```

```
(5)
function coherence(πβ :: TFAnalyticSignalVector,
πβ :: TFAnalyticSignalVector,
frange :: fInterval,
trange :: tInterval;
nonlinear :: Bool = false,
allkinds :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
check :: Bool = true)
(6)
function coherence(π±β :: Vector{Vector{T}},
π±β :: Vector{Vector{T}},
sr :: Int,
wl :: Int,
frange :: fInterval,
trange :: tInterval,
bandwidth :: IntOrReal = 2;
nonlinear :: Bool = false,
allkinds :: Bool = false,
mode :: Function = extract,
func :: Function = identity,
w :: Vector = [],
filtkind :: FilterDesign = Butterworth(2),
fmin :: IntOrReal = bandwidth,
fmax :: IntOrReal = srΓ·2,
fsmoothing :: Smoother = noSmoother,
tsmoothing :: Smoother = noSmoother,
planner :: Planner = getplanner,
β© :: Bool = true) where T<:Real
```

**alias**: coh

If optional keyword parameter `nonlinear`

is false (default), estimate linear (weighted) coherence measure, otherwise estimate the the (weighted) phase concentration measure.

If optional keyword argument `allkinds`

is true all five kinds of coherence are returned. In this case the output is a 5-tuple of `Coherence`

matrices, in the order:

*total*coherence,*real*coherence,*instantaneous*coherence*imaginary*coherence,*lagged*coherence.

If `allkinds`

is false (default) only the *total* (classical) coherence is returned as a single `Coherence`

matrix.

(5) The desired measure is obtained averaging across the TFAnalyticSignal objects in `π`

. Since this method uses pre-computed analytic signal objects, their `.nonlinear`

field must agree with the `nonlinear`

argument passed to this function.

`frange`

, `trange`

, `w`

, `mode`

and `func`

have the same meaning as in the `meanAmplitude`

function, however keep in mind that the two possible `mode`

functions, i.e., `extract`

and `mean`

, in this function operate on complex numbers.

The checks performed in the `meanAmplitude`

function are performed here too. In addition, if `check`

is true, also check that

- if
`nonlinear`

is true, all objects in`π`

are nonlinear; - if
`nonlinear`

is false, all objects in`π`

are linear.

If either check fails, print an error message and return `Nothing`

.

(6) Estimate the analytic signal of all data vectors in `π±`

calling the `TFanalyticsignal`

constructor and then use method (5) to obtained the desired measure.

`frange`

, `trange`

, `mode`

, `func`

, `w`

and `check`

have the same meaning as in the `meanAmplitude`

function. The other arguments are passed to the `TFanalyticsignal`

constructor, to which the reader is referred for understanding their action.

**See also**: `meanAmplitude`

, `meanDirection`

, timefrequencybi.jl.

**Examples**: see the examples of `comodulation`

function.

**References**:

M. Congedo (2018), Non-Parametric Synchronization Measures used in EEG and MEG, Technical Report. GIPSA-lab, CNRS, University Grenoble Alpes, Grenoble INP.

G. Nolte *et al.* (2004), Identifying true brain interaction from EEG data using the imaginary part of coherency. Clin Neurophysiol, 115(10), 2292-307.

R. Pascual-Marqui (2007), Instantaneous and lagged measurements of linear and nonlinear dependence between groups of multivariate time series: frequency decomposition, arXiv:0711.1455.