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

<<TableOfContents(2,2)>>

== 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. <<BR>>You can select the three "link to raw file" nodes, the three folders or the entire subject node. <<BR>><<BR>> {{attachment:psd_select_files.gif||width="400",height="121"}}

 * Click on [Run] to open the pipeline editor window.
 * Select the process "'''Frequency > Power spectrum density (Welch)'''"

 * This process evaluates the power of the MEG/EEG signals at different frequencies, using the Welch's method (see [[http://en.wikipedia.org/wiki/Welch's_method|Wikipedia]] or [[http://www.mathworks.com/help/signal/ref/pwelch.html|MathWorks]]). It splits the signals in overlapping windows of a given length, calculates the Fourier Transform (FFT) of each of these short segments, and averages the power of the FFT coefficients for all the overlapping windows. <<BR>><<BR>> {{attachment:psd_options.png||height="388"}}

 * Set the options as follows (click on [Edit] for the additional options):
  * '''Time window''': [All file] <<BR>>Portion of the input file you want to use for estimating the spectrum. <<BR>>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. <<BR>>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 <<BR>>Estimator length = length of the overlapping time windows for which we calculate the FFT. The number of time samples in the estimator is the same as the number of frequency bins in the output file. Increasing this parameter increases the output frequency resolution (distance between two frequency bins) but degrades the stability of the estimator, 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. See forum post: [[https://neuroimage.usc.edu/forums/t/inconsistent-psd-depending-on-window-length/24957|Effect of window length on the PSD]]

  * '''Overlap''': 50% <<BR>>How much overlap do we want between two consecutive time windows.

  * '''Units''': Physical: U^2/Hz <<BR>>Scaling of the spectrum. This only affects the values on the Y axis of the spectrum. Physical units should be used in most cases. <<BR>>"Normalized" gives normalized frequencies from 0 to 2pi (Hz·s). <<BR>>"Before Nov 2020" reproduces the older Brainstorm spectrum scaling (see [[https://neuroimage.usc.edu/forums/t/psd-amplitude-factor/22258|this forum post]]).

  * '''Sensor types or names''': MEG <<BR>>Defines the list of channels (names or types) on which you want to apply the process.

  * '''Frequency definition''': Matlab's FFT default <<BR>>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. <<BR>>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). <<BR>><<BR>> {{attachment:psd_output_files.gif||width="235",height="240"}}

== Interpretation of the PSD ==
==== File: AEF#01 ====
 * Double-click on the new PSD file for run #01 to display it (or right-click > Power spectrum). <<BR>><<BR>> {{attachment:psd_display.gif}}

 * The power spectrum is displayed in a figure similar to the time series, but the X axis represents the frequencies. Most of the shortcuts described for the recordings are also valid here. Clicking on the white parts of the figure or using the arrow keys of the keyboard move the frequency cursor. The current frequency can also be controlled with a new slider displayed in the Brainstorm window, just below the time panel.
 * Each black line 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. <<BR>><<BR>> {{attachment:psd_run1.png||height="194"}}

 * 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).<<BR>>'''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'''.<<BR>>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: <<BR>><<BR>> {{attachment:psd_topo_ncm4.png||width="604",height="149"}}

 * 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. <<BR>>To view the signal units instead of dB, select '''Display Tab > Measure > Power'''.  Then from the display options icon on the right of the figure, select '''Amplitude > Log scale'''<<BR>><<BR>> {{attachment:psd_run2_units.png||height="200"}} <<BR>><<BR>> {{attachment:psd_run2.png||height="154"}}

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

 * This spectrum looks very similar to the run #01: same alpha and power line peaks.
 * Additionally, we observe higher 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. <<BR>><<BR>> {{attachment:psd_neck.png||height="198"}}

==== File: Noise recordings ====
 * Open the spectrum view for the noise recordings. <<BR>><<BR>> {{attachment:psd_noise.png||height="139"}}

 * 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. In the diaply options for the PSD figure > '''Frequency > Log scale'''. It is sometimes better adapted to represent this type of data than a linear scale (especially with higher sampling frequencies). <<BR>><<BR>> {{attachment:psd_log.png||height="220"}}

== 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 "[[http://neuroimage.usc.edu/brainstorm/Tutorials#Other_analysis_scenarios|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.<<BR>>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''' <<BR>><<BR>> {{attachment:notch_process.gif||width="709",height="313"}}
  * Process the entire file at once: '''NO'''

  * Sensor types or names: '''MEG'''

  * Frequencies to remove: '''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). <<BR>><<BR>> {{attachment:notch_output.gif||width="248",height="235"}}

 * 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. <<BR>><<BR>> {{attachment:notch_psd.gif||width="361",height="134"}}

 * 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. <<BR>><<BR>> {{attachment:notch_evaluation.png||height="420"}}

 * Scroll to zoom in and observe what is happening around 60Hz (before / after). <<BR>><<BR>> {{attachment:notch_zoom.png||height="151"}}

 * See below an example of how this filter can affect the time series: top=before, bottom=after.<<BR>> 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. <<BR>>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). <<BR>><<BR>> {{attachment:notch_timeseries2.gif||width="664",height="274"}}

 * If you look in the list of events, you would see a new category "'''transient_notch'''". This corresponds to the time period during which can expect significant edge effects due to the filtering. Brainstorm doesn't mark these blocks as bad by default, you would have to do it manually - you will see how to do this in one of the following tutorials. In the case of this dataset, the transient duration is much shorter than the delay before the first stimulation, so not relevant in our processing pipeline. See the advanced sections below for more details about the estimation of this transient duration. <<BR>><<BR>>{{attachment:notch_transient.gif}}

== 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''').<<BR>><<BR>> {{attachment:notch_delete.gif||width="405",height="194"}}

 * Always read the confirmation messages carefully, you will avoid some bad surprises. <<BR>><<BR>> {{attachment:notch_confirmation.gif||width="479",height="152"}}

 * This is what your database explorer should look like at the end of this tutorial: <<BR>><<BR>> {{attachment:notch_final.gif}}

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

<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?skip=Tutorials/BadChannels")>>

<<TAG(Advanced)>>

== 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''':<<BR>>
  * 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.<<BR>>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.<<BR>>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<<BR>><<BR>> {{attachment:specs_bandpass.gif||width="670",height="338"}}

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

  * '''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 [[https://www.mathworks.com/help/signal/ref/kaiser.html|Kaiser window]] design. The order N is estimated using Matlab's [[https://www.mathworks.com/help/signal/ref/kaiserord.html|kaiserord]] function and the filter generated with [[https://www.mathworks.com/help/signal/ref/fir1.html|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 ([[https://www.mathworks.com/help/signal/ref/fftfilt.html|fftfilt]] function) or in time domain ([[https://www.mathworks.com/help/matlab/ref/filter.html|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. <<BR>><<BR>> {{attachment:bandpass_transient.gif||width="492",height="147"}}

  * '''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.<<BR>><<BR>> {{attachment:bandpass_chirp.gif||width="497",height="344"}}

 * '''Function: '''brainstorm3/toolbox/math/bst_bandpass_hfilter.m

 * '''External call:''' process_bandpass('Compute', x, Fs, HighPass, LowPass, 'bst-hfilter', isMirror, isRelax)

 * '''Code''': <<BR>><<HTML(<div    style="border:1px solid black;  background-color:#EEEEFF;  width:720px;   height:250px; overflow:scroll;  padding:10px;  font-family:   Consolas,Menlo,Monaco,Lucida  Console,Liberation  Mono,DejaVu Sans   Mono,Bitstream Vera Sans  Mono,Courier  New,monospace,sans-serif;   font-size: 13px; white-space:  pre;">)>><<EmbedContent("http://neuroimage.usc.edu/bst/viewcode.php?math=bst_bandpass_hfilter.m")>><<HTML(</div >)>>

== Filter specifications: Notch ==
 * '''Description''': 2nd order IIR notch filter with zero-phase lag (implemented with [[http://www.mathworks.com/help/signal/ref/filtfilt.html|filtfilt]]).

 * '''Reference''': Mitra, Sanjit Kumar, and Yonghong Kuo. Digital signal processing: a computer-based approach. Vol. 2. New York: McGraw-Hill, 2006. [[http://www.mathworks.com/matlabcentral/newsreader/view_thread/292960|MatlabCentral #292960]]

 * '''Edge effects''': It is computed based on the 99% energy of the estimated impulse response.

 * '''Function''': [[https://github.com/brainstorm-tools/brainstorm3/blob/master/toolbox/process/functions/process_notch.m#L142|brainstorm3/toolbox/process/functions/process_notch.m]]

 * '''External call''': [x, FiltSpec, Messages] = Compute(x, sfreq, FreqList, Method, bandWidth)

== Filter specifications: Band-stop ==
 * '''Description''': 4th order Butterworth IIR filter with zero-phase lag (implemented with [[http://www.mathworks.com/help/signal/ref/filtfilt.html|filtfilt]])

 * '''Reference''': FieldTrip: x = [[http://www.fieldtriptoolbox.org/reference/ft_preproc_bandstopfilter|ft_preproc_bandstopfilter]](x, sfreq, FreqBand, [], 'but', 'twopass')

 * '''Edge effects''' It is computed based on the 99% energy of the estimated impulse response.

 * '''Function''': brainstorm3/toolbox/process/functions/process_bandstop.m

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

== 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 <<BR>>'''File > View file contents''' or '''File > Export to Matlab'''.

 . {{attachment:psd_contents.gif||width="579",height="420"}}

==== 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". <<BR>>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.

 * '''Leff''': Effective number of averages = 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.

== Additional documentation ==
'''Forum discussions'''

 * [[https://neuroimage.usc.edu/forums/t/inconsistent-psd-depending-on-window-length/24957|Effect of window length on the PSD]]

<<HTML(<!-- END-PAGE -->)>>

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

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