Tutorial 10: Power spectrum and frequency filters

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

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

We can correct for these artifacts using frequency filters. Usually we prefer to run these notch and band-pass filters before any other type of correction, on the continuous files. They can be applied to the recordings without much supervision, but they may create important artifacts at the beginning and the end of the signals. Processing the entire continuous recordings at once instead of the imported epochs avoids adding these edge effects to all the trials.

Evaluation of the noise level

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

• Clear the list of files in the Process1 tab.

• Select the three datasets we have linked to our protocol.
  You can select the three "link to raw file" nodes, the three folders or the entire subject node.
  
  

• Click on [Run] to open the pipeline editor window.

• Select the process "Frequency > Power spectrum density (Welch)"

• This process evaluates the power of the MEG/EEG signals at different frequencies, using the Welch's method (see Wikipedia or MathWorks). It splits the signals in overlapping windows of a given length, calculates the Fourier Transform (FFT) of each of these short segments, and averages the power of the FFT coefficients for all the overlapping windows.
  
  

• Set the options as follows (click on [Edit] for the additional options):

  • Time window: [All file]
    Portion of the input file you want to use for estimating the spectrum.
    It is common to observe huge artifacts at the beginning or the end of the recordings, in this case you should exclude these segments from the calculation of the PSD.
    In practice, using just the first 100s or 200s of the file can give you a good enough impression of the quality of the recordings.

  • Window length: 4 seconds
    Estimator length = length of the overlapping time windows for which we calculate the FFT. The number of time samples in the estimator is the same as the number of frequency bins in the output file. Increasing this parameter increases the output frequency resolution (distance between two frequency bins) but degrades the stability of the estimator, as it also decreases the total number of averaged time windows.
    A hamming window is applied to each estimator window before the computation of the FFT.

  • Overlap: 50%
    How much overlap do we want between two consecutive time windows.

  • Sensor types or names: MEG
    Defines the list of channels (names or types) on which you want to apply the process.

  • Frequency definition: Matlab's FFT default
    You have the option to directly use the frequency binning returned by the FFT, or run an additional step of averaging these bins in larger frequency bands. Note that you can freely edit these frequency bands.

  • Output: Save individual PSD value.
    This option will separately estimate the PSD for each of the three files in input, and create three files in output. If you select the other option (save average), it calculates the same three files but averages them on the fly and saves only one file in the database.

  • Implementation details: See function brainstorm3/toolbox/timefreq/bst_psd.m

• Click on [Run] to start the execution of the process.

• Troubleshooting: If you get "Out of memory" errors, try to run this PSD estimation on a shorter time segment. For example, set the time window to [0,100s] instead of the full file. This process starts by loading all the needed recordings in memory, you might not have enough memory available on your system to fit the entire dataset.

• It produces three new files, that appear as "depending" on the three datasets. The comments of the files indicate how many overlapping windows could be used to estimate the PSD. "179/4000ms" means 179 windows of 4s each (total 716s). With the 50% overlap, it sums up to a little less than 2x the file length (360s).
  
  

Interpretation of the PSD

File: AEF#01

• Double-click on the new PSD file for run #01 to display it (or right-click > Power spectrum).
  
  

• The power spectrum is displayed in a figure similar to the time series, but the X axis represents the frequencies. Most of the shortcuts described for the recordings are also valid here. Clicking on the white parts of the figure or using the arrow keys of the keyboard move the frequency cursor. The current frequency can also be controlled with a new slider displayed in the Brainstorm window, just below the time panel.

• Each black line represents the power spectrum for one channel. If you click on a channel, it gets highlighted in red. Click again for deselecting. Right-click on a channel to read its name.
  
  

• The frequency resolution of this power spectrum, ie. the distance between two frequency bins represented in this figure, is about 0.15Hz. This precision depends on the length of the estimator window you used. The FFT is computed on the number of time samples per window (4s*600Hz), rounded to the next power of 2 (nextpow2), and represents the full spectrum of the file (0-600Hz).
  Frequency resolution = sampling_freq / 2^nextpow2(estimator_length*sampling_freq) = 0.1465 Hz

• The shape of this graph is normal, it does not indicate anything unexpected:

  • Peaks related with the subject's alpha rhythms: around 10Hz and 20Hz.

  • Peaks related with the power lines: 60Hz, 120Hz and 180Hz.
    These datasets were recorded in Canada, where the alternating powerline current is delivered at 60Hz. In Europe you would observe similar peaks at 50Hz, 100Hz and 150Hz.

• Add a topography view for the same file, with one of the two methods below:

  • Right-click on the PSD file > 2D Sensor cap.

  • Right-click on the spectrum figure > 2D Sensor cap (shortcut: Ctrl+T)

• Scroll in frequencies to see the spatial distribution of each frequency bin:
  
  

• We have already identified two artifacts we will need to remove: the eye movements and the 60Hz+harmonics from the power lines.

File: AEF#02

• Open the spectrum view for the run AEF#02.
  
  

• Add a 2D sensor cap view for the same file. Scroll to zoom in/out.
  To display the sensor names, right-click on the figure > Channels > Display sensors.

• This spectrum looks very similar to the run #01: same alpha and power line peaks.

• Additionally, we observe higher signal power between 30Hz and 100Hz on many occipital sensors. This is probably related to some tension in the neck muscles due to an uncomfortable seating position in the MEG. We will see later whether these channels need to be tagged as bad or not.
  
  

File: Noise recordings

• Open the spectrum view for the noise recordings.
  
  

• This shows the power spectrum of the signals that are recorded when there is no subject in the MEG room. It gives a good and simple representation of the instrumental noise. If you had one bad MEG sensor, you would see it immediately in this graph. Here everything looks good.

X Log-scale

• One option is worth mentioning when displaying power spectra: the logarithmic scale for the X axis. Right-click on a PSD figure > Figure > X scale: Log. It is sometimes better adapted to represent this type of data than a linear scale (especially with higher sampling frequencies).
  
  

Elekta-Neuromag and EEG users

The Elekta-Neuromag MEG systems combine different types of sensors with very different amplitude ranges, therefore you would not observe the same types of figures. Same thing for EEG users, this might not look like what you observe on your recordings.

For now, keep on following these tutorials with the example dataset to learn how to use all the Brainstorm basic features. Once you're done, read additional tutorials in the section "Other analysis scenarios" to learn about the specificities related with your own acquisition system.

Apply a notch filter

For illustration purposes, we will now run a frequency filter to remove the 60Hz+harmonics from the continuous files. Notch filters are adapted for removing well identified contaminations from systems oscillating at very stable frequencies.

• Keep all the three datasets selected in the Process1 box.
  Remember to always apply the same filters on the subject recordings and the noise recordings.

• Click on [Run] to open the Pipeline editor.

• Run the process: Pre-process > Notch filter
  
  

  • Process the entire file at once: NO

  • Select the frequencies: 60, 120, 180 Hz

  • Sensor types or names: MEG

  • The higher harmonics (240Hz) are not clearly visible, and too high to bother us in this analysis.

• This process creates three new datasets, with additional "notch" tags. These output files are saved directly in the Brainstorm database, in a binary format specific to Brainstorm (.bst).
  
  

• If you delete the folders corresponding to the original files (before the filter), your original recordings in the .ds folders are not altered. If you delete the folders corresponding to the filtered files, you will lose your filtered recordings in .bst format.

• To check where the file corresponding to a "Link to raw file" is actually saved on the hard drive, right-click on it > File > Show in file explorer.

• Important: This is an optional processing step. Whether you need this on your own recordings depends on the analysis you are planning to run on the recordings (see advanced sections below).

Evaluation of the filter

• Right-click on the Process1 list > Clear list.

• Drag and drop the three filtered files in Process1.
  
  

• Run again the PSD process "Frequency > Power spectrum density (Welch)" on the new files, with the same parameters as before, to evaluate the quality of the correction.

• Double-click on the new PSD files to open them.
  
  

• Scroll to zoom in and observe what is happening around 60Hz (before / after).
  
  

• See below an example of how this filter can affect the time series: top=before, bottom=after.
  We show the reference sensor BR2 because it shows a lot more 60Hz than any MEG sensor (sensor type "MEG REF"), ie. oscillations with a period of 16.7ms.
  Note the edge effect at the beginning of the signal: the signals below are 1.5s long, the notch filter at 60Hz is visibly not performing well during the first 500ms (blue window).
  
  

Some cleaning

To avoid any confusion later, delete the links to the original files:

• Select the folders containing the original files and press Delete (or right-click > File > Delete).
  
  

• Always read the confirmation messages carefully, you will avoid some bad surprises.
  
  

• This is what your database explorer should look like at the end of this tutorial:
  
  

Note for beginners

Everything below is advanced documentation, you can skip it for now.

What filters to apply?

The frequency filters you should apply depend on the noise present in your recordings, but also on the type of analysis you are planning to use them for. This sections provides some general recommendations.

High-pass filter

• Purpose: Remove the low frequencies from the signals. Typically used for:

  • Removing the arbitrary DC offset and slow drifts of MEG sensors (< 0.2Hz),

  • Removing the artifacts occurring at low frequencies (< 1Hz, e.g. breathing or eye movements).

• Warnings:

  • Edge effects: Transient effects you should discard at the start and end of each filtered signal because the filtering window extends into time periods outside those for which you have data.
  • Avoid using on epoched data: You need long segments of recordings to run a high-pass filter.
  • Be careful with the frequency you choose if you are studying cognitive processes that may include sustained activity in some brain regions (eg. n-back memory task).

Low-pass filter

• Purpose: Remove the high frequencies from the signals. Typically used for:

  • If the components of interest are below for example 40Hz, you may discard the faster components in the signal by applying a low-pass filter with a frequency cutoff below 40Hz.
  • Removing strong noise occurring at high frequencies (eg. muscle contractions, stimulators).
  • Display averages: In an event-related design, you will import and average multiple trials. You may low-pass filter these averages for display and interpretation purposes.
  • Statistics: In an event-related study with multiple subjects, the latency of the brain response of interest may vary between subjects. Smoothing the subject-level averages before computing a statistic across subjects may help reveal the effect of interest.

• Warnings:

  • Edge effects: Transient effects you should discard at the start and end of each filtered signal because the filtering window extends into time periods outside those for which you have data.
  • It is always better to filter continuous (non-epoched data) when possible.
  • When filtering averages: Import longer epochs, average them, filter, then remove the beginning and the end of the average to keep only the signals that could be filtered properly

Band-pass filter

• Purpose: A band-pass filter is the combination of a low-pass filter and a high-pass filter, it removes all the frequencies outside of the frequency band of interest.

• Warnings: The same considerations and warnings as for high and low pass filtering apply here.

Notch filter

• Purpose: Remove a sinusoidal signal at a specific frequency (power lines noise, head tracking coils).

• Use only if needed:
  • It is not always recommended to remove the 50-60Hz power lines peaks. If you don't have a clear reason to think that these frequencies will cause a problem in your analysis, you don't need to filter them out.
  • In an ERP analysis, the averaging of multiple trials will get rid of the 50-60Hz power line oscillations because they are not time-locked to the stimulus.
  • If you are using a low-pass filter, do not a apply a notch filter at a higher frequency (useless).

• Alternatives: If the notch filter is not giving satisfying result, you have two other options.

  • Band-stop filter: Similar to the notch filter, but more aggressive on the data.
    Useful for removing larger segments of the spectrum, in case the power line peaks are spread over numerous frequency bins or for suppressing other types of artifacts.

  • Sinusoid removal: This process can do a better job at removing precise frequencies by identifying the sinusoidal components and then subtracting them from the signals in the time domain. This is not a frequency filter and works best on short segments of recordings.
    Run it on the imported epochs rather than on the continuous files.

When to apply these filters?

• Continuous files: Frequency filters used for pre-processing purposes should be applied before epoching the recordings. In general, filters will introduce transient effects at the beginning and the end of each time series, which make these parts of the data unreliable and they should be discarded. If possible, it is safer and more efficient to filter the entire recordings from the original continuous file at once.

• Before SSP/ICA cleaning: Artifact cleaning with SSP/ICA projectors require all the channels of data to be loaded in memory at the same time. Applying a frequency filter on a file that contains projectors requires all the file to be loaded and processed at once, which may cause memory issues. Pre-processing filters are rarely changed, whereas you may want to redo the SSP/ICA cleaning multiple times. Therefore it is more convenient to apply the filters first.

• Imported epochs: Filtering epochs after importing them in the database is possible but requires extra attention: you may need to import longer epochs to be able to deal with the edge effects.

• After averaging: You may low-pass filter the averaged epochs for display or statistical purposes but again be aware of the edge effects.

• Empty room measurements: In principle, all the filters that are applied to the experimental data also need to be applied, with the same settings, to the noise recordings. In the source estimation process, we will need all the files to have similar levels of noise, especially for the calculation of the noise covariance matrix. This applies in particular when some channels are noisy.

• Think first: Never apply a frequency filter without a clear reason (artifacts, predefined frequency ranges of interest, etc.) and without keeping the side effects under control. Avoid when possible.

Filter specifications: Band-pass

• Process: Pre-process > Band-pass filter

Warning: mirroring just masks the edge effects, it doesn't improve the quality of the filter.

There are two options in this process. The first one relaxes the stopband attenuation from -60db (default value) to -40db. This results to a lower order filter with a smaller edge effect, faster filtering, but with a lower accuracy.

Another option is the filtering method. After building a FIR band-pass filter, we can perfrom filtering in freuqency domain(using FFTs of input signal and filter impulse response) or in time domain (by convultion). The first approach is much faster for most cases and filters, while the later might be more practical for few cases. Anyway, both methods have a same result.

Relax: "this option reduces the transient length, but at the expense of reduced performance in the filter in terms of stopband attenuation (40dB rather than 60dB)"

The HPF is implemented as a zero-phase (zero delay) FIR filter and based on Kaiser window design.

The edge effect region lasts for a number of samples equal to half of the filter order. If the edge effect affects too much of your data, adjust the filter parameters to reduce filter order

. BrainStorm will generate a warning if your choice of filter parameters results in an edge effect of two seconds or more.BrainStorm will generate a warning if the combined edge effects at the start and end of your data represents 10 percent or more of the total number of samples.

• Description: these filters are implemented using Matlab’s FIR1 function. This uses a Kaiser window to design an even-order linear phase FIR filter, order N. Because the filters are linear phase, we can (and do) compensate for the filter delay by shifting the sequence backward in time by M=N/2 samples. This effectively makes the filters zero-phase and zero-delay.

• The user selects passband edges and indicates filter type (LPF, HPF or BPF) through the function arguments. The allowed ripple in pass and attenuation in stop band are set by default to 10^(-3) and 60dB respectively (note that with Kaiser window design, errors in pass and stopband will always be equal). Transitions between pass and stopbands are set to 15 percent of the upper and lower passband edges. However, when the lower edge of the passband is 5hz or lower we set the transition width to 50 percent of the lower passband edge. Using all of these specs we use the Kaiserord function to determine the minimum filter order required, and then design the filter (see warnings below).

• Reference: Kaiser window design is covered in almost all introductory signal processing textbooks. See also mathworks documentation for Kaiserord and FIR1 for examples and further references.

• Edge effects: With any filtering operation there will always be a transient effect at the beginning of the filtered data. For our filter, this effect will last for M=N/2 samples. We strongly recommend that your data records are sufficiently long that you can discard these M=N/2 samples. Because we are using zero-phase filtering, there is a similar N/2 effect at the end of the sampled data – these samples should also be discarded. We have included an option that allows you to mirror the data at the start and end of the record. This will reduce the apparent N/2 transients at the start and end of your data record, but you should be aware that these samples are still unreliable and we do not recommend using them. When you filter and display your data, the duration of this edge effect will be indicated. If you are losing too much of your data in these edge regions consider modifying your processing as discussed next.

• Guidance and Warnings: The key issue to be aware of when filtering is that the specification you choose for your filter will determine the length of the impulse response (or the filter order) which in turn will affect the fraction of your data that fall into the ‘edge’ region. The most likely factor that will contribute to a high filter order and large edge effects is if you specify a very low frequency at the lower edge of the passband (i.e. the high pass cut-off frequency). If your goal is to remove the DC signal we recommend you first try detrending the data (removes average and best linear fit) to see if this is sufficient. If you still need to remove other low frequency components, then pick the highest cut-off frequency that will fit your needs. When you run the filtering, BrainStorm will automatically generate a warning if the edge effect is longer than 2 seconds. Whether or not this warning is generated, you can also optionally plot the frequency and impulse response of the filter – this plot will also tell you the specifications of the filter and its order. If the edge effects are too long, you can try using the relaxed specification option – this reduces the stopband attenuation from 60dB to 40dB. If that still does produce an acceptable result in terms of edge effect, your will need to increase the lower cut-off frequency.

• If you are performing bandpass filtering and are not satisfied with the results, you can investigate filtering your data twice, once with a low pass filter and once with a high pass filter. The advantage of this is that you can now separately control the transition widths and stop band attenuation of the two filters. When designing a single BPF using the Kaiser (or any other) window, the maximum deviation from the desired response will be equal in all bands, and the transition widths will be equal at the lower and upper edges of the pass band. By instead using a LPF and a HPF you can optimize each of these processes separately using our filtering function.

• Function: brainstorm3/toolbox/math/bst_bandpass_filter.m

• External call: x = process_bandpass('Compute', X, Fs, HighPass, LowPass, Method, isMirror, isRelax)

• Code:

Filter specifications: Notch

• Description: IIR notch filter with zero-phase lag (implemented with filtfilt), 4th order.

• Reference: MatlabCentral #292960

• Edge effects: Length of startup effect ? [TODO]

• Function: brainstorm3/toolbox/process/functions/process_notch.m

• External call: x = process_notch('Compute', x, sfreq, FreqList)

• Code:
  

      % Remove the mean of the data before filtering
      xmean = mean(x,2);
      x = bst_bsxfun(@minus, x, xmean);
  
      % Define a default width
      FreqWidth = 1;
      % Define coefficients of an IIR notch filter
      delta = FreqWidth/2;
  
      % Pole radius
      r = 1 - (delta * pi / sfreq);
      theta = 2 * pi * FreqList(i) / sfreq;
      % Gain factor
      B0 = abs(1 - 2*r*cos(theta) + r^2) / (2*abs(1-cos(theta)));
      % Numerator coefficients
      B = B0 * [1, -2*cos(theta), 1];
      % Denominator coefficients
      A = [1, -2*r*cos(theta), r^2];
      % Filter signal
      x = filtfilt(B,A,x')';
  
      % Restore the mean of the signal
      x = bst_bsxfun(@plus, x, xmean);

Filter specifications: Band-stop

• Description: Butterworth IIR filter with zero-phase lag (implemented with filtfilt), 4th order.

• Reference: FieldTrip: x = ft_preproc_bandstopfilter(x, sfreq, FreqBand, [], 'but', 'twopass')

• Edge effects: ? [TODO]

• Function: brainstorm3/toolbox/process/functions/process_bandstop.m

• External call: x = process_bandstop('Compute', x, sfreq, FreqList, FreqWidth)

• Code:
  

      % Remove the mean of the data before filtering
      xmean = mean(x,2);
      x = bst_bsxfun(@minus, x, xmean);
  
      % Nyqist frequency
      Fnyq = sfreq/2;
      % Frequency band to remove (default value FreqWidth=1.5)
      FreqBand = [FreqList(i) - FreqWidth/2, FreqList(i) + FreqWidth/2];
      % Filter order
      N = 4;
  
      % Butterworth filter
      [B,A] = butter(N, FreqBand ./ Fnyq, 'stop');
      % Filter signal
      x = filtfilt(B, A, x')';
  
      % Restore the mean of the signal
      x = bst_bsxfun(@plus, x, xmean);

On the hard drive

The names of the files generated by the process "Power spectrum density" start with the tag timefreq_psd, they share the same structure as all the files that include a frequency dimension.

To explore the contents of a PSD file created in this tutorial, right-click on it and use the popup menus
File > View file contents or File > Export to Matlab.

• 

Structure of the time-frequency files: timefreq_psd_*.mat

• TF: [Nsignals x Ntime x Nfreq] Stores the spectrum information. Nsignals is the number of channels that were selected with the option "MEG" in the PSD process. Nfreq is the number of frequency bins. There is no time dimension (Ntime = 1).

• Comment: String used to represent the file in the database explorer.

• Time: Window of time over which the file was estimated.

• TimeBands: Defined only when you select the option "Group in time bands".
  Always empty for the PSD files because there is no time dimension.

• Freqs: [1 x Nfreq] List of frequencies for which the power spectrum was estimated (in Hz).

• RefRowNames: Only used for connectivity results.

• RowNames: [Nsignals x 1] Describes the rows of the TF matrix (first dimension). Here it corresponds to the name of the MEG sensors, in the same order as is the .TF field.

• Measure: Function currently applied to the FFT coefficients {power, none, magnitude, log, other}

• Method: Function that was used to produce this file {psd, hilbert, morlet, corr, cohere, ...}

• DataFile: File from which this file was calculated = Parent file in the database explorer.

• DataType: Type of file from which this file was calculated (file type of .DataFile).

  • 'data' = Recordings
  • 'cluster' = Recordings grouped by clusters of sensors
  • 'results' = Source activations
  • 'scouts' = Time series associated with a region of interest

• SurfaceFile / GridLoc / GridAtlas / Atlas: Used only when the input file was a source file.

• nAvg: Number of input files averaged to produce this file.

• Options: Most relevant options that were passed to the function bst_timefreq.

• History: Describes all the operations that were performed with Brainstorm on this file. To get a better view of this piece of information, use the menu File > View file history.

Useful functions

• in_bst_timefreq(PsdFile): Read a PSD or time-frequency file.

• in_bst(FileName): Read any Brainstorm file.

• bst_process('LoadInputFile', FileName, Target): The most high-level function for reading Brainstorm files. "Target" is a string with the list of sensor names or types to load (field RowNames).

• bst_psd(F, sfreq, WinLength, WinOverlap): Computation of the Welch's power spectrum density.

---