Feature WPLI

.gitignore

@@ -82,6 +82,8 @@ ignored/
 
 # data
 *.nix
+*.csv
+*.mat
 
 # neo logs
 **/logs

elephant/phase_analysis.py

@@ -16,7 +16,8 @@
     "spike_triggered_phase",
     "phase_locking_value",
     "mean_phase_vector",
-    "phase_difference"
+    "phase_difference",
+    "weighted_phase_lag_index"
 ]
 
 
@@ -307,3 +308,92 @@ def phase_difference(alpha, beta):
     delta = alpha - beta
     phase_diff = np.arctan2(np.sin(delta), np.cos(delta))
     return phase_diff
+
+
+def weighted_phase_lag_index(signal_i, signal_j, sampling_frequency=None,
+                             absolute_value=True):
+    r"""
+    Calculates the Weigthed Phase-Lag Index (WPLI)
+
+    This function estimates the WPLI, which is a measure of phase-synchrony. It
+    describes for two given signals i and j, who is leading/lagging the other
+    signal in the frequency domain across multiple trials.
+
+    Parameters
+    ----------
+    signal_i, signal_j : np.array, Quantity, neo.AnalogSignal
+        Time-series of the first and second signals,
+        with `t` time points and `n` trials.
+    sampling_frequency : quantity (default: None)
+        Sampling frequency of the signals in Hz. Not needed if signal i and j
+        are neo.AnalogSignals.
+    absolute_value : boolean (default: True)
+        Takes the absolute value of the numerator in the WPLI-formula.
+        When set to `False`, the WPLI contains additional directionality
+        information about which signal leads/lags the other signal:
+            - wpli > 0 : first signal leads second one
+            - wpli < 0 : first signal lags second one
+
+    Returns
+    -------
+    freqs : quantity
+        Positive frequencies in Hz associated with the estimates of `wpli`.
+        Range: :math:`[0, sampling frequency/2]`
+    wpli : float
+        Weighted phase-lag index of `signal_i` and `signal_j` across trials.
+        Range: :math:`[0, 1]`
+
+    Raises
+    ------
+    ValueError
+        If trial number or trial length are different for signal i and j.
+
+    Notes
+    -----
+    This implementation is based on the formula taken from [1] (pp.1550) :
+
+    .. math::
+        WPLI = \frac{| E( |Im(X)| * sgn(Im(X)) ) |}{E( |Im(X)| )}
+
+    with:
+        - E{...} : expected value operator
+        - Im{X} : imaginary component of the cross-spectrum
+        - X = Z_i * Z_j_conjugate : as cross-spectrum
+        - Z_i, Z_j : complex-valued matrix, representing the Fourier spectra
+                     of a particular frequency of the signals i and j.
+
+    References
+    ----------
+    ..  [1] Martin Vinck, Robert Oostenveld, Marijn van Wingerden,
+        Franscesco Battaglia, Cyriel M.A. Pennartz; "An improved index of
+        phase-synchronization for electrophysiological data in the presence
+        of volume-conduction, noise and sample-size bias"; NeuroImage, vol. 55,
+        pp. 1548-1565, 2011
+
+    """
+    if sampling_frequency is None:  # neo.AnalogSignal-input
+        sampling_frequency = signal_i.sampling_rate
+        signal_i = signal_i.magnitude
+        signal_j = signal_j.magnitude
+
+    if np.shape(signal_i) != np.shape(signal_j):
+        if len(signal_i) != len(signal_j):
+            raise ValueError("trial number of signal i and j must be equal")
+        raise ValueError("trial length of signal i and j must be equal")
+
+    # ARRAY-approach
+    # calculate Fourier transforms
+    fft1 = np.fft.rfft(signal_i)
+    fft2 = np.fft.rfft(signal_j)
+    freqs = np.fft.rfftfreq(np.shape(signal_i)[1], d=1.0 / sampling_frequency)
+
+    # obtain cross-spectrum
+    cs = fft1 * np.conjugate(fft2)
+    # calculate WPLI
+    wpli_num = np.mean(np.abs(np.imag(cs)) * np.sign(np.imag(cs)), axis=0)
+    if absolute_value:
+        wpli_num = np.abs(wpli_num)
+    wpli_den = np.mean(np.abs(np.imag(cs)), axis=0)
+    wpli = wpli_num / wpli_den
+
+    return freqs, wpli

elephant/test/cross_testing_scripts/generate_artificial_datasets_for_ground_truth_of_wpli.py

@@ -0,0 +1,92 @@
+import os
+
+import numpy as np
+from scipy.io import savemat
+
+
+def _generate_datasets_for_ground_truth_testing(ntrial=40, tlength=2.5,
+                                                srate=250):
+    """
+    Generates simple sinusoidal, artificial LFP-datasets.
+
+    This function simulates the recording of two LFP-signals over several
+    trials, for a certain duration of time and with a specified sampling
+    rate. These datasets will be used to calculate WPLI-ground-truth
+    with the MATlAB package FieldTrip and its function
+    ft_connectivity_wpli(), its wrapper ft_connectivityanalysis() and
+    with MNEs spectral_connectivity(). The last two use both multitaper
+    for FFT, therefore just certain frequencies will be compared.
+
+    Parameters
+    ----------
+    ntrial: int
+        Number of trials in the datasets.
+    tlength: float
+        Time length of the datasets in seconds.
+    srate: float
+        Sampling rate used to 'record' the signals in Hz.
+
+    Returns
+    -------
+    None
+
+    Notes
+    -----
+    Used versions of MATLAB & FieldTrip:
+    MATLAB          Version 9.9         (2020b)
+    FieldTrip       fieldtrip-20210128
+
+    Instead of using the ft_connectivityanalysis() wrapper function, which
+    expects preprocessed data, the ft_connectivity_wpli() is called
+    directly. That's because the preprocessing functions of FieldTrip have
+    no 'conventional' Fourier-Transformation (FFT) method, but among
+    others a multitaper-FFT. Because this python-version of WPLI doesn't
+    use multitaper, but conventional FFT, the utilized MATLAB script also
+    uses the conventional FFT to preprocess the data and calculate the
+    cross-spectrum, which will then be passed to the ft_connectivity_wpli.
+
+    """
+    times = np.arange(0, tlength, 1 / srate)
+
+    # lfps_1 & 2 will be used for
+    # 1) calculating ground-truth with FieldTrips' ft_connectivity_wpli
+    # 2) comparison to mutlitaper-approaches like
+    # FieldTrips' ft_connectivityanalysis()
+    # and MNEs' spectral_connectivity at certain frequencies
+    lfps_1 = [np.cos(2 * 16 * np.pi * times + np.pi / 2) +
+              np.cos(2 * 36 * np.pi * times + np.pi / 2) +
+              np.cos(2 * 52 * np.pi * times) +
+              np.cos(2 * 100 * np.pi * times) +
+              np.cos(2 * 70 * np.pi * times + np.pi / 2 + np.pi * (i % 2)) +
+              np.random.normal(loc=0.0, scale=1, size=len(times))
+              for i in range(ntrial)]
+    lfps_1 = np.stack(lfps_1, axis=0)
+
+    lfps_2 = [np.cos(2 * 16 * np.pi * times) +
+              np.cos(2 * 36 * np.pi * times) +
+              np.cos(2 * 52 * np.pi * times + np.pi / 2) +
+              np.cos(2 * 100 * np.pi * times + np.pi / 2) +
+              np.cos(2 * 70 * np.pi * times) +
+              np.random.normal(loc=0.0, scale=1, size=len(times))
+              for _ in range(ntrial)]
+    lfps_2 = np.stack(lfps_2, axis=0)
+
+    # save artificial LFP-dataset to .mat files
+    mdic_1 = {"lfp_matrix": lfps_1, "time": times, "sf": srate}
+    mdic_2 = {"lfp_matrix": lfps_2, "time": times, "sf": srate}
+    filename1_artificial = os.path.sep.join(['cross_testing_scripts',
+                                             'artificial_LFPs_1.mat'])
+    filename2_artificial = os.path.sep.join(['cross_testing_scripts',
+                                             'artificial_LFPs_2.mat'])
+
+    savemat(filename1_artificial, mdic_1)
+    savemat(filename2_artificial, mdic_2)
+
+
+def main():
+    """Run generation function."""
+    _generate_datasets_for_ground_truth_testing()
+
+
+if __name__ == '__main__':
+    main()

elephant/test/cross_testing_scripts/wpli_get_ground_truth_via_ft_connectivity_wpli.m

@@ -0,0 +1,36 @@
+%% use FieldTrips' ft_connectivity_wpli() to calculate WPLI ground-truth
+%% BEFORE using this script: change current folder to cross_testing_scripts
+
+% ARTIFICIAL DATA more complex
+dataset1 = load('artificial_LFPs_1.mat');
+dataset2 = load('artificial_LFPs_2.mat');
+
+lfp1 = dataset1.lfp_matrix;
+lfp2 = dataset2.lfp_matrix;
+Fs = double(dataset1.sf);
+
+siz = size(lfp1)
+trial = siz(1)
+tlength = siz(2)
+
+fft1 = realfft(lfp1, tlength, 2);
+fft2 = realfft(lfp2, tlength, 2);
+freq = 0:double(Fs/tlength):(Fs/2); %-0.1 to stay lower Fs/2
+
+cs = fft1 .* conj(fft2);
+length_fft = floor(tlength/2) +1
+cs = reshape(cs, [trial, 1, length_fft]);
+
+[wpli_1, v, n] = ft_connectivity_wpli(cs, 'feedback', 'text', 'dojack', 0, 'debias', 0);
+% wpli_1 = abs(wpli_1); %(optional)
+
+writematrix(wpli_1, 'ground_truth_WPLI_from_ft_connectivity_wpli_with_artificial_LFPs.csv');
+
+
+figure(1);
+plot(freq, wpli_1);
+
+function realfft = realfft(a, tlength, dim)
+     fft_all = fft(a, tlength, dim);
+     realfft = fft_all(:, 1:(floor(length(fft_all)/2)+1));
+end

elephant/test/cross_testing_scripts/wpli_get_ground_truth_via_ft_connectivityanalysis.m

@@ -0,0 +1,50 @@
+%% use FieldTrips' ft_connectivityanalysis() to calculate WPLI ground-truth
+%% BEFORE using this script: change current folder to cross_testing_scripts
+
+% ARTIFICIAL DATA
+dataset1 = load('artificial_LFPs_1.mat');
+dataset2 = load('artificial_LFPs_2.mat');
+
+siz = size(dataset1.lfp_matrix);
+ntrials = siz(1);
+nsamples = siz(2);
+
+% create a data-structure, which is compatible with FieldTrip
+joined_data.label = {'sig_1', 'sig_2'};
+joined_data.trial = cell(1, ntrials);
+for t = 1:ntrials
+    joined_data.trial(1, t) = {[dataset1.lfp_matrix(t, :); dataset2.lfp_matrix(t, :)]};
+end
+joined_data.time = cell(1, ntrials);
+for t = 1:ntrials
+    joined_data.time(1, t) = {dataset1.time};
+end
+joined_data.fsample = double(dataset1.sf);
+
+% do frequency-analysis
+cfg_freq            = [];
+cfg_freq.output     = 'fourier';
+% mtmfft
+cfg_freq.method     = 'mtmfft';
+cfg_freq.keeptrials = 'yes';
+cfg_freq.tapsmofrq  = 0.4;
+cfg_freq.taper      = 'dpss';
+cfg_freq.polyremoval= -1;
+
+freqfourier = ft_freqanalysis(cfg_freq, joined_data);
+
+% do WPLI-calculation
+cfg_wpli            = [];
+cfg_wpli.method     = 'wpli';
+ft_wpli      = ft_connectivityanalysis(cfg_wpli, freqfourier);
+
+% take absolute value (optional)
+wpli_2 = ft_wpli.wplispctrm(1, 2, 1:(end-1));
+wpli_2 = reshape(wpli_2, [1, length(wpli_2)]);
+% wpli_2 = abs(wpli_2);
+
+writematrix(wpli_2, 'ground_truth_WPLI_from_ft_connectivityanalysis_with_artificial_LFPs_multitaped.csv');
+
+
+figure(2);
+plot(ft_wpli.freq(1:(end-1)), wpli_2);

elephant/test/cross_testing_scripts/wpli_get_ground_truth_via_mne_spectral_connectivity.py

@@ -0,0 +1,56 @@
+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import quantities as pq
+import scipy
+from mne.connectivity import spectral_connectivity
+from scipy.io import savemat
+
+
+def main():
+    """
+    This script calculates a WPLI ground-truth with
+    MNE's spectral_connectivity(), which uses multitaper for FFT, so that
+    just certain frequencies will be compared to the results of
+    weighted_phase_lag_index(), which uses conventional FFT.
+    """
+    # Load first & second data file
+    filename1 = os.path.sep.join(['artificial_LFPs_1.mat'])
+    dataset1 = scipy.io.loadmat(filename1, squeeze_me=True)
+
+    filename2 = os.path.sep.join(['artificial_LFPs_2.mat'])
+    dataset2 = scipy.io.loadmat(filename2, squeeze_me=True)
+
+    # get the relevant values
+    lfps1 = dataset1['lfp_matrix'] * pq.uV
+    sf1 = dataset1['sf'] * pq.Hz
+    lfps2 = dataset2['lfp_matrix'] * pq.uV
+
+    # data-shape: epochs/trial X signals X times
+    lfp_both_signals = np.hstack((lfps1.magnitude, lfps2.magnitude))
+    lfp_both_signals = lfp_both_signals.reshape(np.shape(lfps1)[0], 2,
+                                                np.shape(lfps1)[1])
+
+    wpli, freqs, times, n_epochs, n_tapers = spectral_connectivity(
+        data=lfp_both_signals, method='wpli', sfreq=sf1.magnitude,
+        verbose=True, mode='fourier', fmin=0)
+
+    # get relevant relation: first signal compared to the second one
+    wpli = wpli[1][0][0:]
+    freqs = freqs[0:]
+
+    fig, ax = plt.subplots(1, 1, figsize=(10, 8), num=1)
+    fig.suptitle("Weighted Phase Lag Index - Ground Truth from MNE", size=20)
+    ax.plot(freqs, wpli, label="WPLI")
+    ax.set_xlabel('f (Hz)', size=16)
+    ax.legend(fontsize=16, framealpha=0)
+    plt.show()
+
+    np.savetxt("ground_truth_WPLI_from_MNE_spectral_connectivity"
+               "_with_artificial_LFPs_multitaped.csv",
+               wpli, delimiter=",", fmt="%20.20e")
+
+
+if __name__ == '__main__':
+    main()