# Tutorial 10: Power spectrum and frequency filters

*Authors: Hossein Shahabi, 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.

Contents

- Evaluation of the noise level
- Interpretation of the PSD
- Elekta-Neuromag and EEG users
- Apply a notch filter
- Evaluation of the filter
- Some cleaning
- Note for beginners
- What filters to apply?
- When to apply these filters?
- Filter specifications: Low-pass, high-pass, band-pass
- Filter specifications: Notch [TODO]
- Filter specifications: Band-stop [TODO]
- On the hard drive

## 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 selected 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

This filter has been updated in 2019. In the new configuration, the user can define the 3-dB bandwidth of the filter. Please consider that a smaller bandwidth means a sharper filter which in some cases makes filter unstable. In case you want to reproduce the old filter, you can check the box "Use old filter implementation". The 3-dB bandwidth is not applicable to the old configuration.

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**

Sensor types or names:

**MEG**Select the frequencies:

**60, 120, 180 Hz**3-dB notch bandwidth:

**2 Hz**- 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).**Warnings**:

- 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: Low-pass, high-pass, band-pass

**Process**: Pre-process > Band-pass filter

**Process options**:**Lower cutoff frequency**: Defines a high-pass filter (enter 0 for a low-pass filter)**Upper cutoff frequency**: Defines a low-pass filter (enter 0 for a high-pass filter)**Stopband attenuation**: The higher the attenuation, the higher the performance of the filter, but longer the transient duration. Use the default (60dB) unless you need shorter edge effects.**Mirror signal**: Replicates the signal with a mirror symmetry on each side, to minimize the amplitude of the edge effects that these filters can generate. Otherwise the signal is padded with zeros on each side before filtering. We do not recommend using this option, because mirroring just masks the transient effects, but does not improve the quality of the filter at the beginning and the end of the signals. You need to discard the same duration of signal whether you are mirroring or zero-padding.**Use old filter**: For replicating results obtained with older versions of Brainstorm.**View filter response**: Click on this button to display the impulse and frequency response of your filter, and confirm that the responses appear reasonable.

**Filter design**:**Description**: Even-order linear phase FIR filter, based on a Kaiser window design. The order N is estimated using Matlab's kaiserord function and the filter generated with fir1. 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.**Ripple and attenuation**: 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.**Filtering function**: The FIR bandpass filter can be performed in frequency domain (fftfilt function) or in time domain (filter function). The two approaches give the same results, but they have different execution times depending on the filter oder. The time-domain filtering is faster for low-order filters and much slower for high-order filters. The process selects automatically what approach to use.

**Edge effects:****Transient (full)**: 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 half of the filter order: 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.**Transient (99% energy)**: For some filters, the full transient window might be longer than your epochs. However, most of the energy is carried by the beginning of the filter, and you can obtain amplitudes acceptable for most analysis after a fraction of this full window. For this reason we also mention a much shorter window in the documentation of the filter, which corresponds to the duration needed to obtain 99% of the total energy in the impulse response. This duration corresponds to the "transient" event markers that are added to the recordings when applying filters.

**Adjust the parameters**: If possible, always discard the full transient window. If the edge effect affects too much of your data, adjust the filter parameters to reduce filter order (increase the lower cut-off frequency or reduce the stopband attenuation). If you cannot get any acceptable compromise, you can consider discarding shorter transient windows, but never go below this "99% energy" window.**Mirroring**: We included an option to mirror the data at the start and end of the record instead of padding the signal with zeros. 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.**[TODO]**Check this "99% energy" criteria in the case of high-pass filters, it does not seem very useful...

**Additional recommendations:****Filter order**: 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).**Detrending**: 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.**Design optimization**: 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.

**Linear phase, no distortion, zero delay:**As described earlier, FIR filters have a linear phase in the frequency domain. It means that all samples of the input signal will have a same delay in the output. This delay is compensated after filtering. Consequently, no distortion happens during the filtering process. To illustrate this property, we considered a chirp signal in which the oscillation frequency grows linearly. The signal is band-pass filtered in two frequency ranges. The following plot represents the original signal and its filtered versions with our proposed filters. Results show that the input and output signals of this filter are completely aligned without any delay or distortion.

**Function:**brainstorm3/toolbox/math/bst_bandpass_hfilter.m**External call:**process_bandpass('Compute', x, Fs, HighPass, LowPass, 'bst-hfilter', isMirror, isRelax)**Code**:

1 function [x, FiltSpec, Messages] = bst_bandpass_hfilter(x, Fs, HighPass, LowPass, isMirror, isRelax, Function, Rolloff) 2 % BST_BANDPASS_HFILTER Linear phase FIR bandpass filter. 3 % 4 % USAGE: [x, FiltSpec, Messages] = bst_bandpass_hfilter(x, Fs, HighPass, LowPass, isMirror=0, isRelax=0, Function=[detect], Rolloff=[]) 5 % [~, FiltSpec, Messages] = bst_bandpass_hfilter([], Fs, HighPass, LowPass, isMirror=0, isRelax=0, Function=[detect], Rolloff=[]) 6 % x = bst_bandpass_hfilter(x, Fs, FiltSpec) 7 % 8 % DESCRIPTION: 9 % - A linear phase FIR filter is created. 10 % - Function "kaiserord" and "kaiser" are used to set the necessary order for fir1. 11 % - The transition band is hard-coded. 12 % - Requires Signal Processing Toolbox for the following functions: 13 % kaiserord, kaiser, fir1, fftfilt. If not, using Octave-based alternatives. 14 % 15 % INPUT: 16 % - x : [nChannels,nTime] input signal (empty to only get the filter specs) 17 % - Fs : Sampling frequency 18 % - HighPass : Frequency below this value are filtered in Hz (set to 0 for low-pass filter only) 19 % - LowPass : Frequency above this value are filtered in Hz (set to 0 for high-pass filter only) 20 % - isMirror : isMirror (default = 0 no mirroring) 21 % - isRelax : Change ripple and attenuation coefficients (default=0 no relaxation) 22 % - Function : 'fftfilt', filtering in frequency domain (default) 23 % 'filter', filtering in time domain 24 % If not specified, detects automatically the fastest option based on the filter order 25 % - Rolloff : Width of the transition band in Hz 26 % 27 % OUTPUT: 28 % - x : Filtered signals 29 % - FiltSpec : Filter specifications (coefficients, length, ...) 30 % - Messages : Warning messages, if any 31 32 % @============================================================================= 33 % This function is part of the Brainstorm software: 34 % https://neuroimage.usc.edu/brainstorm 35 % 36 % Copyright (c)2000-2019 University of Southern California & McGill University 37 % This software is distributed under the terms of the GNU General Public License 38 % as published by the Free Software Foundation. Further details on the GPLv3 39 % license can be found at http://www.gnu.org/copyleft/gpl.html. 40 % 41 % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE 42 % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY 43 % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF 44 % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY 45 % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. 46 % 47 % For more information type "brainstorm license" at command prompt. 48 % =============================================================================@ 49 % 50 % Authors: Hossein Shahabi, Francois Tadel, John Mosher, Richard Leahy, 2016 51 52 53 %% ===== PARSE INPUTS ===== 54 % Filter is already computed 55 if (nargin == 3) 56 FiltSpec = HighPass; 57 % Default filter options 58 else 59 if (nargin < 8) || isempty(Rolloff) 60 Rolloff = []; 61 end 62 if (nargin < 7) || isempty(Function) 63 Function = []; % Auto-detection based on the filter order later in the code 64 end 65 if (nargin < 6) || isempty(isRelax) 66 isRelax = 0; 67 end 68 if (nargin < 5) || isempty(isMirror) 69 isMirror = 0; 70 end 71 FiltSpec = []; 72 end 73 Messages = []; 74 75 76 %% ===== CREATE FILTER ===== 77 if isempty(FiltSpec) 78 % ===== FILTER SPECIFICATIONS ===== 79 Nyquist = Fs/2; 80 % High-pass filter 81 if ~isempty(HighPass) && (HighPass ~= 0) 82 f_highpass = HighPass / Nyquist; % Change frequency from Hz to normalized scale (0-1) 83 % Default transition band 84 if isempty(Rolloff) 85 if (HighPass <= 5) % Relax the transition band if HighPass<5 Hz 86 f_highstop = .5 * f_highpass; 87 else 88 f_highstop = .85 * f_highpass; 89 end 90 % Specified manually 91 else 92 f_highstop = max(0.2, Highpass - Rolloff) / Nyquist; 93 end 94 else 95 f_highpass = 0; 96 f_highstop = 0; 97 end 98 % Low-pass filter 99 if ~isempty(LowPass) && (LowPass ~= 0) 100 f_lowpass = LowPass / Nyquist; 101 % Default transition band 102 if isempty(Rolloff) 103 if f_highpass==0 % If this is a low-pass filter 104 f_lowstop = 1.05 * f_lowpass; 105 else 106 f_lowstop = 1.15 * f_lowpass; 107 end 108 % Specified manually 109 else 110 f_lowstop = min(Fs/2 - 0.2, Lowpass + Rolloff) / Nyquist; 111 end 112 else 113 f_lowpass = 0; 114 f_lowstop = 0; 115 end 116 % If both high-pass and low-pass are zero 117 if (f_highpass == 0) && (f_lowpass == 0) 118 Messages = ['No frequency band in input.' 10]; 119 return; 120 % Input frequencies are too high 121 elseif (f_highpass >= 1) || (f_lowpass >= 1) 122 Messages = sprintf('Cannot filter above %dHz.\n', Nyquist); 123 return; 124 end 125 % Transition parameters 126 if isRelax 127 Ripple = 10^(-2); 128 Atten = 10^(-2); % Equals 40db 129 else 130 Ripple = 10^(-3); % pass band ripple 131 Atten = 10^(-3); % Equals 60db 132 end 133 134 % ===== DESIGN FILTER ===== 135 % Build the general case first 136 fcuts = [f_highstop, f_highpass, f_lowpass, f_lowstop]; 137 mags = [0 1 0]; % filter magnitudes 138 devs = [Atten Ripple Atten]; % deviations 139 % Now adjust for desired properties 140 fcuts = max(0,fcuts); % Can't go below zero 141 fcuts = min(1-eps, fcuts); % Can't go above or equal to 1 142 143 % We have implicitly created a bandpass, but now adjust for desired filter 144 if (f_lowpass == 0) % User didn't want a lowpass 145 fcuts(3:4) = []; 146 mags(3) = []; 147 devs(3) = []; 148 end 149 if (f_highpass == 0) % User didn't want a highpass 150 fcuts(1:2) = []; 151 mags(1) = []; 152 devs(1) = []; 153 end 154 155 % Generate FIR filter 156 % Using Matlab's Signal Processing toolbox 157 if bst_get('UseSigProcToolbox') 158 [n,Wn,beta,ftype] = kaiserord(fcuts, mags, devs, 2); 159 n = n + rem(n,2); % ensure even order 160 b = fir1(n, Wn, ftype, kaiser(n+1,beta), 'noscale'); 161 % Using Octave-based functions 162 else 163 [n,Wn,beta,ftype] = oc_kaiserord(fcuts, mags, devs, 2); 164 n = n + rem(n,2); % ensure even order 165 b = oc_fir1(n, Wn, ftype, oc_kaiser(n+1,beta), 'noscale'); 166 end 167 168 % Filtering function: Detect the fastest option, if not explicitely defined 169 if isempty(Function) 170 % The filter() function is a bit faster for low-order filters, but much slower for high-order filters 171 if (n > 800) % Empirical threshold 172 Function = 'fftfilt'; 173 else 174 Function = 'filter'; 175 end 176 end 177 178 % Compute the cumulative energy of the impulse response 179 E = b((n/2)+1:end) .^ 2 ; 180 E = cumsum(E) ; 181 E = E ./ max(E) ; 182 % Compute the effective transient: Number of samples necessary for having 99% of the impulse response energy 183 [tmp, iE99] = min(abs(E - 0.99)) ; 184 185 % Output structure 186 FiltSpec.b = b; 187 FiltSpec.a = 1; 188 FiltSpec.order = n; 189 FiltSpec.transient = iE99 / Fs ; % Start up and end transients in seconds (Effective) 190 % FiltSpec.transient_full = n / (2*Fs) ; % Start up and end transients in seconds (Actual) 191 FiltSpec.f_highpass = f_highpass; 192 FiltSpec.f_lowpass = f_lowpass; 193 FiltSpec.fcuts = fcuts * Nyquist ; % Stop and pass bands in Hz (instead of normalized) 194 FiltSpec.function = Function; 195 FiltSpec.mirror = isMirror; 196 % If empty input: just return the filter specs 197 if isempty(x) 198 return; 199 end 200 end 201 202 %% ===== FILTER SIGNALS ===== 203 % Transpose signal: [time,channels] 204 [nChan, nTime] = size(x); 205 % Half of filter length 206 M = FiltSpec.order / 2; 207 % If filter length > 10% of data length 208 edgePercent = 2*FiltSpec.transient / (nTime / Fs); 209 if (edgePercent > 0.1) 210 Messages = [Messages, sprintf('Start up and end transients (%.2fs) represent %.1f%% of your data.\n', 2*FiltSpec.transient, 100*edgePercent)]; 211 end 212 213 % Remove the mean of the data before filtering 214 xmean = mean(x,2); 215 x = bst_bsxfun(@minus, x, xmean); 216 217 % Mirroring requires the data to be longer than the filter 218 if (FiltSpec.mirror) && (nTime < M) 219 Messages = [Messages, 'Warning: Data is too short for mirroring. Option is ignored...' 10]; 220 FiltSpec.mirror = 0; 221 end 222 % Mirror signals 223 if (FiltSpec.mirror) 224 x = [fliplr(x(:,1:M)), x, fliplr(x(:,end-M+1:end))]; 225 % Zero-padding 226 else 227 x = [zeros(nChan,M), x, zeros(nChan,M)] ; 228 end 229 230 % Filter signals 231 switch (FiltSpec.function) 232 case 'fftfilt' 233 if bst_get('UseSigProcToolbox') 234 x = fftfilt(FiltSpec.b, x')'; 235 else 236 x = oc_fftfilt(FiltSpec.b, x')'; 237 end 238 case 'filter' 239 x = filter(FiltSpec.b, FiltSpec.a, x, [], 2); 240 end 241 242 % Remove extra data 243 x = x(:,2*M+1:end); 244 % Restore the mean of the signal (only if there is no high-pass filter) 245 if (FiltSpec.f_highpass == 0) 246 x = bst_bsxfun(@plus, x, xmean); 247 end 248 249

## Filter specifications: Notch [TODO]

**Description**: 2nd order IIR notch filter with zero-phase lag (implemented with filtfilt).**Reference**:Mitra, Sanjit Kumar, and Yonghong Kuo. Digital signal processing: a computer-based approach. Vol. 2. New York: McGraw-Hill, 2006.

**Edge effects**: It is computed based on the 99% energy of the estimated impulse response.**Function**: brainstorm3/toolbox/process/functions/process_notch.m**External call**: [x, FiltSpec, Messages] = Compute(x, sfreq, FreqList, Method, bandWidth)**Code**

brainstorm3/process_notch.m at master · brainstorm-tools/brainstorm3 · GitHub Permalink### Join GitHub today

GitHub is home to over 31 million developers working together to host and review code, manage projects, and build software together.

Sign upfunction varargout = process_notch( varargin ) % PROCESS_NOTCH: Remove one or more sinusoids from a signal % % USAGE: sProcess = process_notch('GetDescription') % sInput = process_notch('Run', sProcess, sInput) % x = process_notch('Compute', x, sfreq, FreqList) % @============================================================================= % This function is part of the Brainstorm software: % https://neuroimage.usc.edu/brainstorm % % Copyright (c)2000-2019 University of Southern California & McGill University % This software is distributed under the terms of the GNU General Public License % as published by the Free Software Foundation. Further details on the GPLv3 % license can be found at http://www.gnu.org/copyleft/gpl.html. % % FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE % UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY % WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF % MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY % LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE. % % For more information type "brainstorm license" at command prompt. % =============================================================================@ % % Authors: Francois Tadel, 2014-2015 % % Code inspired from MatlabCentral post: % http://www.mathworks.com/matlabcentral/newsreader/view_thread/292960 eval(macro_method); end %% ===== GET DESCRIPTION ===== function sProcess = GetDescription() %#ok<DEFNU> % Description the process sProcess.Comment = 'Notch filter'; sProcess.FileTag = 'notch'; sProcess.Category = 'Filter'; sProcess.SubGroup = 'Pre-process'; sProcess.Index = 66; sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/ArtifactsFilter#Filter_specifications:_Notch'; % Definition of the input accepted by this process sProcess.InputTypes = {'data', 'results', 'raw', 'matrix'}; sProcess.OutputTypes = {'data', 'results', 'raw', 'matrix'}; sProcess.nInputs = 1; sProcess.nMinFiles = 1; sProcess.processDim = 1; % Process channel by channel % Definition of the options % === Freq list sProcess.options.freqlist.Comment = 'Frequencies to remove (Hz):'; sProcess.options.freqlist.Type = 'value'; sProcess.options.freqlist.Value = {[], 'list', 2}; % === Sensor types sProcess.options.sensortypes.Comment = 'Sensor types or names (empty=all): '; sProcess.options.sensortypes.Type = 'text'; sProcess.options.sensortypes.Value = 'MEG, EEG'; sProcess.options.sensortypes.InputTypes = {'data', 'raw'}; end %% ===== FORMAT COMMENT ===== function Comment = FormatComment(sProcess) %#ok<DEFNU> if isempty(sProcess.options.freqlist.Value{1}) Comment = 'Notch filter: No frequency selected'; else strValue = sprintf('%1.0fHz ', sProcess.options.freqlist.Value{1}); Comment = ['Notch filter: ' strValue(1:end-1)]; end end %% ===== RUN ===== function sInput = Run(sProcess, sInput) %#ok<DEFNU> % Get options FreqList = sProcess.options.freqlist.Value{1}; if isempty(FreqList) || isequal(FreqList, 0) bst_report('Error', sProcess, [], 'No frequency in input.'); sInput = []; return; end % Get sampling frequency sfreq = 1 ./ (sInput.TimeVector(2)-sInput.TimeVector(1)); % % Test length of the signal % if (size(sInput.A,2) < round(sfreq)) % bst_report('Warning', sProcess, [], 'Signal is too short for performing a proper filtering. Minimum duration = 1s'); % end % Filter data sInput.A = Compute(sInput.A, sfreq, FreqList); % Comment strValue = sprintf('%1.0fHz ', FreqList); sInput.CommentTag = [sProcess.FileTag '(' strValue(1:end-1) ')']; % Do not keep the Std field in the output if isfield(sInput, 'Std') && ~isempty(sInput.Std) sInput.Std = []; end end %% ===== EXTERNAL CALL ===== % USAGE: x = process_notch('Compute', x, sfreq, FreqList) function x = Compute(x, sfreq, FreqList) % Use the signal processing toolbox? UseSigProcToolbox = bst_get('UseSigProcToolbox'); % Check list of freq to remove if isempty(FreqList) || isequal(FreqList, 0) return; end % Define a default width FreqWidth = 1; % Remove the mean of the data before filtering xmean = mean(x,2); x = bst_bsxfun(@minus, x, xmean); % Remove all the frequencies sequencially for ifreq = 1:length(FreqList) % Define coefficients of an IIR notch filter delta = FreqWidth/2; % Pole radius r = 1 - (delta * pi / sfreq); theta = 2 * pi * FreqList(ifreq) / 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 if UseSigProcToolbox x = filtfilt(B,A,x')'; else x = filter(B,A,x')'; x(:,end:-1:1) = filter(B,A,x(:,end:-1:1)')'; end end % Restore the mean of the signal x = bst_bsxfun(@plus, x, xmean); end Press h to open a hovercard with more details.

## Filter specifications: Band-stop [TODO]

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