Eeg dataset
===========

The EegDataset class of nnsa provides a very useful tool for handling EEG data.
It relies on the TimeSeries class and serves is an extension of the BaseDataset class.
This script provides a short introduction on how to use the EegDataset class.

Author: Tim Hermans (tim-hermans@hotmail.com).

Link to script: `containers/eeg_dataset.py <https://gitlab.com/timhermans/nnsa_public/-/blob/main/examples/containers/eeg_dataset.py>`_



.. code-block:: python

    import os
    
    from nnsa import EdfReader, EegDataset, NotchIIR, Butterworth, moving_std, RemezFIR
    import matplotlib.pyplot as plt
    import numpy as np
    
    plt.close('all')


Specify the path to an EDF file to process. 



.. code-block:: python

    
    # Filepath to .EDF file with EEG.
    fp_edf = r'C:/data_temp/test.edf'


Read EEG from the EDF file (if the file exists, otherwise create some dummy data). 



.. code-block:: python

    if os.path.exists(fp_edf):
        # Load the EEG from the EDF.
        with EdfReader(fp_edf) as r:
            # Returns a nnsa.EegDataset
            ds = r.read_eeg_dataset()
    else:
        print(f'File {fp_edf} not found. Creating dummy data... '
              f'Note that the artefact detection model will recognize that its fake and predict all artefacts.')
        fs = 250
        channel_labels = ['EEG Fp1', 'EEG Fp2', 'EEG C3', 'EEG C4', 'EEG Cz', 'EEG T3', 'EEG T4', 'EEG O1', 'EEG O2']
        eeg = (np.random.rand(len(channel_labels), fs*40) - 0.5)*300
        ds = EegDataset.from_array(eeg=eeg, fs=fs, channel_labels=channel_labels, unit='uV')
    
    # Print info string for EegDataset.
    print(ds)


We can do several types of prepocessing, such as referencing, filtering and resampling. For most of such functions, you can typically choose whether to do this inplace or not. 



.. code-block:: python

    # Let's look at referencing first.
    
    # Reference to Cz. If not inplace, then a new object is returned. In contrast, if inplace is True, then the data is
    # processed inplace and nothing is returned.
    ds_ref = ds.reference('Cz', inplace=False)
    print(ds_ref)
    
    # Create bipolar channels.
    ds_bp = ds.create_bipolar_channels(
        channels_1=['Fp2', 'Fp2', 'C4', 'Fp1', 'Fp1', 'C3'],
        channels_2=['C4',  'O2',  'O2', 'C3',  'O1',  'O1']
    )
    print(ds_bp)


Now let's look at filtering. To filter the data, you first need to define a filter based on the nnsa.FilterBase() class. Common filters are already implemented, see the nnsa.preprocessing.filter. For example, a Notch and a Butterworth filter (see also their documentations for details on their usage). 



.. code-block:: python

    notch_filt = NotchIIR(f0=50)
    but_filt = Butterworth(fn=[0.27, 30], order=1)


Now its very easy to apply these filters to the eeg data. You can choose whether you want to do zero-phase filtering (`filtfilt`) or just one-way filtering (`filter`). 



.. code-block:: python

    
    # Filter with notch using zero-phase filtering.
    ds_notch = ds_ref.filtfilt(notch_filt)
    
    # Filter with Butterworth using normal (one-way) filtering.
    ds_notch_but = ds_notch.filter(but_filt)


Resampling is also easy. You can choose two methods for resampling: (anti-aliasing) "polyphase_filtering" or "interpolation". If you there is a risk of aliasing, the polyphase_filtering might be the safest option. 



.. code-block:: python

    
    # Resample using interpolation (the high frequencies have already been filtered out and interpolation is faster).
    ds_notch_but_res = ds_notch_but.resample(fs_new=128, method='interpolation')


By default, these functions do not operate inplace, which makes it easy to chain these functions to make the code shorter: 



.. code-block:: python

    ds_preproc = ds.reference('Cz').filtfilt(NotchIIR(f0=50)).filter(Butterworth(fn=[0.27, 30], order=1)).resample(
        fs_new=128, method='interpolation')


The transform method may also come in handy at times, if you want to transform each channel according to a certain function. 



.. code-block:: python

    
    # For example, compute the moving standard deviation in a moving 2 second window.
    ds_movmean = ds.transform(fun=lambda x: moving_std(x, n=int(ds.fs*2)))


Of course, you can also obtain the EEG data in a numpy array to do other processing of the EEG data (and optionally put is back into an EegDataset format, if convenient). 



.. code-block:: python

    
    # Get the EEG data as a numpy array.
    eeg_array = ds.asarray()  # or np.asarray(ds)
    
    # The corresponding sampling frequency, time array, and channel labels can be obtained, too:
    fs = ds.fs
    time = ds.time
    channel_labels = ds.channel_labels


Plotting EEG data is easy using the plot() method (see the method's documentation for more details and options): 



.. code-block:: python

    fig, ax = plt.subplots(tight_layout=True)
    ds_bp.plot(ax=ax)


The EegDataset class has a bunch of methods to easily compute some features or do some analysis. 



.. code-block:: python

    
    # For example, detect artefacts (see also artefact detection examples).
    af_ds = ds.detect_artefacts()  # Returns an EegDataset with True and False values.
    print(af_ds)
    
    # Burst detection (see also burst detection examples).
    ds_preprocessed = ds.reference('Cz').resample(fs_new=250).filtfilt(RemezFIR(passband=[1, 20], stopband=[0.5, 21]))
    bd_result = ds_preprocessed.burst_detection(method='line_length')  # Returns a BurstDetectionResult object.
    print(bd_result)


For a list of all functions, see the documentation of EegDataset: https://nnsa.readthedocs.io/en/latest/nnsa.containers.html#nnsa.containers.datasets.EegDataset 




If you have multi-channel data which is not all EEG, you might check out the BaseDataset class, which is the parent class of the EegDataset and therefore shares some functionality, see https://nnsa.readthedocs.io/en/latest/nnsa.containers.html#nnsa.containers.datasets.BaseDataset