Tutorial 9: Artifacts & Signal Space Projections (SSP)

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

Contents

  1. Signal Space Projection
  2. Identify the artifacts
  3. Calculate the projectors
  4. Evaluation
  5. Additional discussions on the forum

From Auditory

Saccades

From Continuous

Signal Space Projection

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:

Further, if we know that bn(t) is well characterized by a few field patterns b1...bm, we can express the disturbance as

where the columns of U constitute an orthonormal basis for b1...bm, cn(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., Ce = E{eeT} = 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 b1...bm such that the conditions described above are satisfied. We can now construct the orthogonal complement operator

and apply it to b(t) yielding

since P⊥bn(t) = P⊥Ucn(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:

  1. The basis set b1...bm should be able to characterize the disturbance field patterns completely and

  2. The angles between the noise subspace space spanned by b1...bm and the signal vectors bs(t) should be as close to π/2 as possible.

If the first requirement is not satisfied, some noise will leak through because P⊥bn(t) ≠ 0. If the any of the brain signal vectors bs(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.

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

Identify the artifacts

The first step is to identify several repetitions of the artifact (the vectors b1...bm). 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

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.

[ATTACH]

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:

[ATTACH] [ATTACH]

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.

[ATTACH]

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.

[ATTACH] [ATTACH]

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.

[ATTACH]

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:

[ATTACH]

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:

  1. 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.
  2. Filter the MEG or EEG signals in a frequency band of interest, in which the artifact is the most visible.

  3. Concatenate all those time blocks into a big matrix A = [b1, ..., bm]

  4. Compute the singular value decomposition of this matrix A: [U,S,V] = svd(A, 'econ')

  5. 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.

  6. Calculate the projection operator: P⊥i = I - UiUiT

  7. Apply this projection on the MEG or EEG recordings F: F = P⊥iF

  8. 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 Ui, are saved in the database link ("Link to raw file") and in the channel file.

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.

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).

  2. Calculate the SSP to remove the cardiac artifact.

  3. 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:

[ATTACH] [ATTACH]

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:

[ATTACH] [ATTACH]

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 Ui (singular vector of the SVD decomposition), and the percentage is the singular value for this vector, normalized for this decomposition (percentage = Si / sum(Si)). 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.

[ATTACH]

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.

[ATTACH]

[ATTACH]

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.

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.

[ATTACH] [ATTACH]

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.

[ATTACH]

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.

[ATTACH]

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.

[ATTACH]

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:

[ATTACH]

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:

Evaluation

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:

[ATTACH]

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.

[ATTACH]

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.

[ATTACH]

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):

[ATTACH]

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

[ATTACH]

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

[ATTACH]

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







Feedback: Comments, bug reports, suggestions, questions
Email address (if you expect an answer):


Tutorials/ArtifactsSsp (last edited 2015-02-03 23:47:06 by FrancoisTadel)