= Tutorial 10: Frequency filters =
''Authors: Francois Tadel, Elizabeth Bock, John C Mosher, Sylvain Baillet''

<<TableOfContents(2,2)>>

== Evaluation of the noise level ==
We are now going to process our continuous recordings to remove the main sources of noise. Typically, we can expect contaminations coming from the environment (power lines, stimulation equipment, building vibrations) and from the subject (movements, blinks, heartbeats, breathing, teeth clenching, muscle tension, metal in the mouth or the body). In this tutorial will focus first on the contaminations that occur continuously, at specific frequencies.

Before running any type of cleaning procedure on MEG/EEG recordings, we recommend to always start with a quick evaluation of the noise level. An easy way to do this is to estimate the power spectrum of all the sensors over the entire recordings.

 * Clear the list of files in the Process1 tab.
 * Select all the three datasets we have linked to our protocol. <<BR>>You can select the three "link to raw file" nodes, the three folders or the entire subject node. <<BR>><<BR>> {{attachment:psd_select_files.gif||height="121",width="400"}}
 * Click on [Run] to open the pipeline editor window.
 * Select the process "'''Frequency > Power spectrum density (Welch)'''"
 * This process evaluates the power of the MEG/EEG signals at different frequencies, using the Welch's method (see [[http://en.wikipedia.org/wiki/Welch's_method|Wikipedia]] or [[http://www.mathworks.com/help/signal/ref/pwelch.html|MathWorks]]). It splits the signals in overlapping windows of a given length, calculates the Fourier Transform (FFT) of each of those short segments, and average the power of the FFT coefficients for all the overlapping windows. <<BR>><<BR>> {{attachment:psd_options.gif||height="388",width="607"}}
 * Set the options as follows (click on [Edit] for the additional options):
  * '''Time window''': [All file] <<BR>>Portion of the input file you want to use the for estimating the spectrum.
  * '''Window length''': 4 seconds <<BR>>Estimator length = length of the overlapping time windows for which we calculate the FFT. The number of time samples in the estimator is the same as the number of frequency bins in the output file. Increasing this parameter increases the output frequency resolution (distance between two frequency bins) but degrades the stability of the estimator, because we decrease the total number of averaged time windows.
  * '''Overlap''': 50% <<BR>>How much overlap do we want between two consecutive two time windows.
  * '''Sensor types or names''': MEG <<BR>>In some processes, you can specify the names or types of sensors on which you want to apply the process. This way you can for instance apply different filters on the EEG and the MEG, if you have both in the same files.
  * '''Frequency definition''': Matlab's FFT default <<BR>>You have the option to use directly the frequency binning returned by the FFT, or do an additional step of averaging those bins in larger frequency bands. Note that you can edit those frequency bands.
  * '''Output''': Save individual PSD value. <<BR>>This option will estimate separately the PSD for each of the three files in input, and create three files in output. If you select the other option (save average), it calculates the same three files but averages them on the fly and saves only one file in the database.
 * Click on [Run] to start the execution of the process.
 * '''Troubleshooting''': If you get "Out of memory" errors, try to run this PSD estimation on a shorter time segment. For example, set the time window to [0,100s] instead of the full file. This process starts by loading all the needed recordings in memory, you might not have enough memory available on your system to fit the entire dataset.
 * It produces three new files, that appear as "depending" on the three datasets we have linked. The comment of the files indicate how many overlapping windows could be used to estimate the PSD. "179/4000ms" means 179 windows of 4s each (total 716s). With the 50% overlap, it sums up to a little less than 2x the file length (360s). <<BR>><<BR>> {{attachment:psd_output_files.gif||height="240",width="235"}}

== Interpretation of the PSD ==
==== File: AEF#01 ====
 * Double-click on the new PSD file for run #01 to display it (or right-click > Power spectrum). <<BR>><<BR>> {{attachment:psd_display.gif}}
 * The power spectrum is displayed in a figure similar to the time series, but the X axis represents the frequencies. Most of the shortcuts described for the recordings are also valid here. Clicking on the white parts of the figure or using the arrow keys of the keyboard move the frequency cursor. The current frequency can also be controlled with a new slider displayed in the Brainstorm window, just below the time panel.
 * Each black line represent the spectrum for one channel. If you click on a channel, it is selected and highlighted in red. Click again for deselecting. Right-click on a channel to read its name. <<BR>><<BR>> {{attachment:psd_run1.gif||height="194",width="515"}}
 * The shape of this graph is normal, it does not indicate anything unexpected:
  * Peaks related with the subject's alpha rhythms: around '''10Hz''' and '''20Hz'''
  * Peaks related with the power lines: '''60Hz, 120Hz''' and '''180Hz'''
 * Add a topography view for the same file, with one of the two methods below:
  * Right-click on the PSD file > '''2D Sensor cap'''.
  * Right-click on the spectrum figure > '''2D Sensor cap''' (shortcut: '''Ctrl+T''')
 * Scroll in frequencies to see the spatial distribution of each frequency bin: <<BR>><<BR>> {{attachment:psd_topo.gif||height="149",width="604"}}
 * We have already identified two artifacts we will need to remove: the eye movements and the contamination from the power lines.

==== File: AEF#02 ====
 * Open the spectrum view for the run AEF#02. <<BR>><<BR>> {{attachment:psd_run2.gif||height="154",width="389"}}
 * Add a 2D sensor cap view for the same file. Scroll to zoom in/out.<<BR>>To display the sensor names, right-click on the figure > Channels > Display sensors.
 * This spectrum looks very similar to the run #01: same alpha and power line peaks.
 * Additionally, we observe higher level of noise in frequencies in the range of '''30Hz to 100Hz''' on '''many occipital sensors'''.   This is probably due to some tension in the '''neck muscles''' due to an   uncomfortable position. We will see later whether those channels need to   be tagged as bad. <<BR>><<BR>> {{attachment:psd_neck.gif||height="198",width="565"}}

==== File: Noise recordings ====
 * Open the spectrum view for the noise recordings. <<BR>><<BR>> {{attachment:psd_noise.gif||height="139",width="346"}}
 * This shows the power spectrum of the signals that are recorded when there is no subject in the MEG room. It gives a good and simple  representation of the instrumental noise. If you had one bad MEG sensor, you would see it immediately in this graph. Here everything looks good.

== Power line contamination ==
[Start from here]

Notch  filters are adapted to remove some well  identified  contaminations from  oscillating systems, such as the power  lines 50Hz or  60Hz. Here is an  example of how to evaluate and correct  fixed frequency  artifacts.

One of the typical  pre-processing steps consist in getting rid of the contamination due to  the power lines (50 Hz or 60Hz).

Usually we prefer to apply the frequency filters '''before '''the any other type of correction.

 1. Put '''ALL''' the "Link to raw file" into the Process1 box (or directly the Subject01 folder)

 1. Run the process: '''Pre-process > Notch filter''' {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=sinRemoval.gif|sinRemoval.gif|height="272",width="333",class="attachment"}}
  * Select the frequencies: '''60, 120, 180 Hz'''
  * Sensor types or names: '''MEG'''
  * The higher harmonics are too high to bother us in this analysis, plus they are not clearly visible in all the recordings.
  * In  output, this process creates new .ds folders in the same folder as the  original files, and links the new files to the database.
  {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=psd3.gif|psd3.gif|class="attachment"}}
 1. Run again the PSD process "'''Frequency > Power spectrum density (Welch)'''" on those new files, with the same parameters, to evaluate the quality of the correction.
 1. Double-click on the new PSD files to open them.
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=psd5.gif|psd5.gif|height="177",width="416",class="attachment"}}
 1. Zoom in with the mouse wheel to observe what is happening around 60Hz (before / after).
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=psd6.gif|psd6.gif|height="143",width="453",class="attachment"}}
 1. To  avoid the confusion later, delete the links to the original files:  Select the folders containing the original unfiltered files and press  the Delete key (or right-click > File > Delete).
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=psd7.gif|psd7.gif|height="184",width="356",class="attachment"}}

The   operation creates a new "Link to raw file" entry in the database,   pointing to a new CTF dataset in the same folder as the original   continuous file. To check where this file was created: right-click on   the file "'''Raw | notch(60Hz 120Hz 180Hz)'''" > '''File '''> '''Show in file explorer'''.   This file has the exact same properties as the original file,  including  the SSP projectors, only the values of the MEG sensors were  updated.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=sinRemovalDb.gif|sinRemovalDb.gif|height="144",width="326",class="attachment"}}

== Alternatives ==
If the notch filter is not giving satisfying result, you can use two other processes:

 * '''Sinusoid removal''':  This process can do a better job at removing precise frequencies by    identifying the sinusoidal components and then subtracting them from    the signals in time domain (not a frequency filter).
 * '''Band-stop filter''':  This process can be useful for removing larger segments of the  spectrum, in case the power line peaks are spread over numerous  frequency bins or for suppressing other types of artifacts.

== Evaluate the results ==
Now run the same frequency analysis: process "'''Frequency > Power spectrum density (Welch)'''", time window = '''[0,50]s''', window length = '''4000ms''', overlap = '''50%'''.  Double click on the new "Power (MEG)" file. As expected, the three  first peaks due to the power lines are gone.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=psdAfter.gif|psdAfter.gif|height="177",width="411",class="attachment"}}

You   can run a few other pre-processing operations directly on the   continuous files, including a band-pass filter. Frequency filters and   sinusoid removals are operations that you should apply at very early   stages of the analysis, before epoching the recordings. Those operations  do not perform well next to the beginning and the end of the signals,  they may generate important artifacts. It is  therefore much more  accurate to filter the entire recordings from the  original continuous  file at once, rather than filtering small epochs  after importing them  in the database.

From now on, we are only going to work on this clean file "Raw | notch(60Hz 120Hz 180Hz)".

<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/PipelineEditor&next=Tutorials/ArtifactsDetect")>>

<<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/ArtifactsFilter)>>