Differences between revisions 6 and 52 (spanning 46 versions)

Tutorial 12: Signal Space Projections (SSP)

Authors: Francois Tadel, Elizabeth Bock, John C Mosher, Sylvain Baillet

As previously said, the frequency filters are not adapted to remove artifacts that are transient or overlapping in frequency domain with the brain signals of interest. Other approaches exist to correct for those artifacts, based on the spatial signature of the artifacts.

If an event is very reproducible and occurs always at the same position (eg. eye blinks and heartbeats), the sensors will always record the same values when it occurs. We can identify the topographies corresponding to this artifact (ie. the spatial distributions of values at one time point) and remove them from the recordings. This spatial decomposition is the basic idea behind two widely used approaches: the SSP (Signal-Space Projection) and ICA (Independent Component Analysis) methods.

This introduction tutorial will focus on the SSP approach, as it is a lot simpler and faster but still very efficient for removing blinks and heartbeats from MEG and high-density EEG recordings. The interface for running ICA decompositions is very similar and will be described in an advanced tutorial.

Contents

Overview
The order matters
SSP: Heartbeats
Evaluate the components
Evaluate the correction
SSP: Eye blinks
Run #02
SSP: Generic
Averaged artifact
Troubleshooting
SSP Theory
References
SSP Algorithm
On the hard drive
Additional documentation

Overview

The general SSP objective is to identify the sensor toporagraphies that are typical of a specific artifact, then to create spatial projectors to remove the contributions of these topographies from the recordings.

We start by identifying many examples of the artifact we are trying to remove. This is what we've been doing in the previous tutorial with the creation of the "cardiac" and "blink" events.
We extract a short time window around each of these event markers and concatenate in time all those small blocks of recordings.
We run a principle components analysis (PCA) on the concatenated artifacts in order to get a decomposition in various spatial components (number of components = number of sensors).
If it works well, we can find in the first few principal components some topographies that are very specific of the type of artifact we are targeting. We select those components to remove.
We compute a linear projector for each spatial component to remove and save them in the database (in the "Link to raw file"). They are not applied immediately to the recordings.
Whenever some recordings are read from this file, the SSP projectors are applied on the fly to remove the artifact contributions. This approach is fast and memory efficient.

The order matters

This procedure has to be repeated separately for each artifact type. The order in which you process the artifacts matters, because for removing the second artifact we typically use the recordings cleaned with the first set of SSP projectors. We have to decide which one to process first.

It works best if each artifact is defined precisely and as independently as possible from the other artifacts. If the two artifacts happen simultaneously, the SSP projectors calculated for the blink may contain some of the heartbeat topography and vice versa. When trying to remove the second artifact, we might not be able to isolate it clearly anymore.

Because the heart beats every second or so, there is a high chance that when the subject blinks there is a heartbeat not too far away in the recordings. Therefore a significant number of the blinks will be contaminated with heartbeats. But we have usually a lot of "clean" heartbeats, we can start by removing those ones. To isolate correctly these two common artifacts, we recommend the following procedure:

Remove the markers "cardiac" that are occurring during a blink (done in the previous tutorial),
Compute the cardiac SSP (with no eye movements, because we removed the co-occurring events),
Compute the blink SSP (with no heartbeats, because they've already been taken care of).

If you have multiple modalities recorded simultaneously, for example MEG and EEG, you should run twice this entire procedure, once for the EEG only and once for the MEG only. You will always get better results if you process the different types of sensors separately. Same thing when processing Elekta-Neuromag recordings: process separately the magnetometers (MEG MAG) and the gradiometers (MEG GRAD).

SSP: Heartbeats

Double-click on the link to show the MEG sensors for Run #01.
In the Record tab, select the menu: "Artifacts > SSP: Heartbeats".

Event name: Name of the event to use to calculate the projectors, enter "cardiac".
Sensor types: Type of sensors for which the projection should be calculated ("MEG"). Note that you will always get better results if you process the different types of sensors separately.
Compute using existing SSP projectors: You have the option to calculate the projectors from the raw recordings, or from the recordings filtered with the previously computed SSP projectors. Unless you have a good reason for not considering the existing projectors, you should select this option.

After the computation is done, a new figure is displayed, that lets you select the active projectors.

On the left: The projector categories where each row represents the result of an execution of this process (usually one for each sensor type and each artifact).
On the right: The spatial components returned by the PCA decomposition. The percentage indicated between brackets is the singular value for this each component, normalized for this decomposition (percentage = Si / sum(Si), see technical details at the end of this page).
Percentage: More practically, it indicates the amount of signal that was captured by the component during the decomposition. The highest it is, the more the component is representative of the artifact recordings that were used to calculate it. In the good cases, you would typically see one to three components with values that are significantly higher than the others.
When a component is selected, it means that it is removed from the recordings. A spatial projector is computed and applied to the recordings on the fly when reading from the continuous file.
Default selection: The software selects the components with values superior to 12% and leaves the others unselected. This threshold is completely empirical and depends on the acquisition device, you should always review manually the components that you want to remove.

Evaluate the components

For the cardiac event, none of the components show a value superior to 12%, therefore the entire "cardiac" projector category was unselected. This doesn't mean that none of those components is actually interesting for us. The percentage indicated for the first value (9%) is much higher than the following ones (5%, 5%, 4%, 3%...), this could indicate that it targets relatively well the cardiac artifact. Let's investigate this in more details.

Select the category "cardiac" to enable the list on the right.
Click on the first component, then click on the toolbar button [Display component topography]. This menu shows the spatial distribution of the sensor values for this component. Note that you don't have to select the component (ie. check the box) to display it. This topography seems to correspond to a strong dipolar activity located relatively far from the sensor array, it matches the type of artifact we expect from the heart activity.
The second button "Display component topography [No magnetic interpolation]" produces the same figure but without the reinterpolation of the magnetic fields that is typically applied to the MEG recordings in Brainstorm, it may help understand some difficult cases. This magnetic interpolation will be detailed later in the introduction tutorials.
You can display multiple components in the same figure: select them at the same time in the list (holding the Shift/Ctrl/Cmd button of your keyboard) and then click on the button "Display topography".
No other strong components looks like it could be related with the heartbeats.
The last button in the toolbar [Display component time series], opens a figure that represents the evolution of the contribution of this component over time. The higher the amplitude, the more present the selected topography in the recordings. Click on it to show the component #1, then display the ECG signal at the same time (right-click on the file > ECG > Display time series).
We observe that the "SSP1" trace correlates relatively well with the ECG trace, in the sense that we captured most the ECG peaks with this component. However, the component seems also to capture much more signal than just the heartbeats: many alpha oscillations and some of the ocular activity. The example below shows a blink in the EOG, ECG and SSP component #1.
If you remove this component from the recordings, you can expect to see most of the artifacts related with the cardiac activity to go away, but you will also remove additional signal elements that were not really well identified. The job is done but it causes some unwanted side effects.
It is in general possible to refine the SSP decomposition by going back to the selection of "cardiac" markers that we used to compute it. You could look at all the ECG peaks individually and remove the markers located in segments of recordings that are noisier or that contain a lot of alpha activity (~10Hz). You would need to delete this SSP decomposition and run again the same process.
Alternatively, or if you don't manage to extract a clean cardiac component with a PCA/SSP decomposition, you could try to run an ICA decomposition instead. Note that for some subjects, the cardiac artifact is not very strong and could be simply ignored in the analysis.

Evaluate the correction

The topography of the component #1 looks like it represents a cardiac topography and its temporal evolution shows peaks where we identified heartbeats. It is therefore a good candidate for removal, we just need to make sure the signals look good after the correction before validating this choice.

Show the left-temporal MEG sensors (CTF LT) and select/unselect the first SSP component.
Repeat this for different time windows, to make sure that the cardiac peaks in the MEG sensors really disappear when the projector #1 is selected and that the rest is not altered too much.
No correction:
Cardiac component #1 removed:
In this example we will consider that the current decomposition is good enough.
Make sure you select the component #1, then click on [Save] to validate the modifications.
After this window is closed, you can always open it again from the Record tab with the menu Artifacts > Select active projectors. At this stage of the analysis, you can modify the list of projectors applied to the recordings at any time.

SSP: Eye blinks

Let's try the same thing with the eye blinks.

Select the process "SSP > Compute SSP: Eye blinks"
Run it on the event type "blink", that indicates the peaks of the EOG signal.
You see now a new category of projectors. It is selected by default because its first component shows a percentage higher than 12% (empirical and somewhat random threshold). Based on the distribution of values, this first component is most likely a good representation of the artifact we are trying to remove. The second one could be a good candidate as well.
Select the first three components and display their topographies:
- Component #1: Most likely a blink,
- Component #2: Probably a saccade (other type of eye movement we were not targeting),
- Component #3: Not related with the eye movements (maybe related with the alpha activity).
Display the times series for those three components, together with the EOG signals. You have to uncheck temporarily the component #1 to be able to display its signal. When it is checked, it is removed from the signal therefore it is corresponds to a flat trace.
The figure below shows the EOG and SSP values between 318s and 324s. The SSP1 trace matches the blink observed in VEOG and SSP2 matches the saccade observed in HEOG.
Observe the activity when there is no component selected:
With only the component #2 selected (saccade):
With components #1 and #2 selected (blink + saccade):
Keep the components #1 and #2 selected and click on [Save] to validate your changes.

Run #02

Reproduce the same operations on Run #02:

Close everything with the [X] button at the top-right corner of the Brainstorm window.
Open the MEG recordings for run AEF #02 (double-click on the file link).
Artifacts > SSP: Heartbeats: Event name "cardiac", sensors "MEG.
Select component #1, click on [Save].
Artifacts > SSP: Heartbeats: Event name "cardiac", sensors "MEG.
Select component #1, click on [Save].
Note that in this second session, the representation of the saccade was not as clear as in the first file. The distribution of the percentage values does not show any clear component other from the blink one, and the topographies are not as clear. In general, the saccade processing requires a separate step, we will illustrate this in the next tutorial.

Advanced

SSP: Generic

The calculation of the SSP for the heartbeats and the eye blinks are shortcuts to a more generic algorithm "Compute SSP: Generic". You may need this process if the standard parameters do not work of if you want to use this technique to remove other types of artifacts.

Time window: What segment of the file do you want to consider.
Event name: Markers that are used to characterize the artifact. If you don't specify any event name, it will use this entire time window as the input of the PCA decomposition.
Event window: Time segment to consider before and after each event marker. We want this time window to be longer than the artifact effect itself, in order to have a large number of time samples representative of a normal brain activity. This helps the PCA decomposition to separate the artifact from the ongoing brain activity.
Frequency band: Definition of the band-pass filter that is applied to the recordings before calculating the projector. Usually you would use the same frequency band as we used for the detection, but you may want to try to refine this parameter if the results are not satisfying.
Sensor types or names: List of sensors for which to calculate the projector. You can get better results if you process one sensor type at a time.
Compute using existing SSP/ICA projectors: Same as in the heartbeats/blinks processes.
Save averaged artifact in the database: If you check this option with an event name selected, the process will save the average of all the artifact epochs (event marker + event window) before and after the application of the first component of the SSP decomposition. This is illustrated in the next section.
Method to calculate the projector:
- PCA: What was described until now: SVD decomposition to extract spatial components.
- Average: Uses only one spatial component, the average of the time samples at which the select events occur. This has no effect if there are no events selected.

Advanced

Averaged artifact

One efficient way of representing the impact of this artifact correction is to epoch the recordings around the artifacts before and after the correction and compute the average of those epochs.

Run the process "SSP: Generic" with:
- The default blink options: Blink event, [-200,+200]ms, [1.5-15]Hz.
- The option "Computing existing SSP/ICA projectors" disabled.
- The option panel should look like the screen capture in the previous section.
Look at the topography of the first component. You can notice that the percentage value is higher than what we got previously, and that the topography looks cleaner than previously.
This difference comes from the fact that this time we did not use the cardiac SSP to compute the blink SSP ("Compute existing SSP" disabled). This leads to think that there is some form of cross-contamination of the "blink" and "cardiac" events that we designed here. The origin of the common signals between the different segments of artifact is sometimes due to important alpha waves (around 10Hz) that are present for most of the recordings. It don't matter much, you just have to remember that the computation order matters and that you can try variations of the suggested workflow to fit better your recordings.
You should see now two additional files in your database. They are both the average of the 19 blinks identified in the recordings, [-200,+200]ms around the "blink" events. The top row shows the average before the SSP correction, the bottom row the same average but recomputed after removing the first component of the decomposition. The artifact is gone.
Delete this new category, and make sure you get back to the previous settings (first component of both "cardiac" and "blink" selected). Click on save to validate this modification.

Advanced

Troubleshooting

You have calculated your SSP projectors as indicated here but you don't get any good results. No matter what you do, the topographies don't look like the targeted artifact. You can try the following:

Review one by one the events indicating the artifacts, remove the ones that are less clear or that occur close to another artifact.
Select or unselect the option "Compute using existing SSP".
Change the order in which you compute the projectors.
Use the process "Compute SSP: Generic" and modify some parameters:
- Use a narrower frequency band: especially the EOG, if the projectors capture some of the alpha oscillations, you can limit the frequency band to [1.5 - 9] Hz.
- Reduce or increase the time window around the peak of the artifact.
- Change the method: Average / SSP.
If you have multiple acquisition runs, you may try to use all the artifacts from all the runs rather than processing the files one by one. For that, use the Process2 tab instead of Process1. Put the "Link to raw file" of all the runs on both side, Files A (what is used to compute the SSP) and Files B (where the SSP are applied).

Always look at what this procedure gives you in output. Most of the time, the artifact cleaning will be an iterative process where you will need several experiments to adjust the options and the order of the different steps in order to get optimal results.

Advanced

SSP Theory

The Signal-Space Projection (SSP) is one approach to rejection of external disturbances. Here is a short description of the method by Matti Hämäläinen, from the MNE 2.7 reference manual, section 4.16.

Unlike many other noise-cancellation approaches, SSP does not require additional reference sensors to record the disturbance fields. Instead, SSP relies on the fact that the magnetic field distributions generated by the sources in the brain have spatial distributions sufficiently different from those generated by external noise sources. Furthermore, it is implicitly assumed that the linear space spanned by the significant external noise patterns has a low dimension.

Without loss of generality we can always decompose any n-channel measurement b(t) into its signal and noise components as:

b(t) = b_s(t) + b_n(t)

Further, if we know that b_n(t) is well characterized by a few field patterns b₁...b_m, we can express the disturbance as

b_n(t) = Uc_n(t) + e(t) ,

where the columns of U constitute an orthonormal basis for b₁...b_m, c_n(t) is an m-component column vector, and the error term e(t) is small and does not exhibit any consistent spatial distributions over time, i.e., C_e = E{ee^T} = I. Subsequently, we will call the column space of U the noise subspace. The basic idea of SSP is that we can actually find a small basis set b₁...b_m such that the conditions described above are satisfied. We can now construct the orthogonal complement operator

P_⊥ = I - UU^T

and apply it to b(t) yielding

b(t) ≈ P_⊥b_s(t) ,

since P_⊥b_n(t) = P_⊥Uc_n(t) ≈ 0. The projection operator P_⊥ is called the signal-space projection operator and generally provides considerable rejection of noise, suppressing external disturbances by a factor of 10 or more. The effectiveness of SSP depends on two factors:

The basis set b₁...b_m should be able to characterize the disturbance field patterns completely and
The angles between the noise subspace space spanned by b₁...b_m and the signal vectors b_s(t) should be as close to π/2 as possible.

If the first requirement is not satisfied, some noise will leak through because P_⊥b_n(t) ≠ 0. If the any of the brain signal vectors b_s(t) is close to the noise subspace not only the noise but also the signal will be attenuated by the application of P_⊥ and, consequently, there might by little gain in signal-to-noise ratio.

Since the signal-space projection modifies the signal vectors originating in the brain, it is necessary to apply the projection to the forward solution in the course of inverse computations.

Advanced

References

For more information on the SSP method, please consult the following publications:

C. D. Tesche, M. A. Uusitalo, R. J. Ilmoniemi, M. Huotilainen, M. Kajola, and O. Salonen, "Signal-space projections of MEG data characterize both distributed and well-localized neuronal sources," Electroencephalogr Clin Neurophysiol, vol. 95, pp. 189-200, 1995.
M. A. Uusitalo and R. J. Ilmoniemi, "Signal-space projection method for separating MEG or EEG into components," Med Biol Eng Comput, vol. 35, pp. 135-40, 1997.

Advanced

SSP Algorithm

The logic of the SSP computation is the following:

Take a small time window around each marker to capture the full effect of the artifact, plus some clean brain signals before and after. The default time window is [-200,+200]ms for eye blinks, and [-40,+40]ms for the heartbeats.
Filter the signals in a frequency band of interest, in which the artifact is the most visible.
Concatenate all those time blocks into a big matrix A = [b₁, ..., b_m]
Compute the singular value decomposition of this matrix A: [U,S,V] = svd(A, 'econ')
The singular vectors Ui with the highest singular values Si are an orthonormal basis of the artifact subspace that we want to subtract from the recordings. The software selects by default the vectors with eigenvalues above a certain threshold: Si > 12% * sum(Si)
Then it is possible to redefine interactively the selected components.
Calculate the projection operator: P_⊥i = I - UiUi^T
Apply this projection on the MEG or EEG recordings F: F = P_⊥iF
The process has to be repeated separately several times for each sensor type and each artifact.

Steps #1 to #5 are done automatically by the processes "SSP > Compute SSP" in the Record tab: the results, the vectors Ui, are saved in the channel file (field ChannelMat.Projector(i).Components).

Steps #6 and #7 are calculated on the fly when reading a block of recordings from the continuous file: when using the raw viewer, running a process a process on the continuous file, or importing epochs in the database.

Step #8 is the manual control of the process. Take some time to understand what you are trying to remove and how to do it. Never trust blindly any fully automated artifact cleaning algorithm, always check manually what is removed from the recordings, and do not give up if the first results are not satisfying.

Advanced

On the hard drive

The projectors are saved in the channel file associated with the recordings. This means that they will be shared by all the files that share the same channel file. As a consequence, you cannot share the channel files between acquisition runs if you are planning to use different SSP projectors for different runs.

You can find them in the field ChannelMat.Projector (array of structures):

Comment: String representing the projector in the window "Select active projectors".
Components: [Nsensors x Ncomponents], each column is one spatial component.
CompMask: [1 x Ncomponents], Indicates if each component is selected or not (0 or 1).
Status: 0=Category not selected, 1=Category selected, 2=Projectors already applied to the file.
SingVal: [1 x Ncomponents], Singular values of the SVD decomposition for each component.

Advanced

Additional documentation

Forum: Refining the pre-processing pipeline
McGill guidelines
SSP cookbook
Examples of topographies

-  ⇤ ← Revision 6 as of 2015-02-04 00:07:58 → 
  Size: 37379
  Editor: FrancoisTadel
  Comment:
+   ← Revision 52 as of 2015-07-06 22:44:25 → ⇥
  Size: 29003
  Editor: FrancoisTadel
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 1:
-= Tutorial 9: Artifacts & Signal Space Projections (SSP) =
+= Tutorial 12: Signal Space Projections (SSP) =
 Line 4:
+As previously said, the frequency filters are not adapted to remove artifacts that are transient or overlapping in frequency domain with the brain signals of interest. Other approaches exist to correct for those artifacts, based on the  spatial signature of the artifacts.

If an event is very reproducible and  occurs always at the same position (eg. eye blinks and heartbeats), the  sensors will always record the same values when it occurs. We can  identify the topographies corresponding to this artifact (ie. the spatial  distributions of values at one time point) and remove them from the  recordings. This spatial decomposition is the basic idea behind two  widely used approaches: the '''SSP '''(Signal-Space Projection) and '''ICA '''(Independent Component Analysis) methods.

This introduction tutorial will focus on the SSP approach, as it is a lot simpler and faster but still very efficient for removing blinks and heartbeats from MEG and high-density EEG recordings. The interface for running ICA decompositions is very similar and will be described in an advanced tutorial.
-Line 6:
+Line 12:
-= From Auditory =
== Heartbeats and eye blinks ==
 * Select the two AEF runs in the Process1 box.
 * Select successively the following processes, then click on [Run]:
  * '''Events > Detect  heartbeats:''' Select channel '''ECG''', check "All file", event name "cardiac".
  * '''Events > Detect eye blinks:''' Select channel '''VEOG''', check "All file", event name "blink".
  * '''Events > Remove simultaneous''': Remove "'''cardiac'''", too close to "'''blink'''", delay '''250ms'''.
  * '''Compute SSP: Heartbeats''': Event name "cardiac", sensors="MEG", '''do not use existing SSP'''.
  * '''Compute SSP: Eye blinks''': Event name "blink", sensors="MEG", '''do not use existing SSP'''. {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=ssp_pipeline.gif|ssp_pipeline.gif|height="359",width="527",class="attachment"}}
 * Double-click on '''Run01 '''to open the MEG. You can change the color of events "standard" and "deviant" to make the figure more readable.
 * Review the '''EOG '''and '''ECG '''channels and make sure the events detected make sense.
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=events.gif|events.gif|height="369",width="557",class="attachment"}}
 * In the Record tab, menu SSP > '''Select active projectors'''.
  * Blink: The first component is selected and looks good.
  * Cardiac: The category is disabled because no component has a value superior to 12%.
  * Select the first component of the cardiac category and display its topography.
  * It looks exactly like a cardiac topography, keep it selected and click on [Save].
  {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=ssp_result.gif|ssp_result.gif|height="172",width="700",class="attachment"}}
 * Repeat the same operations for '''Run02''':
  * Review the events.
  * Select the first cardiac component.

== Saccades ==
 * Run02 contains a few  saccades that generate a large amount of noise in the MEG recordings.  They are not identified well by the automatic detection process based on  the horizontal EOG. We have marked some of them, you have already  loaded those events together with the bad segments. We are going to use  again the SSP technique to remove the spatial components associated with  those saccades.
 * Open the MEG recordings for '''Run02''' and select the right-frontal sensors (Record tab > CTF RF).
 * In the Record tab, menu SSP > Compute SSP: Generic Event name='''saccade''', Time='''[0,500]ms''', Frequency='''[1,15]Hz''', '''Use existing SSP'''
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=ssp_saccade_process.gif|ssp_saccade_process.gif|height="442",width="329",class="attachment"}}
 * Example of saccade without correction:
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=ssp_saccade_before.gif|ssp_saccade_before.gif|height="259",width="711",class="attachment"}}
 * With the first component of saccade SSP applied:
  . {{http://neuroimage.usc.edu/brainstorm/Tutorials/Auditory?action=AttachFile&do=get&target=ssp_saccade_after.gif|ssp_saccade_after.gif|height="283",width="507",class="attachment"}}
 * This first component removes really well the saccade, keep it selected and click on [Save].

= From Continuous =
== Signal Space Projection ==
+== Overview ==
The general SSP objective is to identify the sensor toporagraphies that are typical of a specific artifact, then to create spatial projectors to remove the contributions of these topographies from the recordings.

 1. We start by identifying many examples of the artifact we are trying to remove. This is what we've been doing in the previous tutorial with the creation of the "cardiac" and "blink" events.
 1. We extract a short time window around each of these event markers and concatenate in time all those small blocks of recordings.
 1. We run a principle components analysis (PCA) on the concatenated artifacts in order to get a  decomposition in various spatial components (number of components = number of sensors).
 1. If it works well, we can find in the first few principal components some topographies that are very specific of the type of artifact we are targeting. We select those components to remove.
 1. We compute a linear projector for each spatial component to remove and save them in the database (in the "Link to raw file"). They are not applied immediately to the recordings.
 1. Whenever some recordings are read from this file, the SSP projectors are applied on the fly to remove the artifact contributions. This approach is fast and memory efficient.<<BR>><<BR>> {{attachment:ssp_intro.gif}}

== The order matters ==
This procedure has to be repeated separately for each artifact type. The order in which you process the artifacts matters, because for removing the  second artifact we typically use the recordings cleaned with  the first set of SSP projectors. We have to decide which one to process first.

It works best if each artifact is defined precisely and as independently as possible from the other artifacts. If the two artifacts happen simultaneously,  the SSP projectors  calculated for the blink may contain some of    the  heartbeat topography and vice  versa. When trying to remove the second  artifact,    we might not be  able to isolate it clearly anymore.

Because the heart beats every second or so, there is a high chance that  when the subject blinks there is a heartbeat not too far away in the  recordings. Therefore a   significant  number of the blinks will be contaminated with  heartbeats. But we have usually a lot of "clean" heartbeats, we can start by removing those ones.   To isolate correctly these two common artifacts, we   recommend the following procedure:

 * Remove the markers "cardiac" that are occurring during a blink (done in the previous tutorial),
 * Compute the '''cardiac SSP''' (with no eye movements, because we removed the co-occurring events),
 * Compute the '''blink SSP''' (with no heartbeats, because they've already been taken care of).

If you have multiple modalities recorded simultaneously, for example MEG and EEG, you should run twice this entire procedure, once for  the EEG only  and once for the MEG only. You will always  get better results if  you  '''process the  different types of sensors  separately'''. Same  thing  when processing  Elekta-Neuromag  recordings: process   separately the  magnetometers (MEG  MAG) and the gradiometers (MEG GRAD).

== SSP: Heartbeats ==
Double-click on the link to show the MEG sensors for '''Run #01'''.<<BR>>In the Record tab, select the menu: '''"Artifacts > SSP: Heartbeats"'''.

 * '''Event name''': Name of the event to use to calculate the projectors, enter "'''cardiac'''".
 * '''Sensor types''': Type of sensors for which the projection should be calculated ("'''MEG'''").    Note that you will always  get better results if  you  process the different types of sensors  separately.
 * '''Compute using existing SSP projectors''':    You have the option to calculate the projectors from the raw   recordings, or from the recordings filtered with the previously    computed SSP projectors. Unless you have a good reason for not   considering the existing projectors, you should select this   option. <<BR>><<BR>> {{attachment:ssp_ecg_process.gif||height="267",width="501"}}

After   the computation is done, a new figure is  displayed, that lets you   select the active projectors.

 * '''On the left''': The projector   categories where each row  represents the  result of  an execution of this process  (usually one for  each sensor  type and each artifact).
 * '''On the right''': The spatial components returned by the PCA decomposition. The percentage indicated between brackets is the  singular value for this each component, normalized  for  this decomposition  (percentage = '''S'''i / sum('''S'''i), see technical details at the end of this page).
 * '''Percentage''': More practically, it indicates  the amount of signal that was captured by the component during the decomposition. The highest it is, the more the component is representative of the  artifact recordings that were used to    calculate it. In the good  cases, you would typically see one to three    components with values  that are significantly higher than the others.
 * When a  component is selected, it means that it is removed from the  recordings. A spatial projector is computed and applied to the   recordings on  the fly  when reading from the continuous file.
 * '''Default selection''': The software selects '''the components with values superior to  12%''' and leaves the others unselected. '''This threshold is completely  empirical''' '''and depends on the acquisition device''', you should always  review manually the components that you want to remove. <<BR>><<BR>> {{attachment:ssp_select1.gif||height="178",width="366"}}

== Evaluate the components ==
For   the cardiac event, none of the components show a value superior to  12%,  therefore  the entire "cardiac" projector category was unselected.  This  doesn't  mean that none of those components is actually  interesting for  us. The percentage indicated for the first value (9%) is much higher than the following ones (5%, 5%, 4%, 3%...), this could indicate that it targets relatively well the cardiac artifact. Let's  investigate this in more details.

 * Select the category "cardiac" to enable the list on the   right.
 * Click on the first component, then click on the toolbar button '''[Display component topography]'''. This menu shows the spatial distribution of the sensor values for this component. Note that you don't have to select the component  (ie. check the box) to display it. This topography seems to correspond to a strong dipolar activity located relatively far from the sensor array, it matches the type of artifact we expect from the heart activity. <<BR>><<BR>> {{attachment:ssp_select2.gif||height="143",width="436"}}
 * The second button "'''Display component topography [No magnetic interpolation]'''" produces the same figure but without the reinterpolation of the magnetic fields that is typically applied to the MEG recordings in Brainstorm, it may help understand some difficult cases. This magnetic interpolation will be detailed later in the introduction tutorials.
 * You can display multiple components in the same figure: select them at the same time in the list (holding the Shift/Ctrl/Cmd button of your keyboard) and then click on the button "Display topography".<<BR>>No other strong components looks like it could be related with the heartbeats. <<BR>><<BR>> {{attachment:ssp_multiple.gif||height="197",width="286"}}
 * The last button in the toolbar '''[Display component time series]''', opens a figure that represents the evolution of the contribution of this component over time. The higher the amplitude, the more present the selected topography in the recordings. Click on it to show the''' component #1''', then display the '''ECG '''signal at the same time (right-click on the file > ECG > Display time series).   <<BR>><<BR>> {{attachment:ssp_cardiac_good.gif||height="158",width="361"}}
 * We observe that the "SSP1" trace correlates relatively well with the ECG  trace, in the sense that we captured most the ECG peaks with this  component. However, the component seems also to capture much more signal than just the heartbeats: many alpha oscillations and some of the ocular activity. The example below shows a blink in the EOG, ECG and SSP component #1. <<BR>><<BR>> {{attachment:ssp_cardiac_bad.gif||height="208",width="361"}}

 * If you remove this component from the recordings, you can expect to see most of the artifacts related with '''the cardiac activity to go away''', but you will also '''remove additional signal''' elements that were not really well identified. The job is done but it causes some unwanted side effects.

 * It is in general possible to refine the SSP decomposition by going back to the selection of "cardiac" markers that we used to compute it. You could look at all the ECG peaks individually and remove the markers located in segments of recordings that are noisier or that contain a lot of alpha activity (~10Hz). You would need to delete this SSP decomposition and run again the same process.

 * Alternatively, or if you don't manage to extract a clean cardiac component with a PCA/SSP decomposition, you could try to run an ICA decomposition instead. Note that for some subjects, the cardiac artifact is not very strong and could be simply ignored in the analysis.

== Evaluate the correction ==
The topography of the component #1 looks like it represents a cardiac topography and its temporal evolution shows peaks where we identified heartbeats. It is therefore a good candidate for removal, we just need to make sure the signals look good after the correction before validating this choice.

 * Show the left-temporal MEG sensors (CTF LT) and select/unselect the first SSP component.
 * Repeat this for different time windows, to make sure that the cardiac peaks in the MEG sensors really disappear when the projector #1 is selected and that the rest is not altered too much.
 * No correction: <<BR>><<BR>> {{attachment:ssp_ecg_off.gif||height="235",width="550"}}
 * Cardiac component #1 removed: <<BR>><<BR>> {{attachment:ssp_ecg_on.gif||height="235",width="549"}}

 * In this example we will consider that the current decomposition is good enough. <<BR>>Make sure you '''select the component #1''', then click on '''[Save]''' to validate the modifications.
 * After this window is closed, you can always open it again from the Record tab with the menu '''Artifacts > Select active projectors'''. At this stage of the analysis, you can modify the list of projectors applied to the recordings at any time.

== SSP: Eye blinks ==
Let's try the same thing with the eye blinks.

 * Select the process "'''SSP > Compute SSP: Eye blinks'''"
 * Run it on the event type "blink", that indicates the peaks of the EOG  signal. <<BR>><<BR>> {{attachment:ssp_eog_process.gif||height="268",width="501"}}
 * You see now a new category of projectors. It is selected by default because its first component shows a percentage higher than 12% (empirical and somewhat random threshold). Based on the distribution of values, this first component is most likely a good representation of the    artifact we are trying to remove. The second one could be a good candidate as well. <<BR>><<BR>> {{attachment:ssp_select_eog1.gif||height="169",width="329"}}
 * Select the first three components and display their topographies:<<BR>><<BR>> {{attachment:ssp_eog_topo.gif||height="158",width="300"}}
  * Component #1: Most likely a '''blink''',
  * Component #2: Probably a '''saccade''' (other type of eye movement we were not targeting),
  * Component #3: Not related with the eye movements (maybe related with the alpha activity).
 * Display the times series for those three components, together with the EOG signals. You have to uncheck temporarily the component #1 to be able to display its signal. When it is checked, it is removed from the signal therefore it is corresponds to a flat trace. <<BR>>The figure below shows the EOG and SSP values between 318s and 324s. The SSP1 trace matches the blink observed in VEOG and SSP2 matches the saccade observed in HEOG.<<BR>><<BR>> {{attachment:ssp_eog_ts.gif||height="157",width="367"}}
 * Observe the activity when there is no component selected: <<BR>><<BR>> {{attachment:ssp_eog_off.gif||height="205",width="550"}}
 * With only the component #2 selected (saccade): <<BR>><<BR>> {{attachment:ssp_eog_on1.gif||height="205",width="551"}}
 * With components #1 and #2 selected (blink + saccade): <<BR>><<BR>> {{attachment:ssp_eog_on2.gif||height="205",width="551"}}
 * Keep the '''components #1 and #2 selected''' and click on [Save] to validate your changes.

== Run #02 ==
Reproduce the same operations on '''Run #02''':

 * Close everything with the [X] button at the top-right corner of the Brainstorm window.
 * Open the MEG recordings for run '''AEF #02''' (double-click on the file link).
 * '''Artifacts > SSP: Heartbeats:''' Event name "cardiac", sensors "MEG. <<BR>>Select '''component #1''', click on [Save]. <<BR>><<BR>> {{attachment:ssp_run02_ecg.gif||height="155",width="450"}}
 * '''Artifacts > SSP: Heartbeats:''' Event name "cardiac", sensors "MEG. <<BR>>Select '''component #1''', click on [Save]. <<BR>><<BR>> {{attachment:ssp_run02_eog.gif||height="155",width="450"}}

 * Note that in this second session, the representation of the saccade was not as clear as in the first file. The distribution of the percentage values does not show any clear component other from the blink one, and the topographies are not as clear. In general, the saccade processing requires a separate step, we will illustrate this in the next tutorial.

<<TAG(Advanced)>>

== SSP: Generic ==
The calculation of the SSP for the heartbeats and the eye blinks are shortcuts to a more generic algorithm "'''Compute SSP: Generic'''".   You may need this process if the standard parameters do not work of if you want to use this technique to remove other types of artifacts.

 * '''Time window''': What segment of the file do you want to consider.
 * '''Event name''': Markers that are used to characterize the artifact. If you don't specify any event name, it will use this entire time window as the input of the PCA decomposition.
 * '''Event window''':   Time segment to consider before and after each event marker. We want   this time window to be longer than the artifact effect itself, in order   to have a large number of time samples representative of a normal brain   activity. This helps the PCA decomposition to separate the artifact  from  the ongoing brain activity.
 * '''Frequency band''':   Definition of the band-pass filter that is applied to the recordings   before calculating the projector. Usually you would use the same   frequency band as we used for the detection, but you may want to try to   refine this parameter if the results are not satisfying.
 * '''Sensor types or names''': List of sensors for which to calculate the projector. You can get better results if you process one sensor type at a time.
 * '''Compute using existing SSP/ICA projectors''':  Same as in the heartbeats/blinks processes.
 * '''Save averaged artifact in the database''': If you check this option with an event name selected, the process will save the average of all the artifact epochs (event marker + event window) before and after the application of the first component of the SSP decomposition. This is illustrated in the next section.
 * '''Method to calculate the projector''':
  * '''PCA''': What was described until now: SVD decomposition to extract spatial components.
  * '''Average''': Uses only one spatial component, the average of the time samples at which the select events occur. This has no effect if there are no events selected. <<BR>><<BR>> {{attachment:ssp_generic.gif||height="429",width="316"}}

<<TAG(Advanced)>>

== Averaged artifact ==
One   efficient way of representing the impact of this artifact correction  is  to epoch the recordings around the artifacts before and after the   correction and compute the average of those epochs.

 * Run the process "SSP: Generic" with:
  * The default blink options: Blink event, [-200,+200]ms, [1.5-15]Hz.
  * The option "Computing existing SSP/ICA projectors" '''disabled'''.
  * The option panel should look like the screen capture in the previous section.
 * Look at the topography of the first component. You can notice that the percentage value is higher than what we got previously, and that the topography looks cleaner than previously.  <<BR>><<BR>> {{attachment:ssp_generic_topo.gif||height="181",width="499"}}
 * This difference comes from the fact that this time we did not use the cardiac SSP to compute the blink SSP  ("Compute existing SSP" disabled). This leads to think that there is some form of cross-contamination of the "blink" and "cardiac" events that we designed here. The origin of the common signals between the different segments of artifact is sometimes due to important alpha waves (around 10Hz) that are present for most of the recordings. It don't matter much, you just have to remember that the computation order matters and that you can try variations of the suggested workflow to fit better your recordings.
 * You should see now two additional files in your database. They are both the average of the 19 blinks identified in the recordings, [-200,+200]ms around the "blink" events. The top row shows the average before the SSP correction, the bottom row the same average but recomputed after removing the first component of the decomposition.  The artifact is gone. <<BR>><<BR>> {{attachment:ssp_eog_average.gif||height="277",width="599"}}
 * Delete this new category, and make sure you get back to the previous settings (first component of both "cardiac" and "blink" selected). Click on save to validate this modification.

<<TAG(Advanced)>>

== Troubleshooting ==
You    have calculated your SSP projectors as indicated here but you don't get any good    results. No matter what you do, the topographies don't look like the    targeted artifact. You can try the following:

 * Review one by    one the events indicating the artifacts, remove the ones that are  less   clear or that occur close to another artifact.
 * Select or unselect the option "Compute using existing SSP".
 * Change the order in which you compute the projectors.
 * Use the process "Compute SSP: Generic" and modify some parameters:
  * Use    a narrower frequency band: especially the EOG, if the projectors     capture some of the alpha oscillations, you can limit the frequency  band    to [1.5 - 9] Hz.
  * Reduce or increase the time window around the peak of the artifact.
  * Change the method: Average / SSP.
 * If    you have multiple acquisition runs, you may try to use all the    artifacts from all the runs rather than processing the files one by one.    For that, use the Process2 tab instead of Process1. Put the "Link to    raw file" of all the runs on both side, Files A (what is used to  compute   the SSP) and Files B (where the SSP are applied).

Always   look at what this procedure gives you in output. Most of the time, the  artifact  cleaning will  be an iterative process where you will need several experiments to  adjust the options  and the order of the different steps in  order to get optimal results.

<<TAG(Advanced)>>

== SSP Theory ==
-Line 70:
+Line 181:
+<<TAG(Advanced)>>

== References ==
-Line 72:
+Line 186:
- * C.  D. Tesche, M. A. Uusitalo, R. J. Ilmoniemi, M. Huotilainen, M. Kajola,  and O. Salonen, "Signal-space projections of MEG data characterize both  distributed and well-localized neuronal sources," Electroencephalogr  Clin Neurophysiol, vol. 95, pp. 189-200, 1995.
 * M. A. Uusitalo  and R. J. Ilmoniemi, "Signal-space projection method for separating MEG  or EEG into components," Med Biol Eng Comput, vol. 35, pp. 135-40, 1997.

== Identify the artifacts ==
The first step is to identify several repetitions of the artifact (the vectors '''b'''<<HTML(<sub>)>>''1''<<HTML(</sub>)>>...'''b'''<<HTML(<sub>)>>''m''<<HTML(</sub>)>>).  We need to set markers in the recordings that indicate when the events  that we want to correct for occur. To help with this task, it is  recommended to always record with bipolar electrodes the activity of the  eyes (electro-oculogram or EOG, vertical and horizontal), the heart  (electro-cardiogram or ECG), and possibly other sources of muscular  contaminations (electromyogram or EMG). In this example, we are going to  use the ECG and vertical EOG traces to mark the cardiac activity and  the eye blinks.  Two methods can be used, manual or automatic.

=== EOG/ECG channels ===
 * Select the protocol '''TutorialRaw'''  created in the previous tutorial, and select the "Functional data" view  (second button in the toolbar on top of the database explorer).
 * Double-click on the '''clean''' continuous recordings ("'''Raw | notch(60Hz 120Hz 180Hz)'''") to open the MEG recordings.
 * Set the length of the reviewed time window to '''3 seconds''' (in the Record tab, text box "Duration")
 * Right-click on the continuous recordings again > '''Misc > Display time series'''.

  {{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=tsMisc.gif|tsMisc.gif|height="148",width="577",class="attachment"}}
 * In this file the channel type "'''Misc'''" groups the two channels '''EEG057 '''('''ECG''', in green) and '''EEG058 '''('''vertical EOG''',  in red). This configuration depends on the acquisition setup, and can  be redefined afterwards in Brainstorm (right-click on the channel file  > Edit channel file, and then change manually the string in the  column Type for any channel).
 * Use  the shortcuts introduced in the previous tutorial to adjust the  vertical scale of this display: Shift+mouse wheel, +/- keys, or the  buttons on the right side of the figure.
 * Then go further in time to see what is happening on those channels over the time: use the "'''>>>'''" buttons, or the associated shortcuts (read the tooltips of the button)
 * '''ECG''':  On the green trace, you can recognize the very typical shape of the  electric activity of the heart (P, QRS and T waves). This is a very good  example, the signal is not always that clean.
 * '''EOG''':  On the red trace, there is not much happening for most of the  recordings except for a few bumps, typical of eye blinks, eg. at  33.590s. This subject is sitting very still and not blinking much. We  can expect MEG recordings of a very good quality.
 * You can observe the contamination from a blink on the left-frontal sensors: move to '''33.590s'''  (you can use the text boxes in the time panel, the Record tab, or the  scrollbar in the figure), and select a subset of sensors (Shift+B, or  right-click on the figure > Display setup > CTF LF)

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

=== Manual marking ===
Create a new category of markers "'''blink_manual'''",  using the menu Events > Add group. Select this new group. Review the  file, and mark the peaks you observe on the vertical EOG trace, using  the '''Ctrl+E''' keyboard shortcut. Do that for a few eye blinks.

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

You  could repeat the same operation for all the blinks, then for all the  ECG peaks and jump to the next chapter of the tutorial and compute the  SSP. It would be uselessly time consuming, as there is a process that  does it for you automatically. However, it is good to remember how to do  it manually because you may face some cases where you don't have clean  ECG/EOG, or if you want to correct for another type of artifact.

=== Automatic detection: EOG ===
In the Record tab, select the menu: '''SSP > Detect eye blinks'''. It opens automatically the pipeline editor, with the process "Detect eye blinks" selected:

 * '''Channel name''': Name of the channel that is used to perform the detection. Select or type "'''EEG058'''" as it is the name of the EOG channel
 * '''Time window''':  Time range that the algorithm should scan for the selected artifact.  Leave the default values to process the entire file.
 * '''Event name''': Name of the event group that is created for saving all the detected events. Leave the default "blink".

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

Click on Run. After the process stops, you can see two new event categories "'''blink'''" and "'''blink2'''"     in the Record tab. You can review a few of them,   to make sure that   they really indicate the EOG events. In the Record tab, click on the   "blink" event category, then on a time occurrence to   jump to it in the   MEG and Misc time series figures.

Two  types of  events are created because this algorithm not only detects  specific  events in a signal, it also classifies them by shape. If you  go through  all the events that were detected in the two categories, you  would see  that the "blink" are all round bumps, typical of the '''eye blinks'''.   In the category "blink2", the morphologies don't look as uniform; it   mixes small blinks, and ramps or step functions followed by sharp drops   that could indicate '''eye saccades'''. The saccades can be  observed on  the vertical EOG, but if you want a better characterization  of them you  should also record the horizontal EOG. The detection of  the saccades  should be performed with a different set of parameters,  using the  process "Detect custom events", introduced later in this  chapter.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=detectEogDone.gif|detectEogDone.gif|class="attachment"}}

=== Automatic detection: ECG ===
Now do the same thing for the heartbeats. In the Record tab, select the  menu "'''SSP > Detect heartbeats'''". Configure the process to use the  channel '''EEG057''' (name of the ECG channel), and leave the other options to the default values.

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

Click on Run. After the process stops, you can see a new event category "'''cardiac'''"   in the Record tab, with 346 occurrences. You can check a few of them,   to make sure that the "cardiac" markers really indicate the ECG peaks,   and that there are not too many peaks that are skipped.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=detectEcgDone.gif|detectEcgDone.gif|class="attachment"}}

=== Automatic detection: Custom ===
Those two previous processes are shortcuts for a generic process "'''Detect custom events'''".  We are not going to use it here, but it is interesting to introduce it  to understand how the blinks and heartbeats detection work. The logic is  the following:

 * The channel to analyze is read from the continuous file, for a given time window.
 * '''Frequency band''': The signal is filtered in a frequency band where the artifact is easy to detect. For EOG: 1.5-15Hz ; for ECG: 10-40Hz.
 * '''Threshold''':  An event of interest is detected if the absolute value of the filtered   signal value goes over a given number of times the standard deviation.  For EOG: 2xStd, for ECG: 4xStd
 * '''Minimum duration between two events''':  If the filtered signal crosses the threshold several times in relation  with the same artifact (like it would be the case for muscular activity  recordings on an EMG channel), we don't want to trigger several events  but just one at the beginning of the activity. This parameter would  indicate the algorithm to take only the maximum value over the given  time window; it also prevents from detecting other events immediately  after a successful detection.
 For the ECG, this value is set to 500ms, because it is very unlikely  that the heart rate of the subject goes over 120 beats per minute.
 * '''Ignore the noisy segments''':  If this option is selected, the detection is not performed on the  segments that are much noisier than the rest of the recordings.
 * '''Enable classification''':  If this option is selected, the events are classified by shape, based  on correlation measure. In the end, only the categories that have more  than 5 occurrences are kept, all the other successful detections are  ignored.

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

== Calculate the projectors ==
=== Methodology ===
We  have now 346 examples of heartbeats, 18 examples of blinks, and 12  examples of saccades. Those are sufficient numbers of repetitions to  calculate reliable signal-space projections. The logic of the  computation is the following:
+ * C.    D. Tesche, M. A. Uusitalo, R. J. Ilmoniemi, M. Huotilainen, M.  Kajola,   and O. Salonen, "Signal-space projections of MEG data  characterize  both  distributed and well-localized neuronal sources,"   Electroencephalogr  Clin Neurophysiol, vol. 95, pp. 189-200, 1995.
 * M.   A. Uusitalo  and R. J. Ilmoniemi, "Signal-space projection method for   separating MEG  or EEG into components," Med Biol Eng Comput, vol. 35,   pp. 135-40, 1997.

<<TAG(Advanced)>>

== SSP Algorithm ==
The logic of the SSP  computation is the following:
-Line 143:
+Line 195:
-. Filter the MEG or EEG signals in a frequency band of interest, in which the artifact is the most visible.
+. Filter the signals in a frequency band of interest, in which the artifact is the most visible.
-Line 149:
+Line 201:
-. The process has to be repeated '''separately '''several times:
  * for each sensor type (EEG, MEG magnetometers, MEG gradiometers) and
  * for  each artifact, starting with the one that leads to the strongest  contamination in amplitude on the sensors (in the case: the eye blinks)
  * If  you have difficulties removing one artifact or the other, you may try  to process them in a different order. You may also try removing some of  the artifact markers in the case of co-occurring artifacts. If a lot of  the blinks happen at the same time as heartbeats, you may end up  calculating projectors that mix both effects, but that do not remove  efficiently one or the other. In this case, remove all the markers that  happen in the segments contaminated by multiple artifacts.

Steps #1 to #5 are done automatically by the processes "SSP > Compute SSP" in the Record tab: the results, the vectors '''U'''i, are saved in the database link ("Link to raw file") and in the channel file.
+. The process has to be repeated '''separately '''several times for each sensor type and each artifact.

Steps #1 to #5 are done automatically by the processes "SSP > Compute SSP" in the Record tab: the results, the vectors '''U'''i, are saved in the channel file (field ChannelMat.Projector(i).Components).
-Line 160:
+Line 209:
-=== Recommended procedure ===
We  have multiple types of artifact to correct for and it can be tricky to  know what artifact to correct first. If you correct for artifact A and  then for artifact B, you may get different results than if you correct  first B and then A. We propose here a default order for the two  artifacts that we have in all the recordings: the blinks and the  heartbeats.

If two artifacts happen at the same  time, the SSP projectors calculated for the blink may contain some of  the heartbeat and vice versa. When trying to remove the second artifact,  we might not be able to isolate it clearly anymore. To avoid this  behavior, we should try to avoid having overlapping artifacts. Because  the subject's heart is beating at least once per second, there is a high  chance to record a heartbeat next to any eye blink. Therefore a  significant number of the blinks will be contaminated with heartbeats.  If we decide to remove all the blink events that include some heart  activity, we may not have enough blink samples left to build correct SSP  projectors. To isolate correctly these two common contaminations we  recommend the following procedure:

 1. Remove all the  cardiac peaks that are occurring next to a blink (we can remove half of  the heartbeats, we will still have enough data for calculating the  projectors).
 1. Calculate the SSP to remove the '''cardiac '''artifact.
 1. Calculate the SSP to remove the '''eye blinks''', based on the recordings already corrected for the cardiac artifact.

=== Remove simultaneous events ===
If you closed the file viewer, open the MEG recordings again:
 double click on the clean file "Raw | notch(60Hz 120Hz 180Hz)".

In the Record tab, select the menu "'''SSP > Remove simultaneous'''". Set the options:

 * Remove events named: "'''cardiac'''"
 * When too close to events: "'''blink'''"
 * Minimum delay between events: '''250ms'''

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

After executing this process, the number of "cardiac" events goes from '''346''' to '''333'''. The deleted heartbeats were all less than 250ms away from a blink.

=== Compute SSP: Heartbeats ===
In the Record tab, select the menu "'''SSP > Compute SSP: Heartbeats'''". Set the options:

 * '''Event name''': Name of the event to use to calculate the projectors, enter "'''cardiac'''"
 * '''Sensor types''': Type of sensors for which the projection should be calculated ("'''MEG'''").   It doesn't matter if there are other sensor types indicated in the  text  box: the sensor names or types that do not exist in the recordings  are  just ignored. Note that you will always  get better results if you  process the different types of sensors  separately.
 If you have MEG and EEG in the same file, you should run twice this  process, once for the EEG only and once for the MEG only. Same  thing  when processing Elekta-Neuromag recordings, you should process   separately the magnetometers (MEG MAG) and the gradiometers (MEG GRAD).
 * '''Compute using existing SSP projectors''':   You have the option to calculate the projectors from the raw  recordings, or from the recordings filtered with the previously   computed SSP projectors. Unless you have a good reason for not  considering the existing SSP projectors, you should always select this  option.

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

After  the computation is done, a new figure is  displayed, that lets you  select the active projectors. On the left you  have the projector  categories (matrix '''U''') where each row  represents the  result of an execution of a "Compute SSP " process  (usually one for  each sensor type and each artifact). On the right, you  can activate  independently the different components for the selected  projector  category. Each entry represents the projector created from a  column  vector '''U'''i (singular vector of the SVD decomposition),  and  the percentage is the singular value for this vector, normalized  for  this decomposition (percentage = '''S'''i / sum('''S'''i)).  When a  projector is selected, it means that it is applied to the  recordings on  the fly when reading from the continuous file, ie. that  the  corresponding spatial component is removed from the recordings.

The   percentage indicates the amount of signal that was captured by this   projector during the decomposition. The highest it is, the more the   projector is representative of the artifact recordings that were used to   calculate it. In the good cases, you would typically see one to three   components with values that are significantly higher that the others.

By  default, the software selects '''the components with values superior to  12%''' and leaves all the others unselected. '''This threshold is completely  empirical''' '''and depends on the acquisition device''', you should always  review manually the components that you want to remove.

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

For  the cardiac event, none of the components show a value superior to 12%,  therefore  the entire "cardiac" projector category was unselected. This  doesn't  mean that none of those components is actually interesting for  us, let's  investigate this.

A list of buttons is available to manipulate the projector categories: load, save, rename, delete, and '''display component topography'''.   The last one is one of the most useful tool to understand those   results. Select the category ("cardiac") to enable the list on the  right, click on the first component, then click on the toolbar button to  display its topography. Do it  again for components #2 and #3 (the  components don't need to be  active for that).

The  other  button [No interp] displays the same results without the  magnetic field  smoothing that is applied by default to MEG  topographies. It may help  understand some difficult cases.

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

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=sspEcgComp.gif|sspEcgComp.gif|class="attachment"}}

None  of those topographies can be clearly related to a  cardiac component.  This can have two interpretations: the cardiac  artifact is not very  strong for this subject and  the influence of the  heart activity over  the MEG sensors is completely buried in the noise or  the brain signals,  or the characterization of the artifact that we did  was not so good  and we should refine it by selecting for instance  smaller time windows  around the cardiac peaks. Here, it's probably due  to the subject's  morphology. Some people generate strong artifacts in  the MEG, some  don't.

In this case, it is safer to '''unselect '''this "cardiac" category, rather than removing randomly spatial  components that we are not sure to identify.

=== Compute SSP: Eye blinks ===
Let's try the same thing with the eye blinks. Select the process "'''SSP > Compute SSP: Eye blinks'''".  Run it on the event type "blink", that indicate the peaks of the EOG  signal.

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

In   this specific case, there is one value that is much higher than the   others (21%), and it is selected by default because it is superior to   12%. The first component is most likely a good representation of the   artifact we are trying to remove.

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

Select  the cardiac category, and display the topographies for the first three  components. Component #1 is clearly related with the ocular  activity  occurring during the blinks, the others are not as clear.  Conclusion:  the default proposition was good (only component #1  selected).  Setting  this first projector as active removes from the MEG  signals the  contributions related with the spatial topography we can see  in the  first figure, ie. a spatial distribution of sensor values  related  specifically to the eye blinks.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=sspEogComp.gif|sspEogComp.gif|class="attachment"}}

While  the recordings are displayed in the file viewer, the selection of the  active projectors interactively updates the recordings. Go to the first  detected "blink" event ('''t=33.600s'''), and select the left frontal sensors ('''Shift+B'''  or right-click > Display setup > CTF LF). Then play with the  checkboxes to change the active projectors, and see how well this first  component of the blink category removes the ocular artifact.

{{http://neuroimage.usc.edu/brainstorm/Tutorials/TutRawSsp?action=AttachFile&do=get&target=sspEogMeg.gif|sspEogMeg.gif|class="attachment"}}

Select  the projector category "blink", select the first component, and click  on Save, to validate your changes. Starting from now, all the successive  operations you apply on this file will integrate this eye blink  correction, including the calculation of other projectors. After you  closed this window, you can get back to it by selecting the menu "'''SSP > Select active projector'''" in the Record tab, and change again the selection of active projectors.

We are not going to calculate the SSP for the other category of events identified on the EOG, the group "'''blink2'''".  The events classified in this group could represent eye saccades or  other types of blinks, it is not very clear that they correspond to the  same physiological event. Most of them do not seem to create any obvious  contamination on the MEG recordings, there is probably no need to  consider them now.

=== Compute SSP : Generic ===
The calculation of the SSP for the heartbeats and the eye blinks are shortcuts to a more generic algorithm "'''Compute SSP: Generic'''".  We do not need it here, but the illustration of its properties is  interesting to understand the correction of ocular and cardiac  artifacts:

 * '''Event name''': Markers that are used to characterize the artifact
 * '''Time window''':  Time segment to consider before and after each event marker. We want  this time window to be longer than the artifact effect itself, in order  to have a large number of time samples representative of a normal brain  activity. This helps the SVD decomposition to separate the artifact from  the ongoing brain activity.
 * '''Frequency band''':  Definition of the band-pass filter that is applied to the recordings  before calculating the projector. Usually you would use the same  frequency band as we used for the detection, but you may want to try to  refine this parameter if the results are not satisfying.
 * '''Sensor types or names''': List of sensors for which to calculate the projector. You can get better results if you process one sensor type at a time.

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

=== Troubleshooting ===
You  have calculated the SSP as indicated here but you don't get any good  results. No matter what you do, the topographies don't look like the  targeted artifact. You can try the following:

 * Review one by  one the events indicating the artifacts, remove the ones that are less  clear or that occur close to another artifact.
 * Select or unselect the option "Compute using existing SSP".
 * Use the process "Compute SSP: Generic" and modify some parameters:
  * Reduce the time window around the peak of the artifact.
  * Use  a narrower frequency band: especially the EOG, if the projectors  capture some of the alpha oscillations, you can limit the frequency band  to [1.5 - 9] Hz.
  * Change the method: Average / SSP.
 * If  you have multiple acquisition runs, you may try to use all the  artifacts from all the runs rather than processing the files one by one.  For that, use the Process2 tab instead of Process1. Put the "Link to  raw file" of all the runs on both side, Files A (what is used to compute  the SSP) and Files B (where the SSP are applied).

== Evaluation ==

=== Eye blinks ===
One  efficient way of representing the impact of this artifact correction is  to epoch the recordings around the artifacts before and after the  correction. Use the same time window as the one used in the SSP options  around each marker, and average all the segments of recordings. Those  operations are detailed in the next tutorial, we are just presenting the  results. You don't need to reproduce them at this time, just remember  that it is doable, in case you need it with your own data. This is what  the database would look like after those operations:

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

The following image represents the file '''"blink_uncorrected / Avg: blink"''', which is the average of the 18 identified blinks '''before'''  the SSP correction, [-200,+200]ms around the "blink" events. The time  series, the 2D topography at the peak (t = 0ms), and the mapping on the  cortex of those fields, using the basic inverse model calculated in the  previous tutorial.

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

The next image represents the file  '''"blink_ssp / Avg: blink"''', which is the exact same thing, but '''after '''the  SSP correction. The artifact is gone. If you do this and you can still  see some clear evoked component in the time series figure, the SSP  correction was not efficient: the artifact is not properly identified,  or you should select different components using the "Select active  projectors" window.

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

=== Heartbeats ===
We can do  the same thing with the heartbeats: epoch [-40,+40]ms around the   "cardiac" events, average the trials, and review the event-related  fields at three instants:

The peak of the '''P wave''' (t = '''-30ms'''):

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

The peak of the '''QRS complex''' (t = '''0ms'''):

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

The peak of the '''T wave''' (t = '''25ms'''):

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

Those  peaks may look big, but they are in fact much smaller (300 fT for the  average of 346 repetitions) than what we observed for the eye blinks  (1500 fT for 18 repetitions). You can notice that none of those  topographies are similar to the components we obtained after calculating  the SSP for the cardiac artifact. We don't know how to correct this  artifact, it doesn't look too bad in terms of recordings contamination,  so we just leave it uncorrected.

== Additional discussions on the forum ==
 * Refining the pre-processing pipeline: <<BR>>http://neuroimage.usc.edu/forums/showthread.php?1272

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

<<EmbedContent(http://neuroimage.usc.edu/bst/get_feedback.php?Tutorials/Epoching)>>
+<<TAG(Advanced)>>

== On the hard drive ==
The projectors are saved in the channel file associated with the recordings. This means that they will be shared by all the files that share the same channel file. As a consequence, you cannot share the channel files between acquisition runs if you are planning to use different SSP projectors for different runs.

You can find them in the field ChannelMat.Projector (array of structures):

 * '''Comment''': String representing the projector in the window "Select active projectors".
 * '''Components''': [Nsensors x Ncomponents], each column is one spatial component.
 * '''CompMask''': [1 x Ncomponents], Indicates if each component is selected or not (0 or 1).
 * '''Status''': 0=Category not selected, 1=Category selected, 2=Projectors already applied to the file.
 * '''SingVal''': [1 x Ncomponents], Singular values of the SVD decomposition for each component.

{{attachment:ssp_storage.gif||height="288",width="675"}}

<<TAG(Advanced)>>

== Additional documentation ==
 * Forum: [[http://neuroimage.usc.edu/forums/showthread.php?1272|Refining the pre-processing pipeline]]
 * McGill guidelines
 * SSP cookbook
 * Examples of topographies

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

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

Software

Users

Development

Tutorial 12: Signal Space Projections (SSP)

Overview

The order matters

SSP: Heartbeats

Evaluate the components

Evaluate the correction

SSP: Eye blinks

Run #02

SSP: Generic

Averaged artifact

Troubleshooting

SSP Theory

References

SSP Algorithm

On the hard drive

Additional documentation

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):