bci_essentials.resting_state

  1"""
  2A module for processing resting state data.
  3
  4The EEG data inputs for each function are either a single trials or
  5a set of trials.
  6- For single trials, inputs are of the shape `n_channels x n_samples`, where:
  7    - n_channels = number of channels
  8    - n_samples = number of samples
  9- For multiple trials, inputs are of the shape `n_trials x n_channels x n_samples`, where:
 10    - n_trials = number of trials
 11    - n_channels = number of channels
 12    - n_samples = number of samples
 13
 14"""
 15
 16import numpy as np
 17import scipy.signal
 18
 19from matplotlib import pyplot as plt
 20
 21from .utils.logger import Logger  # Logger wrapper
 22
 23# Instantiate a logger for the module at the default level of logging.INFO
 24# Logs to bci_essentials.__module__) where __module__ is the name of the module
 25logger = Logger(name=__name__)
 26
 27
 28def get_shape(data):
 29    """Get the shape of the input data.
 30
 31    Parameters
 32    ----------
 33    data : numpy.ndarray
 34        Trial(s) of resting state EEG data.
 35        2D or 3D array containing data with `float` type.
 36
 37        shape = (`n_channels`,`n_samples`) OR
 38        (`n_trials`,`n_channels`,`n_samples`)
 39
 40    Returns
 41    -------
 42    n_trials : int
 43        Number of trials.
 44    n_channels : int
 45        Number of channels.
 46    n_samples : int
 47        Number of samples.
 48
 49    """
 50    try:
 51        n_trials, n_channels, n_samples = np.shape(data)
 52    except Exception:
 53        n_channels, n_samples = np.shape(data)
 54        n_trials = 1
 55
 56    return n_trials, n_channels, n_samples
 57
 58
 59# This function is never used at the moment, anywhere in the code. Renaming to more clear convention on it being a public function.
 60def get_bandpower(data, fs, fmin, fmax, normalization=None):
 61    """Get the bandpower of a trial of EEG.
 62
 63    Parameters
 64    ----------
 65    data : numpy.ndarray
 66        A single resting state EEG trial
 67        2D array containing data with `float` type.
 68
 69        shape = (`n_channels`,`n_samples`)
 70    fs : float
 71        EEG sampling frequency.
 72    fmin : float
 73        Lower frequency bound.
 74    fmax : float
 75        Upper frequency bound.
 76    normalization : string, *optional*
 77        Method for normalization.
 78        - `"norm"` for divide by norm.
 79        - `"sum"` for divide by sum.
 80        - Default is `None` and no normalization is done.
 81
 82    Returns
 83    -------
 84    power : numpy.ndarray
 85        Bandpower of frequency band.
 86
 87        shape = (`n_channels`)
 88
 89    """
 90    n_channels, n_samples = data.shape
 91
 92    f, Pxx = scipy.signal.welch(data, fs=fs)
 93
 94    # Normalization
 95    if normalization == "norm":
 96        Pxx = np.divide(Pxx, np.tile(np.linalg.norm(Pxx, axis=1), (len(f), 1)).T)
 97    if normalization == "sum":
 98        Pxx = Pxx / Pxx.sum(axis=1).reshape((Pxx.shape[0], 1))
 99
100    ind_local_min = np.argmax(f > fmin) - 1
101    ind_local_max = np.argmax(f > fmax) - 1
102
103    power = np.zeros([n_channels])
104    for channel in range(n_channels):
105        power[channel] = np.trapz(
106            Pxx[channel, ind_local_min:ind_local_max], f[ind_local_min:ind_local_max]
107        )
108
109    return power
110
111
112# Alpha Peak
113def get_alpha_peak(data, alpha_min=8, alpha_max=12, plot_psd=False):
114    """Get the alpha peak based on the all channel median PSD.
115
116    Parameters
117    ----------
118    data : numpy.ndarray
119        Resting state EEG trial with eyes closed.
120        3D array containing data with `float` type.
121
122        shape = (`n_trials`,`n_channels`,`n_samples`)
123    alpha_min : float, *optional*
124        Lowest possible value of alpha peak (Hz)
125        - Default is `8`.
126    alpha_max : float, *optional*
127        Highest possible value of alpha peak (Hz)
128        - Default is `12`.
129    plot_psd : bool, *optional*
130        Plot the PSD for inspection.
131        - Default is `False`.
132
133    Returns
134    -------
135    alpha_peaks : numpy.ndarray
136        The peak alpha frequency (in Hz) for each trial.
137    """
138
139    fs = 256
140
141    n_trials, n_channels, n_samples = get_shape(data)
142
143    # Create alpha_peaks of length n_trials
144    alpha_peaks = np.zeros(n_trials)
145
146    for trial in range(n_trials):
147        # Get the current trial
148        current_trial = data[trial, :, :]
149
150        # Calculate PSD using Welch's method, nfft = n_samples
151        f, Pxx = scipy.signal.welch(current_trial, fs=fs, nperseg=n_samples)
152
153        # Limit f, Pxx to the band of interest
154        ind_min = scipy.argmax(f > alpha_min) - 1
155        ind_max = scipy.argmax(f > alpha_max) - 1
156
157        f = f[ind_min:ind_max]
158        Pxx = Pxx[:, ind_min:ind_max]
159
160        alpha_peaks[trial] = f[np.argmax(np.median(Pxx, axis=0))]
161        logger.info("Alpha peak of trial %s is %s", trial, alpha_peaks[trial])
162
163        if plot_psd:
164            nrows = int(np.ceil(np.sqrt(n_channels)))
165            ncols = int(np.ceil(np.sqrt(n_channels)))
166
167            fig, axs = plt.subplots(nrows, ncols, figsize=(10, 8))
168            for r in range(nrows):
169                for c in range(ncols):
170                    ch = (ncols * r) + c
171                    axs[r, c].set_title(ch)
172                    axs[r, c].plot(f, Pxx[ch, :])
173
174            plt.show()
175
176            plt.figure()
177            plt.plot(f, np.median(Pxx, axis=0))
178
179            plt.show()
180
181    return alpha_peaks
182
183
184# Bandpower features
185def get_bandpower_features(data, fs, transition_freqs=[0, 4, 8, 12, 30]):
186    """Get bandpower features.
187
188    Parameters
189    ----------
190    data : numpy.ndarray
191        Trials of resting state EEG data.
192        3D array containing data with `float` type.
193
194        shape = (`n_trials`,`n_channels`,`n_samples`)
195    fs : float
196        Sampling frequency (Hz).
197    transition_freqs : array-like, *optional*
198        The transition frequencies of the desired power bands.
199        The first value is the lower cutoff, the last value is the
200        upper cutoff, the middle values are band transition frequencies.
201        - Default is `[0, 4, 8, 12, 30]`.
202
203    Returns
204    -------
205    abs_bandpower : numpy.ndarray
206        A numpy array of the absolute bandpower of provided bands.
207        The last item is the total bandpower. The array length is equal
208        to the length of `transition_freqs`.
209    rel_bandpower : numpy.ndarray
210        A numpy array of the relative bandpower of provided bands with
211        respect to the entire region of interest from `transition_freqs[0]`
212        to `transition_freqs[-1]`. The last item is the total relative
213        bandpower and should always equal `1`. The array length is equal to
214        the length of `transition_freqs`.
215    rel_bandpower_mat : array_like
216        A numpy array of the relative bandpower of provided bands such that
217        location `(R,C)` corresponds to the power of `R` relative to `C`.
218        The final row and column correspond to the cumulative power of
219        all bands such that `(R,-1)` is `R` relative to all bands
220
221    """
222    # Get Shape
223    n_trials, n_channels, n_samples = get_shape(data)
224
225    # Initialize
226    abs_bandpower = np.zeros((len(transition_freqs), n_trials))
227    rel_bandpower = np.zeros((len(transition_freqs), n_trials))
228    rel_bandpower_mat = np.zeros(
229        (len(transition_freqs), len(transition_freqs), n_trials)
230    )
231
232    # for each trial
233    for trial in range(n_trials):
234        # Get the current trial
235        current_trial = data[trial, :, :]
236
237        # Calculate PSD using Welch's method
238        f, Pxx = scipy.signal.welch(current_trial, fs=fs)
239
240        # Limit f, Pxx to the band of interest
241        ind_global_min = scipy.argmax(f > min(transition_freqs)) - 1
242        ind_global_max = scipy.argmax(f > max(transition_freqs)) - 1
243
244        f = f[ind_global_min:ind_global_max]
245        Pxx = Pxx[:, ind_global_min:ind_global_max]
246
247        # Calculate the absolute power of each band
248        for tf in range(len(transition_freqs)):
249            # The last item is the total
250            if tf == len(transition_freqs) - 1:
251                abs_bandpower[tf, trial] = np.sum(abs_bandpower[:tf, trial])
252                continue
253
254            fmin = transition_freqs[tf]
255            fmax = transition_freqs[tf + 1]
256
257            ind_local_min = np.argmax(f > fmin) - 1
258            ind_local_max = np.argmax(f > fmax) - 1
259
260            # Get power for each channel
261            abs_power = np.zeros([n_channels])
262            for ch in range(n_channels):
263                abs_power[ch] = np.trapz(
264                    Pxx[ch, ind_local_min:ind_local_max], f[ind_local_min:ind_local_max]
265                )
266
267            # Median across all channels
268            abs_bandpower[tf, trial] = np.median(abs_power)
269
270        rel_bandpower[:, trial] = abs_bandpower[:, trial] / abs_bandpower[-1, trial]
271
272        # Calculate the relative power of each band
273        for tf1 in range(len(transition_freqs)):
274            for tf2 in range(len(transition_freqs)):
275                rel_bandpower_mat[tf1, tf2, trial] = (
276                    abs_bandpower[tf1, trial] / abs_bandpower[tf2, trial]
277                )
278
279    return abs_bandpower, rel_bandpower, rel_bandpower_mat