NeuroAnalyzer tutorials: Analyze: Functional connectivity: Coherence
Initialize NeuroAnalyzer
using NeuroAnalyzer
using DSP
using Plots
using Statistics
eeg = load("files/eeg.hdf")
e10 = epoch(eeg, ep_len=10);Cross Power Spectrum Density (CPSD)
Coherence is calculated from Cross Power Spectrum Density (CPSD).
CPSD is a measure used in signal processing to analyze the relationship between two signals in the frequency domain. It shows how the power of one signal relates to another at different frequencies.
\[ P_{xy} = X × Y^* \] \(P_{xy}\): CPSD between two signals, \(x\) and \(y\)
\(X\): Fourier transform of \(x\)
\(Y^*\): complex conjugate of the Fourier transform of \(y\)
The magnitude of \(P_{xy}\) indicates the strength of the relationship between \(x\) and \(y\) at each frequency.
The phase of \(P_{xy}\) gives information about the phase difference between the signals at each frequency.
Calculate CPSD:
cp, f = cpsd(e10,
e10,
ch1="Fp1",
ch2="Fp2",
ep1=10,
ep2=10,
method=:mt)(pxy = ComplexF64[437.0236942544301 + 0.0im 868.9236051592155 - 12.508792647975682im … -8.276320194032056e-7 + 5.675676411344001e-7im -3.4888346683607675e-7 + 0.0im;;;],
f = [0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, 0.5, 0.5625 … 127.4375, 127.5, 127.5625, 127.625, 127.6875, 127.75, 127.8125, 127.875, 127.9375, 128.0],)Available methods: :mt (multi-taper), :fft (Fast Fourier transform), :stft (short time Fourier transform).
Plot magnitudes:
plot_psd(f,
abs.(cp[1, :, 1]),
xlabel="Frequency [Hz]",
ylabel="Magnitude")Magnitude peaks: indicate strong connectivity at specific frequencies.
Plot phases:
plot_phase(DSP.angle.(cp[1, :, 1]),
f,
ylims=(-pi, pi))Phase information: shows if one region leads or lags another, which can imply directionality.
CPSD may be also compared between conditions (e.g. rest vs. task, healthy vs. patient) at the same location.
Coherence
Both coherence and magnitude squared coherence (MSC) are used to measure the linear relationship between two signals in the frequency domain.
Coherence is a complex-valued function that represents both the magnitude and phase of the relationship between two signals at each frequency.
\[ C_{xy} = \frac{P_{xy}}{\sqrt{P_{xx} P_{yy}}} \] \(P_{xy}\): cross power spectral density
\(P_{xx}\), \(P_{yy}\): auto-power spectral densities
The magnitude of coherence ranges from 0 to 1.
Magnitude: Strength of the linear relationship (0 = no relationship, 1 = perfect linear relationship).
Phase: Phase difference between the signals at each frequency.
Magnitude Squared Coherence (MSC) is the square of the magnitude of coherence. It is always a real number between 0 and 1.
\[ MSC_{xy} = \frac{\vert P_{xy} \vert ^2}{P_{xx} P_{yy}} = \vert C_{xy} \vert ^2 \]
MSC ranges from 0 to 1. Interpretation: 0: no linear relationship at that frequency, 1: perfect linear relationship at that frequency.
Use coherence if you need to know both the strength and the phase lag/lead between signals (e.g., for directional connectivity analysis).
Use MSC if you only care about the strength of the linear relationship (e.g., for detecting frequency-specific coupling).
Calculate coherence, imaginary part of coherence and magnitude-squared coherence using multi-taper method:
coh, imc, msc, f = NeuroAnalyzer.coherence(e10,
e10,
ch1="Fp1",
ch2="Fp2",
method=:mt);Get phases:
coh_ph = DSP.angle.(coh)1×2049×120 Array{Float64, 3}:
[:, :, 1] =
0.0 0.081198 -0.265702 -0.922933 … -2.92992 -2.93967 3.14159
[:, :, 2] =
0.0 -0.000313519 0.0210012 … 1.79804 0.108879 -0.0781269 0.0
[:, :, 3] =
0.0 -0.121932 -0.192776 -0.389349 … 2.41434 2.25677 1.95497 0.0
;;; …
[:, :, 118] =
0.0 0.118422 0.195129 0.223914 … 2.77749 2.81521 2.99892 3.14159
[:, :, 119] =
0.0 -0.413611 -0.864246 -1.1039 … 2.91877 2.85507 2.86219 3.14159
[:, :, 120] =
0.0 -0.000211489 -0.0221209 … 3.11961 3.13855 -3.10562 3.14159Get magnitudes:
coh = abs.(coh)1×2049×120 Array{Float64, 3}:
[:, :, 1] =
0.0788964 0.0634694 0.0392488 0.11019 … 0.558604 0.577596 0.585778
[:, :, 2] =
0.890396 0.892747 0.901109 0.906805 … 0.0422334 0.204124 0.294254
[:, :, 3] =
0.605836 0.582083 0.552101 0.561181 … 0.367192 0.298474 0.0126522
;;; …
[:, :, 118] =
0.944702 0.945704 0.942293 0.933107 … 0.288824 0.345377 0.396036
[:, :, 119] =
0.295663 0.322462 0.452344 0.654096 … 0.323764 0.280346 0.23787
[:, :, 120] =
0.79531 0.788859 0.766273 0.727637 … 0.328292 0.414555 0.460571Available methods: :mt (multi-taper), :fft (Fast Fourier transform), :stft (short time Fourier transform).
Plot coherence between O1 vs O2 in the alpha band:
eeg_alpha = NeuroAnalyzer.filter(e10,
ch="all",
fprototype=:butterworth,
ftype=:bp,
cutoff=band_frq(eeg, band=:alpha),
order=8)
coh, imc, msc, f = NeuroAnalyzer.coherence(eeg_alpha,
eeg_alpha,
ch1="O1",
ch2="O2",
ep1=1,
ep2=1,
frq_lim=band_frq(eeg, band=:alpha))
# get magnitude
coh = abs.(coh)
Plots.plot(f, coh[1, :, 1], label="coherence (magnitude)", xlabel="Frequency [Hz]")
Plots.plot!(f, msc[1, :, 1], label="MSC")
Plots.plot!(f, imc[1, :, 1], label="imaginary part")Also, plot_coherence() function is available:
coh, imc, mscoh, f = NeuroAnalyzer.coherence(eeg,
eeg,
ch1=get_channel(eeg, type="eeg")[1:4],
ch2=get_channel(eeg, type="eeg")[5:8],
ep1=1,
ep2=1,
frq_lim=band_frq(eeg, band=:alpha))
# get magnitude
coh = abs.(coh)
# simplify values
coh, _ = areduce(coh, f, n=0.1)
mscoh, f = areduce(mscoh, f, n=0.1)
# generate label pairs
l = paired_labels(get_channel(eeg, type="eeg")[1:4],
get_channel(eeg, type="eeg")[5:8])
# plot coherence values along frequencies
plot_coherence(coh[1, :, 1],
f,
title="Coherence $(l[1])")For multi-channel data:
plot_coherence(mscoh[:, :, 1],
f,
clabels=l,
title="Magnitude-squared coherence")plot_coherence_avg(mscoh[:, :, 1],
f,
clabels=l,
title="Averaged magnitude-squared coherence")plot_coherence_butterfly(mscoh[:, :, 1],
f,
clabels=l,
title="Magnitude-squared coherence")Plot averaged coherence values as connections:
c = zeros(19, 19)
coh_avg = round.(mean(coh[:, :, 1], dims=2), digits=2)
c[1, 5] = coh_avg[1]
c[2, 6] = coh_avg[2]
c[3, 7] = coh_avg[3]
c[4, 8] = coh_avg[4]
plot_locs(eeg,
ch="eeg",
connections=c)(!) Plotting connections requires locs to be added to the object.
(!) Line thickness represents connectivity value. If mono=true, positive values are represented with gray lines, negative – with dotted gray lines. If mono=false, positive values are represented with red lines, negative – with blue lines.
To presents values numerically, set the option weights to false.
plot_locs(eeg,
ch="eeg",
connections=c,
weights=false)(!) If mono=false, positive values are in red and negative values are in blue.