elephant.spectral.welch_coherence

elephant.spectral.welch_coherence(signal_i, signal_j, n_segments=8, len_segment=None, frequency_resolution=None, overlap=0.5, fs=1.0, window='hann', nfft=None, detrend='constant', scaling='density', axis=-1)[source]

Estimates coherence between a given pair of analog signals.

The estimation is performed with Welch’s method: the given pair of data are cut into short segments, cross-spectra are calculated for each pair of segments, and the cross-spectra are averaged and normalized by respective auto-spectra.

By default, the data are cut into 8 segments with 50% overlap between neighboring segments. These numbers can be changed through respective parameters.

Parameters:
signal_ineo.AnalogSignal or pq.Quantity or np.ndarray

First time series data of the pair between which coherence is computed.

signal_jneo.AnalogSignal or pq.Quantity or np.ndarray

Second time series data of the pair between which coherence is computed. The shapes and the sampling frequencies of signal_i and signal_j must be identical. When signal_i and signal_j are not neo.AnalogSignal, sampling frequency should be specified through the keyword argument fs. Otherwise, the default value is used (fs = 1.0).

n_segmentsint, optional

Number of segments. The length of segments is adjusted so that overlapping segments cover the entire stretch of the given data. This parameter is ignored if len_seg or frequency_resolution is given. Default: 8

len_segmentint, optional

Length of segments. This parameter is ignored if frequency_resolution is given. If None, it is determined from other parameters. Default: None

frequency_resolutionpq.Quantity or float, optional

Desired frequency resolution of the obtained coherence estimate in terms of the interval between adjacent frequency bins. When given as a float, it is taken as frequency in Hz. If None, it is determined from other parameters. Default: None

overlapfloat, optional

Overlap between segments represented as a float number between 0 (no overlap) and 1 (complete overlap). Default: 0.5 (half-overlapped)

fspq.Quantity or float, optional

Specifies the sampling frequency of the input time series. When the input time series are given as neo.AnalogSignal, the sampling frequency is taken from their attribute and this parameter is ignored. Default: 1.0

windowstr or tuple or np.ndarray, optional

Desired window to use. See Notes [1]. Default: ‘hann’

nfftint, optional

Length of the FFT used. See Notes [1]. Default: None

detrendstr or function or False, optional

Specifies how to detrend each segment. See Notes [1]. Default: ‘constant’

scaling{‘density’, ‘spectrum’}, optional

If ‘density’, computes the power spectral density where Pxx has units of V**2/Hz. If ‘spectrum’, computes the power spectrum where Pxx has units of V**2, if signal is measured in V and fs is measured in Hz. See Notes [1]. Default: ‘density’

axisint, optional

Axis along which the periodogram is computed. See Notes [1]. Default: last axis (-1)

Returns:
freqspq.Quantity or np.ndarray

Frequencies associated with the estimates of coherency and phase lag. freqs is always a vector irrespective of the shape of the input data. If signal_i and signal_j are neo.AnalogSignal or pq.Quantity, a pq.Quantity array is returned. Otherwise, a np.ndarray containing frequency in Hz is returned.

coherencynp.ndarray

Estimate of coherency between the input time series. For each frequency, coherency takes a value between 0 and 1, with 0 or 1 representing no or perfect coherence, respectively. When the input arrays signal_i and signal_j are multidimensional, coherency is of the same shape as the inputs, and the frequency is indexed depending on the type of the input. If the input is neo.AnalogSignal, the first axis indexes frequency. Otherwise, frequency is indexed by the last axis.

phase_lagpq.Quantity or np.ndarray

Estimate of phase lag in radian between the input time series. For each frequency, phase lag takes a value between \(-\pi\) and \(\pi\), with positive values meaning phase precession of signal_i ahead of signal_j, and vice versa. If signal_i and signal_j are neo.AnalogSignal or pq.Quantity, a pq.Quantity array is returned. Otherwise, a np.ndarray containing phase lag in radian is returned. The axis for frequency index is determined in the same way as for coherency.

Raises:
ValueError

Same as in welch_psd().

See also

welch_psd

Notes

  1. The computation steps used in this function are implemented in scipy.signal module, and this function is a wrapper which provides a proper set of parameters to scipy.signal.welch function.

  2. The parameters window, nfft, detrend, return_onesided, scaling, and axis are directly passed to the scipy.signal.welch function. See the respective descriptions in the docstring of scipy.signal.welch for usage.

  3. When only n_segments is given, parameter nperseg of scipy.signal.welch function is determined according to the expression

    signal.shape[axis] / (n_segments - overlap * (n_segments - 1))

    converted to integer.

Examples

>>> import neo
>>> import numpy as np
>>> import quantities as pq
>>> from elephant.spectral import welch_coherence
>>> signal = neo.AnalogSignal(np.cos(np.linspace(0, 2 * np.pi, num=100)),
...     sampling_rate=20 * pq.Hz, units='mV')

Sampling frequency will be taken as signal.sampling_rate.

>>> freq, coherency, phase_lag = welch_coherence(signal, signal)
>>> freq
array([ 0.        ,  0.90909091,  1.81818182,  2.72727273,  3.63636364,
        4.54545455,  5.45454545,  6.36363636,  7.27272727,  8.18181818,
        9.09090909, 10.        ]) * Hz
>>> coherency.flatten()
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
>>> phase_lag.flatten() 
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]) * rad