Skip to content

Commit

Permalink
update of exercises
Browse files Browse the repository at this point in the history
  • Loading branch information
behinger committed Oct 23, 2024
1 parent e198858 commit 2f09d0d
Show file tree
Hide file tree
Showing 5 changed files with 391 additions and 513 deletions.
32 changes: 18 additions & 14 deletions exercises/ex3_filter.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -19,20 +19,22 @@ x = cos(2*pi*0.5*t) + 0.2*sin(2*pi*2.5*t+0.1) + \
```


## Task
Plot the signal against time.
**Task:** Plot the signal against time.



## Transform to fourier space

Next run the FFT algorithm.
:::: callout-tip
## Task

plot the magnitude response ($log10(abs(fft(x))))$).
- run the fourier analysis via `fft`.
- plot the full log10-magnitude
::::

::: callout-note
If you plot e.g. the magnitude you will notice it is mirrorsymmetric around the middle! This is because FFT is actually defined for complex signals. Because we only have real signals, practically half of the FFT is redundant. But this also means, all of our manual-filters in the next steps need to take the mirror-symmetry into account. For this reason, typically a `rfft` function exists, returning only the real part.
:::

## A simple filter

A very simple frequency filter is to set the unwanted amplitude of the fourier components to something close to zero. We have to be a tad careful not to remove the phase information though.
Expand All @@ -49,9 +51,7 @@ A very simple frequency filter is to set the unwanted amplitude of the fourier c
- plot the signal with what you started out
::::

::: callout-note
If you plot e.g. the magnitude you will notice it is mirrorsymmetric around the middle! This is because FFT is actually defined for complex signals. Because we only have real signals, practically half of the FFT is redundant. But this also means, all of our manual-filters need to take the mirror-symmetry into account.
:::


## Highpass instead of lowpass
:::: callout-tip
Expand All @@ -62,21 +62,24 @@ Repeat the steps from above, but this time, remove the low frequency components


## What happens to the frequency and time response if we add "artefacts"?
Let's see how the FFT reacts to steps and spikes. A step-function could look like this:

![](https://upload.wikimedia.org/wikipedia/commons/thumb/d/d9/Dirac_distribution_CDF.svg/325px-Dirac_distribution_CDF.svg.png)

:::: callout-tip
## Task

Add a DC-offset (a step-function) starting from `x[200:]` and investigate the fourier space. Filter it again (low or high pass) and transfer it back to the time domain and investigate the signal around the spike.
Add a DC-offset (a step-function) starting from `x[200:]` and investigate the fourier space. Filter it again (low or high pass) and transfer it back to the time domain and investigate the signal around the step.
::::

## Impulse Response Function
To get a bit deeper understanding of what is going on, have a look at the fourier transform of a new impulse signal (e.g. `1:400 => 0`. and `200 => 1`.).
To get a bit deeper understanding of what is going on, have a look at the fourier transform of a impulse signal (e.g. `1:400 => 0`. and `200 => 1`.).

:::: callout-tip
## Question

What do you observe?

Why would we see ringing if we put most of the coefficients to 0?
- What do you observe?
- Why did we see ringing in the previous filtered step-function, if when filtering we put coefficients to 0 - after all, we are not adding something?
::::

## Filtering EEG data
Expand All @@ -103,7 +106,7 @@ raw.load_data()
:::: callout-tip
## Task

- Choose the channel "Pz", plot the channel (previous HW it was "Cz")
- Choose the channel "Pz", plot the channel
- Plot the fourier space using `raw.plot_psd()`
::::

Expand All @@ -118,6 +121,7 @@ Now we filter using `raw.filter()`, specify a highpass of 0.5Hz and a lowpass of

::: callout-note
## Bonus-task
**Bonus: comparing over-filtered ERPs**

If you want, you can compare the ERP with and without filtering. You can also use "invalid" filter settings - HP up to 2-5Hz, lowpass until 10-20Hz. I say invalid here, because usually with such ranges, you would filter out results that you are actually interested in.

Expand Down
35 changes: 28 additions & 7 deletions exercises/ex4_cleaning.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -2,24 +2,45 @@
## Cleaning Data
We will re-use the dataset we used in the first exercise.

**T:** Go through the dataset using the MNE explorer and clean it. You can use `raw.plot()` for this. If you are working from a jupyter notebook, try to use `%matplotlib qt` for better support of the cleaning window. To get an understanding how the tool works, press `help` or type `?` in the window. (Hint: You first have to add a new annotation by pressing `a`)
**Task:** Go through the dataset using the MNE explorer and clean it. You can use `raw.plot()` for this.

**T:** While going through the dataset, mark what you observe as bad electrodes. Those are saved in `raw.info['bads']`. The channels can be interpolated with `raw.interpolate_bads()` or `epoch.interpolate_bads()`. Compare the channel + neighbours before and after. Did the interpolation succeed? (If you are interested in the mathematical details of spline interpolation, checkout this https://mne.tools/dev/overview/implementation.html#id26)
Hint: You need channel locations to run the interpolation which you can get by using the default-standardized channel locations `raw.set_montage('standard_1020',match_case=False)`
:::callout-tip
If you are working from a jupyter notebook, try to use `%matplotlib qt` for better support of the cleaning window. To get an understanding how the tool works, press `help` or type `?` in the window. (Hint: You first have to add a new annotation by pressing `a`)
:::

**T:** Save the annotations to a file. Best practice is to use a csv file of sorts e.g. using `raw.annotations.save(filename)`. Try to subsetting for `BAD_` first. Bonus: Save it in a BIDS derivate folder according to the BIDS guidelines.
## Marking electrodes
**Task:** While going through the dataset, mark what you observe as bad electrodes. Those are saved in `raw.info['bads']`.

**T:** In MNE, if annotations with `BAD_` exist, the epoching function automatically removes them. We are now ready to compare ERP results with and without removal of bad segments. Epoch the data, the respective bad segments will be removed automatically. Compare the two ERPs for the channel `Cz`

**T:** In the epoching step, we can also specify rejection criterion for a peak-to-peak rejection method
The channels can later be interpolated with `raw.interpolate_bads()` or `epoch.interpolate_bads()`. Compare the channel + neighbours before and after. Did the interpolation succeed?

You need channel locations to run the interpolation which you can get by using the default-standardized channel locations `raw.set_montage('standard_1020',match_case=False)`

:::callout-tip
If you are interested in the mathematical details of spline interpolation, checkout [this link](https://mne.tools/dev/overview/implementation.html#id26)
:::

**Task:** Save the annotations to a file. Best practice is to use a csv file of sorts e.g. using `raw.annotations.save(filename)`. Try to subsetting for `BAD_` first. **Bonus:** Save it in a BIDS derivate folder according to the BIDS guidelines.

:::callout-tip
In the semester project you will use mne-bids-pipeline, which will automatically save many results - you might nevertheless come to a point where saving your own logs might be a good idea
:::

## Remove bad trials
In MNE, if annotations with `BAD_` exist, the epoching function automatically removes them. We are now ready to compare ERP results with and without removal of bad segments. Epoch the data, the respective bad segments will be removed automatically.

**Task:** Compare the two ERPs for the channel `Cz`

## Automatic rejection
In the epoching step, we can also specify rejection criterion for a peak-to-peak rejection method

```python
reject_criteria = dict(eeg=100e-6, # 100 µV
eog=200e-6) # 200 µV
epochs = mne.Epochs(raw, events, reject=reject_criteria,reject_by_annotation=False)
```

Compare these epochs with your manual rejection and with the ERPs without rejection. Plot a single channel `Cz` overlaying all three "solutions".
**Task:** Compare these epochs with your manual rejection and with the ERPs without rejection. Plot a single channel `Cz` overlaying all three "solutions".

**Bonus:** We are going to add to our comparison by using **autoreject** in MNE python. Have a look at https://autoreject.github.io/ for installation and usage examples. Hint: You need channel locations to run autoreject which you can get by using the default-standardized channel locations `raw.set_montage('standard_1020',match_case=False)`

Expand Down
41 changes: 27 additions & 14 deletions exercises/ex5_ICA.qmd
Original file line number Diff line number Diff line change
@@ -1,19 +1,24 @@
# Signal processing and analysis of human brain potentials (EEG) [Exercise 5]
## Overview
**T:** Load simulation dataset by `x = ccs_eeg_utils.simulate_ICA(dims=2)` and plot it. We will write simple ICA algorithm to try to return the original signals
We will write simple ICA algorithm to try to return the original unmixed signals based on a simulation example.

**Task:** load the simulation dataset by `x = ccs_eeg_utils.simulate_ICA(dims=2)` and plot it.

Remember that for ICA we have to undo a scaling and a rotation. The scaling we can do by whitening our signal `x`. The sphering (=whitening) matrix is the inverse of the matrix-squareroot of the covariance matrix of x.

In order to whiten (`x_white`), we have to remove the mean of `x` and then calculate the matrix product with the sphering matrix `sph*x`.
Now we are ready to try a extermely simple, naive, restricted and stupid ICA algorithm:

We simply rotate the signal using a 2D rotation matrix and for each rotatio calculate the kurtosis
**Task:** Whiten the data!

Now we are ready to try a extermely simple, naive, restricted and inefficient ICA algorithm:

We will rotate the signal using a 2D rotation matrix and for each rotatio calculate the sum of the absolute valued kurtosis (use e.g. `scp.stats.kurtosis`). Later we will try to find the maximum value, which will be the rotation for the most non-normal signal - hopefully the unmixed one

```python
turn =lambda k: np.reshape(np.array([np.cos(k), np.sin(k), -np.sin(k), np.cos(k)]),(2,2))
```

(pseudocode)
**Task:** Implement the following pseudocode

```
for t in 0:pi
Expand All @@ -23,17 +28,25 @@ for t in 0:pi

The solution to ICA would be then `argmax(kurtosis_list)`.

**T:** Plot the solution!
**Task:** Plot the solution!

**Bonus** If you want to be fancy, you can also generate an animation linking kurtosis + rotating the datapoints + marginal histograms together

**T:** Next we will use the infomax algorithm implemented in `mne` (you could try implementing it yourself, but I found an effective adjustment of learning rate quite tricky).
you can call the infomax algorithm using `mne.preprocessing.infomax(x.T,verbose=True)`. Run it on `x = ccs_eeg_utils.simulate_ICA(dims=4)` and plot the signals before and after.
## Infomax
Next we will use the infomax algorithm implemented in `mne` (you could try implementing it yourself, but I found an effective adjustment of learning rate quite tricky).
you can call the infomax algorithm using `mne.preprocessing.infomax(x.T,verbose=True)`.

**Task:** Run it on `x = ccs_eeg_utils.simulate_ICA(dims=4)` and plot the signals before and after.


## ICA on EEG data
**T:** Now we are ready to run ICA on real EEG data. Select a dataset from a previous exercise and start the ICA (you can downsample the data to speed up ICA calculation
`raw.resample(100)`).
Now we are ready to run ICA on real EEG data.

**Task:** Select the dataset from a previous exercise and start the ICA

:::callout-hint
You could downsample the data to speed up ICA calculation via `raw.resample(100)`.
:::
```python
from mne_bids import (BIDSPath,read_raw_bids)
import mne_bids
Expand All @@ -53,17 +66,17 @@ raw.set_montage('standard_1020',match_case=False)

```

Generate an ICA object
1. Generate an ICA object
`ica= mne.preprocessing.ICA(method="fastica")`

fit the ICA on the data.
2. fit the ICA on the data.
`ica.fit(raw,verbose=True)`

After the fit (can take some minutes) you can plot the component using `ica.plot_components()`
3. After the fit (can take some minutes) you can plot the component using `ica.plot_components()`

If you want to inspect a specific component, you can select it using `ica.plot_properties(raw,picks=[3])`.
4. If you want to inspect a specific component, you can select it using `ica.plot_properties(raw,picks=[3])`.

In order to exclude components, write: ica.exclude = [3,4,16]
5. In order to exclude components, write: `ica.exclude = [3,4,16]`

You can then compare the EEG before & after using:
```python
Expand Down
366 changes: 103 additions & 263 deletions exercises/solutions/ex1_overview.ipynb

Large diffs are not rendered by default.

Loading

0 comments on commit 2f09d0d

Please sign in to comment.