<jupyter_start><jupyter_text>Analyse des signaux médicaux - Séance 1 de TP--At the end of this session, you will be able to : - Create and manage the Jupyter Notebooks environment to run code, insert text and math equations- Perform basic matrix manipulations using Numpy - Create signals and perform basic scientific computing using Scipy and Numpy- Produce simple data visulisation using Matplotlib - Open, manipulate and visualize an EEG dataset used in the course<jupyter_code>!pip install mne<jupyter_output><empty_output><jupyter_text>Part 1 - Intro to Jupyter Notebook--Here, we will only cover the basics. You can skip this part if you are already familiar with jupyter notebook. Jupyter Notebook is based on the .ipynb format (iPython Notebook), and is essentially a way to do rapid prototyping / demon of code / visualisation of scientific python. The basic idea is to define *cells*. Cells can be of several types, including python code, or rich text (using [markdown formatting](https://www.markdownguide.org/basic-syntax/)).When a code cell is evaluated (i.e. the python code will be executed), the output of this evaluation will show up right below the cell. When a text cell is evaluated, the text will be formatted. Please now do the "User Interface Tour" from the Help menu.  Done ? When working with Jupyter Notebook, you will essentially switch between two modes : - The Edit mode in which you edit the content of the cells - The Command mode, that enables you to change the cell types. When in Command mode, you can select cells. If you select a single cell, you can edit it by simply pressing enter, or double clicking on it.  For example, try editing THIS CELL and change its content.  Now, edit the cell below, change the code, and when you're done, press Shift+Enter to evalute the code.<jupyter_code>### CELL TO BE EDITED

a=32
b= 2*a
print("%d + %d"%(a,b))<jupyter_output><empty_output><jupyter_text>Text cells can contain math expressions that use the Markdown formatting, in which you can use LaTEx expressions for maths (enclosed between two dollar signs). For example : $A(k) \triangleq \sum_{\mathbf{n} =1}^{k}{n^2}$Now : - Edit the current cell to show the code that displays the math expression,- Create a code cell below that defines a function that calculates $A(k)$ given k, and evaluate this cell,- Create another cell and use the function to display $A(k)$ for a few values of $\mathbf{k}$ (eg 10 and 20).<jupyter_code>### CELL TO BE COMPLETED 

def A(k):
    result = 1
    ###  your code here
    return result
### CELL TO BE COMPLETED 
###  your code here<jupyter_output><empty_output><jupyter_text>Note that using Jupyter Notebook, if you evaluate a cell with a function followed by a "?" sign, the help of the function will pop up. Example :<jupyter_code>import os

os.listdir?<jupyter_output><empty_output><jupyter_text>You can also display the code of a function using the syntax "??"<jupyter_code>A??<jupyter_output><empty_output><jupyter_text>The popup can be closed by pressing the Escape key. Use the listdir function to browse the content of some directories...<jupyter_code>os.listdir('../')
os.listdir('./')<jupyter_output><empty_output><jupyter_text>Part 2 - Introduction to Numpy, Scipy and Matplotlib -- A code cell can contain any python code, including imports. Let's start by importing the Numpy package.<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>Numpy can be used to generate pseudo-random values from various distributions. In particular, a very useful distribution is the standard normal (zero mean and unit variance). Let's generate two vectors sampled from the normal distribution, using a length parameter that we'll be able to change if needed.<jupyter_code>length = 50

vecA = np.random.randn(length)
vecB = np.random.randn(2*length)<jupyter_output><empty_output><jupyter_text>vecA and vecB are numpy *arrays*. One of their attributes can be fetched to check their *shape*<jupyter_code>print(vecA.shape)
print(vecB.shape)
print(vecA)<jupyter_output><empty_output><jupyter_text>Numpy arrays can be vectors as well as matrices, or any tensor. For example the following code will create a 3 x length matrix, using the standard normal<jupyter_code>matC = np.random.randn(3,length)
matC.shape<jupyter_output><empty_output><jupyter_text>This could be used for instance for EEG signals in the time frequency domain. Let's imagine we have Power Spectral Density in 10 frequency bands, 1 second at 500 Hz, and 5 electrodes. We can generate a random signal with such shape by using numpy random generators.Check the syntax of the np.random.randn function, and generate an array called EEG_timefreq, with the specifications  (10 frequency bands, 1 second at 500 Hz, and 5 electrodes)<jupyter_code>### EEG_timefreq =  ## TO BE COMPLETED (2-1)<jupyter_output><empty_output><jupyter_text>Note that the random package of Numpy has several other interesting functions. Try to test the two functions proposed in the cell below. Try uncommenting the two functions below one by one, look up their help page, and try to use them.<jupyter_code>### CELL TO BE COMPLETED  (2-2)

#np.random.randint
#np.random.permutation<jupyter_output><empty_output><jupyter_text>A very important features of arrays is the fact they can be used as *iterables*. For example, you can iterate over the dimensions of an array by simply "looping" over it using a *for* loop<jupyter_code>for freqband in EEG_timefreq:
    print(freqband.shape)<jupyter_output><empty_output><jupyter_text>This way, we can loop over frequency bands, and process each frequency band separately.  Also possible to enumerate along the dimension in order to get the index of the current "smaller" array<jupyter_code>print('Initial shape is %d %d %d' % (EEG_timefreq.shape[0],EEG_timefreq.shape[1],EEG_timefreq.shape[2]))
print('Iterating over the first dimension using an index k')
for k,curdim in enumerate(EEG_timefreq):
    print('k = %d, shape is %d %d' % (k,curdim.shape[0],curdim.shape[1]))<jupyter_output><empty_output><jupyter_text>Now let's imagine we want to calculate the average power over time and electrode, for each couple frequency band.  Use the previous principle in order to calculate the average of each frequency band subvector, using the function np.mean()<jupyter_code>### CELL TO BE COMPLETED (2-3)<jupyter_output><empty_output><jupyter_text>Check that you obtain the same result when directly computing the average over the two axis 1 and 2 (look up the arguments of np.mean), which correspond to averaging in time and electrode.<jupyter_code>### CELL TO BE COMPLETED (2-4)<jupyter_output><empty_output><jupyter_text>These features will prove to be very useful when manipulate large arrays.  Another important operation when working with Numpy Arrays is *reshaping*. Essentially, *reshaping* consists in changing the organisation of the array (in terms of dimension), while keeping the same number of elements. For example, a 20x10 2D array can be converted into a 4x5x10 array<jupyter_code>A = np.random.randint(1,5,(10,20))
print('Initial shape of A is %d x %d' % (A.shape[0],A.shape[1]))
print(A)
B = A.reshape((4,5,10))
print('B is A reshaped to %d x %d x %d' % (B.shape[0],B.shape[1],B.shape[2]))
print(B)<jupyter_output><empty_output><jupyter_text>Finally, a property of numpy that is really important is broadcasting. The goal of broadcasting is to simplify the vectorization of certain operations when the vectors do not have the same shape. For example you can easily perform element-wise multiplication.To test this try doing an element-wise multiplication of the vector x and matrix y below<jupyter_code>### To be completed (2-5)


print("X: ",x)
print("Y: ")
print(y)
print("X shape is: ",x.shape)
print("Y shape is: ",y.shape)
print("Element-wise multiplication shape:", result.shape)
print("Element-wise multiplication result:")
print(result)<jupyter_output><empty_output><jupyter_text>Another very powerful tool in numpy is indexing. You can use either an integer vector or a boolean vector to choose which indexes you want to extract from your numpy tensor.Consider that we want to extract all elements from the first line of your vector y that have a higher value than 1, you would have to do:<jupyter_code>first_row = y[0]
first_row_higher_than_one = first_row > 1
print("Result: ", first_row[first_row_higher_than_one])<jupyter_output><empty_output><jupyter_text>You can also choose specific lines to query, for example if you want to query lines 0 and 2<jupyter_code>rows = [0,2]
rows_result = y[rows]
values_higher_than_one = rows_result > 1
print("Result: ", rows_result[values_higher_than_one])<jupyter_output><empty_output><jupyter_text>You can also save and load your numpy tensors using np.savez and np.load. This will be really important as this enable you to generate your data only once instead of having to do all the calculations every time you need your data.<jupyter_code>filename = "x.npz"
source_tensor = x
np.savez(filename,data=source_tensor)
loaded_npz = np.load(filename)
loaded_tensor = loaded_npz["data"]
print("Your tensor was loaded and contains: ", loaded_tensor)<jupyter_output><empty_output><jupyter_text>Part 3 - Application to an EEG + ECG dataset-- Let's load some EEG + ECG data acquired with OpenBCI and visualize it<jupyter_code>import matplotlib.pyplot as plt
### The next command is not python, but is specific to jupyter.

%matplotlib inline

from mne.io import read_raw_edf # This is the function to open the raw data

import matplotlib.pyplot as plt # We also import the plotting functions from matplotlib<jupyter_output><empty_output><jupyter_text>The following cell opens the data and fixes the data structure - just run it.  Before running this cell, if you are using Google Colab, you should upload the file "eeg_ecg.edf" on your drive, from the github repository : https://github.com/brain-bzh/health_exg/blob/master/TP/data/eeg_ecg.edf<jupyter_code>filepath = 'data/eeg_ecg.edf' #### Set path to file here, 

Raw = read_raw_edf(filepath,stim_channel=False) 
# The second argument is neeeded, otherwise the function assumes that the last channel 
# is a "stimulation" channel, but in our case we only have EEG channels

#from mne.channels import read_montage
# We read the electrode position from the standard template positions
#montage = read_montage(kind='standard_1020',ch_names=Raw.ch_names)

scales = dict(eeg='auto') ### This is necessary because the data acquired by openBCI 

Raw.load_data()<jupyter_output><empty_output><jupyter_text>Information about the data can be obtained by examining the "info" structure from the Raw object. Let's print it .<jupyter_code>print(Raw.info)<jupyter_output><empty_output><jupyter_text>Most importantly, we can see the channel names and the sampling frequency. The data can be fetched into a numpy array, using the "get_data()" method.<jupyter_code>exg_array = Raw.get_data()<jupyter_output><empty_output><jupyter_text>What is the shape of this array ? What do the respective dimensions correspond to ? What is the duration of the recording (in seconds) ?<jupyter_code>### To be completed (3-1)<jupyter_output><empty_output><jupyter_text>ECG Analysis -- The Raw.ch_names gives the names of all electrodes. Which index corresponds to the ECG ?<jupyter_code>### TO BE COMPLETED  (3-2)<jupyter_output><empty_output><jupyter_text>The Raw.times vector gives the time points of measurements, in seconds. Thus, it can be used to plot the data. In the cell below: - create an time_indices vector, containing the indices of the first five seconds of the Raw.times vector- generate a x vector as using those time indices- create a y vector that selects the first five seconds of ECG- use plt.plot(x,y) to plot the first five seconds of data of ECG<jupyter_code>#### TO BE COMPLETED (3-3)<jupyter_output><empty_output><jupyter_text>Can you identify the PQRST complex ? How long is approximately a whole cycle ? What about the T-wave ?  Now we would be interested in analysis Heart Rate Variability, by extracting the times corresponding to the peak of the R wave. Try to perform this using the [find_peaks](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.find_peaks.html) function from scipy and find the peaks corresponding to the first 5 seconds of ECGA few things to note : - find_peaks works for local maxima. - You have to set a minimal distance between two peaks<jupyter_code>from scipy.signal import find_peaks

### TO BE COMPLETED (3-4)<jupyter_output><empty_output><jupyter_text>Once you have found the times of the peak, plot them aligned with the ECG signal. You can use plt.plot, followed by plt.stem, like this (assuming you have named R_times the vector containing the times correspond to the peaks of the R waves)<jupyter_code>plt.plot(x,y)
plt.stem(R_times,200*np.ones_like(R_times),'r')<jupyter_output><empty_output><jupyter_text>This seems to work fine on the first five seconds. Now, can you do it on the whole recording ? A longer section of the recording ?<jupyter_code>### TO BE COMPLETED (3-5)<jupyter_output><empty_output><jupyter_text>Next step : we are going to estimate heart rate variability. First, compute the durations between peak R times.<jupyter_code>#### TO BE COMPLETED (3-6)<jupyter_output><empty_output><jupyter_text>Visualize all those durations using an histogram, provided by matplotlib plt.hist()<jupyter_code># To be completed (3-7)<jupyter_output><empty_output><jupyter_text>Calculate the average heart rate in Beats Per Minute<jupyter_code>### To be computed (3-8)<jupyter_output><empty_output><jupyter_text>Part 4 - EEG Analysis -- Now we will analyse the first four electrodes of the exg_array, corresponding to EEG electrodes. What are the names of those electrodes ? Look up how they are placed by checking the 10/20 EEG electrode placement.<jupyter_code>### To be completed (4-1)<jupyter_output><empty_output><jupyter_text>Use what we have learned so far to plot a few segments of EEG data. What can we notice when looking at the signal globally ?<jupyter_code>### To be completed (4-2)<jupyter_output><empty_output><jupyter_text>We will preprocess this data using filtering. This can be done using the method "filter" from MNE.More background on filtering can be found [here](https://mne.tools/stable/auto_tutorials/discussions/plot_background_filtering.htmldisc-filtering)The filter method modifies the data in place, so we work on a copy.<jupyter_code>raw_filt = Raw.copy()<jupyter_output><empty_output><jupyter_text>We will remove two things with filtering : slow trend, and high frequency noise due to line noise (> 40 Hz)<jupyter_code>#TO be completed (4-3)<jupyter_output><empty_output><jupyter_text>Plot the result<jupyter_code># To be completed (4-4)<jupyter_output><empty_output><jupyter_text>Estimating the Power Spectral Density of Raw Data-- We are using [MNE](https://mne.tools/stable/auto_tutorials/index.html) to load and process the EEG data. The data is managed using an object, Raw, which includes methods for visualization, filtering, and much more. We can directly plot the Power Spectral density (PSD) using the plot_psd method of the Raw object. A few things : - Use the "picks" argument to select only the first four electrodes, as the fifth one is the ECG.- Compare the filtered and raw version<jupyter_code>### TO BE COMPLETED  (4-5)<jupyter_output><empty_output><jupyter_text>Finding EOG artifacts-- MNE comes with very nice methods to automatically find artefacts. This is explained in depth [here](https://mne.tools/stable/auto_tutorials/preprocessing/plot_10_preprocessing_overview.htmlsphx-glr-auto-tutorials-preprocessing-plot-10-preprocessing-overview-py). Quickly go through this tutorial, and apply the create_eog_epochs to find all EOG artefacts in the provided data.<jupyter_code>from mne.preprocessing import create_eog_epochs
create_eog_epochs?
# TO BE COMPLETED (4-6)<jupyter_output><empty_output><jupyter_text>This technique can be combined with ICA, to automatically identify an ICA component that correspond to eye blinks. This is explained in detail [here](https://mne.tools/stable/auto_tutorials/preprocessing/plot_40_artifact_correction_ica.htmlusing-an-eog-channel-to-select-ica-components). Can you perform ICA on this dataset, and remove the EOG artefacts?<jupyter_code># TO BE COMPLETED (4-7)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07942898975109809, 'fraction_numerical': 0.008113714006832601, 'mean_word_length': 5.669243287225386, 'pattern_counts': {'":': 0, '<': 169, '<?xml version=': 0, '>': 172, 'https://': 7, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '42', 'dir': 'jupyter-structured-clean-dedup', 'id': '264924', 'license': 'MIT', 'n_tokens_mistral': 5058, 'n_tokens_neox': 4829, 'n_words': 2008, 'path': 'TP/SignauxMedicaux_ECG_EEG.ipynb', 'repo_name': 'brain-bzh/health_exg', 'size': '16392'}