Source code for tvb.analyzers.node_coherence

# -*- coding: utf-8 -*-
#
#
# TheVirtualBrain-Scientific Package. This package holds all simulators, and
# analysers necessary to run brain-simulations. You can use it stand alone or
# in conjunction with TheVirtualBrain-Framework Package. See content of the
# documentation-folder for more details. See also http://www.thevirtualbrain.org
#
# (c) 2012-2023, Baycrest Centre for Geriatric Care ("Baycrest") and others
#
# This program is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software Foundation,
# either version 3 of the License, or (at your option) any later version.
# This program is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# You should have received a copy of the GNU General Public License along with this
# program.  If not, see <http://www.gnu.org/licenses/>.
#
#
#   CITATION:
# When using The Virtual Brain for scientific publications, please cite it as explained here:
# https://www.thevirtualbrain.org/tvb/zwei/neuroscience-publications
#
#

"""
Compute cross coherence between all nodes in a time series.

.. moduleauthor:: Stuart A. Knock <Stuart@tvb.invalid>
.. moduleauthor:: Marmaduke Woodman <marmaduke.woodman@univ-amu.fr>

"""

import numpy
import matplotlib.mlab as mlab
from matplotlib.pylab import detrend_linear
import tvb.datatypes.spectral as spectral
from tvb.basic.logger.builder import get_logger
from tvb.basic.neotraits.api import narray_describe

log = get_logger(__name__)


# TODO: Should do this properly, ie not with mlab, returning both coherence and
#      the complex coherence spectra, then supporting magnitude squared
#      coherence, etc in a similar fashion to the FourierSpectrum datatype...


def _hamming(M, sym=True):
    """
    The M-point Hamming window.
    From scipy.signal
    """
    if M < 1:
        return numpy.array([])
    if M == 1:
        return numpy.ones(1, 'd')
    odd = M % 2
    if not sym and not odd:
        M = M + 1
    n = numpy.arange(0, M)
    w = 0.54 - 0.46 * numpy.cos(2.0 * numpy.pi * n / (M - 1))
    if not sym and not odd:
        w = w[:-1]
    return w


[docs] def coherence_mlab(data, sample_rate, nfft=256): _, nsvar, nnode, nmode = data.shape # (frequency, nodes, nodes, state-variables, modes) coh_shape = nfft/2 + 1, nnode, nnode, nsvar, nmode coh = numpy.zeros(coh_shape) for mode in range(nmode): for var in range(nsvar): data = data[:, var, :, mode].copy() data -= data.mean(axis=0)[numpy.newaxis, :] for n1 in range(nnode): for n2 in range(nnode): cxy, freq = mlab.cohere(data[:, n1], data[:, n2], NFFT=nfft, Fs=sample_rate, detrend=detrend_linear, window=mlab.window_none) coh[:, n1, n2, var, mode] = cxy return coh, freq
def _coherence(data, sample_rate, nfft=256, imag=False): "Vectorized coherence calculation by windowed FFT" nt, ns, nn, nm = data.shape nwin = nt // nfft if nwin < 1: raise ValueError( "Not enough time points ({0}) to compute an FFT, given a " "window size of nfft={1}.".format(nt, nfft)) # ignore leftover data; need shape (nn, ... , nwin, nfft) wins = data[:int(nwin * nfft)]\ .copy()\ .transpose((2, 1, 3, 0))\ .reshape((nn, ns, nm, nwin, nfft)) wins *= _hamming(nfft) F = numpy.fft.fft(wins) fs = numpy.fft.fftfreq(nfft, 1e3 / sample_rate) # broadcasts to [node_i, node_j, ..., window, time] G = F[:, numpy.newaxis] * F.conj() if imag: G = G.imag dG = numpy.array([G[i, i] for i in range(nn)]) C = (numpy.abs(G)**2 / (dG[:, numpy.newaxis] * dG)).mean(axis=-2) mask = fs > 0.0 # C_ = numpy.abs(C.mean(axis=0).mean(axis=0)) return numpy.transpose(C[..., mask], (4, 0, 1, 2, 3)), fs[mask]
[docs] def calculate_cross_coherence(time_series, nfft): """ # type: (TimeSeries, int) -> CoherenceSpectrum # Adapter for cross-coherence algorithm(s) # Evaluate coherence on time series. Parameters __________ time_series : TimeSeries The TimeSeries to which the Cross Coherence is to be applied. nfft : int Data-points per block (should be a power of 2). """ srate = time_series.sample_rate coh, freq = _coherence(time_series.data, srate, nfft=nfft) log.debug("coherence") log.debug(narray_describe(coh)) log.debug("freq") log.debug(narray_describe(freq)) spec = spectral.CoherenceSpectrum( source=time_series, nfft=nfft, array_data=coh.astype(numpy.float64), frequency=freq) return spec