Differences between revisions 23 and 127 (spanning 104 versions)

Tutorial 24: Time-frequency

Authors: Francois Tadel, Dimitrios Pantazis, Elizabeth Bock, Sylvain Baillet

This tutorial introduces how to compute time-frequency decomposition of MEG/EEG recordings and cortical currents using complex Morlet wavelets and Hilbert transforms.

Contents

Introduction
Morlet wavelets
Edge effects
Simulation
Process options
Display: Time-frequency map
Display: Mouse and keyboard shortcuts
Display: Power spectrum and time series
Normalized time-frequency maps
Tuning the wavelet parameters
Hilbert transform [TODO]
MEG recordings: Single trials
Display: All channels
Display: Topography
Scouts
Full cortical maps
Unconstrained sources
Getting rid of the edge effects
On the hard drive
References
Additional documentation

Introduction

Some of the MEG/EEG signal properties are difficult to access in time domain (graphs time/amplitude). A lot of the information of interest is carried by oscillations at certain frequencies, but the amplitude of these oscillations is sometimes a lot lower than the amplitude of the slower components of the signal, making them difficult to observe.

The averaging in time domain may also lead to a cancellation of these oscillations when they are not strictly locked in phase across trials. Averaging trials in time-frequency domain allows to extract the power of the oscillation regardless of the phase shifts. For a better understanding of this topic, we recommend the lecture of the following article: Bertrand O, Tallon-Baudry C (2000), Oscillatory gamma activity in humans: a possible role for object representation

In Brainstorm we offer two approaches for computing time-frequency decomposition (TF): the first one is based on the convolution of the signal with series of complex Morlet wavelets, the second filters the signal in different frequency bands and extracts the envelope of the filtered signals using the Hilbert transform.

Morlet wavelets

Complex Morlet wavelets are very popular in EEG/MEG data analysis for time-frequency decomposition. They have the shape of a sinusoid, weighted by a Gaussian kernel, and they can therefore capture local oscillatory components in the time series. An example of this wavelet is shown below, where the blue and red curves represent the real and imaginary part, respectively.

Contrary to the standard short-time Fourier transform, wavelets have variable resolution in time and frequency. For low frequencies, the frequency resolution is high but the time resolution is low. For high frenqucies, it's the opposite. When designing the wavelet, we basically decide a trade-off between temporal and spectral resolution.

To design the wavelet, we first need to choose a central frequency, ie. the frequency where we will define the mother wavelet. All other wavelets will be scaled and shifted versions of the mother wavelet. Unless interested in designing the wavelet at a particular frequency band, the default 1Hz should be fine.

Then, the desirable time resolution for the central frequency should be defined. For example, we may wish to have a temporal resolution of 3 seconds at frequency 1 Hz (default parameters). These two parameters, uniquely define the temporal and spectral resolution of the wavelet for all other frequencies, as shown in the plots below. Resolution is given in units of Full Width Half Maximum of the Gaussian kernel, both in time and frequency.

Edge effects

Users should pay attention to edge effects when applying wavelet analysis. Wavelet coefficients are computed by convolving the wavelet kernel with the time series. Similarly to any convolution of signals, there is zero padding at the edges of the time series and therefore the wavelet coefficients are weaker at the beginning and end of the time series.

From the figure above, which designs the Morlet wavelet, we can see that the default wavelet (central frequency fc=1Hz, FWHM_tc=3sec) has temporal resolution 0.6sec at 5Hz and 0.3sec at 10Hz. In such case, the edge effects are roughly half these times: 300ms in 5Hz and 150ms in 10Hz.

More precisely, if f is your frequency of interest, you can expect the edge effects to span over FWHM_t seconds: FWHM_t = FWHM_tc * fc / f / 2. Examples of such transients are given in the figures below.

Simulation

We will illustrate the time-frequency decomposition process with a simulated signal.

The following code generates a sum of three sinusoids (2Hz, 20Hz, 50Hz) with random white noise. The 50Hz and noise are present everywhere, the 2Hz and 20Hz start only after two seconds.
```
f1 = 2; f2 = 20; f3 = 50;
i =2000:6000;
Data(1,i) = sin(f1*2*pi*t(i)) + 0.4 * cos(f2*2*pi*t(i));
Data = Data + 0.2 * sin(f3*2*pi*t) + 0.4 * rand(1,6000);
```
Empty the Process1 list (right-click > Clear list) then click on [Run].
Run process: Simulate > Simulate generic signals.
Ntime=6000, Sampling frequency=1000Hz (signal duration = 6000/1000 = 6 seconds).
Copy-paste the few lines of code above to generate the sum of three sinusoids.
Double-click on the new file to look at the simulated signal.
In Process1, select the simulated signal.
Run process: Frequency > Time-frequency (Morlet wavelets).
Select the option 1/f compensation: The normalization will be discussed later.
Click on the button [Edit] to see all the process options.
Time definition: Same as input, Frequency definition: Linear 1:1:60, Compute measure: Power.

Process options

Comment: String that will be displayed in the database explorer to represent the output file.

Time definition

Same as input file: The output file has the same time definition as the input file.
In this example, it means: 6000 samples between 0 and 6s.
Group in time bands: This option adds a step of computation. First it computes the TF decomposition for all the input file, then averages the power by time band. To define a time band:
- Enter your own time bands in the text area, one line per time band, with the following format: "name / time definition / function"
- Click on the button [Generate] to automatically create a list of time bands with the same length. You will be asked the maximal length of each time band.
- The function is the measure we take to combine the values for all the individual frequencies into one for the frequency band. Possible values are: mean, max, std, median.

Frequency definition: Frequencies for which the power will be estimated at each time instant.

Linear: You can specify the frequencies with the Matlab syntax start:step:stop.
The default is "1:1:60", which produces 60 values [1, 2, 3, 4, ..., 59, 60].
Log: With the option start:N:stop, produces a list of N frequencies logarithmically scaled between "start" and "stop". For example "1:40:80" is converted to [1, 1.5, 2.1, 2.7, ..., 61.5, 65.8, 75, 80]
Group in frequency bands: As for the time definition, this option leads to a two-step process. First it computes the TF decomposition for several values in the frequency band, then it averages the power of TF coefficients per frequency band. To define a frequency band:
- One line per frequency band, with the format "name / frequency definition / function"
- The frequency definition is a Matlab expression evaluated with an eval() call. If the frequency definition contains only two values, Brainstorm adds two extra values in the middle so that the final averaged value is a bit more robust. Example of valid expressions:
  "2,4": Evaluates to [2,4], and then expands to the frequency vector [2, 2.66, 3.33, 4]
  "2:0.5:4": Evaluates to [2 2.5 3 3.5 4]
  "2, 2.5, 3, 3.5, 4": Evaluates to [2 2.5 3 3.5 4]
- The function is the measure we take to combine the values for all the individual frequencies into one for the frequency band. Possible values are: mean, max, std, median.

Morlet wavelet options

Central frequency: Frequency where the mother wavelet is designed. All other wavelets will be shifted and scaled versions of the mother wavelet
Time resolution (FWHM): Temporal resolution of the wavelet at the central frequency (in units of Full Width Half Maximum). Click [Display] to see the resolution of the wavelet for other frequencies.

Compute the following measure:

The convolution of the signal with complex Morlet wavelets returns the complex coefficients for each frequency/time/sensor. Typically, what we display is the power of the coefficients (square of the amplitude: abs(TF)²). You can choose if you want to apply this transformation or not.
Power: Computes the "power" transformation immediately after the TF decomposition. This discards the phase information, but produces files that are twice smaller and a lot easier to process.
Magnitude: Save the magnitude of the complex values instead of the power: abs(TF).
None: Save the TF coefficients as they are computed (complex values). This can be useful if you plan to use these decompositions for other purposes that require the phase.
Some combinations of options may disable this choice. If you select frequency bands, the program will have to compute the power before averaging the values, therefore "none" is not an option.

Display: Time-frequency map

Right on the new time-frequency file > Time-freq: One matrix (same as double-clicking).
This menu displays the time-frequency decomposition of the first (and only) signal. The Brainstorm window shows two new elements: the tab "Display" and the frequency slider.
We can easily identify the three horizontal bars as the three sinusoids in the simulated signals, and the trade-off between accuracy in time and accuracy in frequency. Click on the figure to move the time-frequency cursor and explore the two axes of the plane.
- 2Hz: High frequency resolution but poor time resolution (supposed to start sharply at 2s)
- 20Hz: Better time resolution but poorer frequency resolution (17-24Hz)
- 50Hz: Continuous over the 6s - Frequency resolution gets even worse (40-60Hz). It looks discontinuous because this oscillation has the same amplitude as the white noise we added in the signal (weight 0.4, relatively to the 2Hz oscillation).
Current frequency: Slider that shows the current frequency selected in all the figures.
Just like the time, the frequency selection is centralized and managed by one control only for all the figures. As a consequence, it is impossible to display TF files with different frequency definitions at the same time. This can be perceived as an annoying limitation, but it allows all the simultaneous displays to be consistent at anytime and makes the interface more intuitive to manipulate, with lower risks of mistakes in the interpretation of the different figures.
List of signals: This drop-down list shows the signal currently displayed in the selected figure. In this case, there is only one channel of data called "s1". It will be more useful later.
Hide edge effects: When this option is selected, the time-frequency coefficients that could not be properly estimated because of a lack of time samples are hidden. It allows you to see only the information that is really reliable. The lower the frequency, the longer the edge effects.
Smooth display: Re-interpolates the time-frequency maps on a finer grid to produce nicer plots.
Measure: Type of measure that is currently represented in the selected figure. The entries that are enabled depend on the type of data that is saved in the file. In this case, we saved directly the power of the wavelet coefficients in the file, we discarded the angle/phase information, so the "phase" option is disabled. The other options are: Magnitude = sqrt(power), Log = 10*log10(power)
Colormap: As explained in the previous tutorials, you can change the colormap by clicking+moving on the colorbar on the right of the figure. Double-click on the colorbar to restore the defaults.

Display: Mouse and keyboard shortcuts

Mouse shortcuts

Left-click: Selection of current time and frequency.
Left-click + move: Select a time/frequency range. The legends of the X/Y axis show the selection.
Mouse wheel: Zoom in time, centered on the current time cursor.
Control + mouse wheel: Zoom in frequencies, centered on the current frequency cursor.
Right-click + move, or Control + left-click + move: Move in the zoomed image.
Double-click: Restore initial view.

Keyboard shortcuts:

Left/right arrows: Change the current time.
Page-up/page-down: Change the current time, 10 time samples at a time.
Up/down arrows: Change the the sensor displayed in this figure.
Control + up/down arrows: Change the current frequency.
Enter: View the original time series for this sensor.
Control + R: View the original MEG recordings.
Control + T: View the time-frequency 2D topography.
Control + I: Save as image.
Control + D: Dock figure in the Matlab environment.

Figure popup menu

Set selection manually: Does the same thing as drawing a time/freq selection square on a figure, but by typing the values for time and frequency manually.
Export to database: Save the selection for the displayed sensor in a new file in the database.
Export to file: Same as "Export to database", but the saved file is not registered in the database.
Export to Matlab: Same as "Export to database", but the output structure is sent to a variable in the Matlab base workspace instead of being saved to a file.

Display: Power spectrum and time series

Right-click on the file in the database or directly on the time-frequency figure to access these menus.

Power spectrum: For the current time, shows the power for all the frequencies.
Time series: For the current frequency, shows the power for all the time samples.
Example: Power spectrum density at 0.5s and power time series at 2Hz.
We see the oscillation at the 50Hz in the PSD plot, and the oscillation at 2Hz in the TS plot.
Example: Power spectrum density at 4s and power time series at 20Hz.
We see all three oscillations in the PSD plot, and the oscillation at 20Hz in the TS plot.
Note that if you right-click on the file in the database explorer and then select one of these menus, it will show all the signals. If you right-click on an existing time-frequency figure, it will show only the selected signal. It doesn't make any difference here because there is only one signal, but it will with the MEG recordings.

Normalized time-frequency maps

The brain is always active, the MEG/EEG recordings are never flat, some oscillations are always present in the signals. Therefore we are often more interested in the transient changes in the power at certain frequencies than at the actual power. A good way to observe these changes is to compute a deviation of the power with respect with a baseline.

There is another reason for which we are usually interested in standardizing the TF values. The power of the time-frequency coefficients are always lower in the higher frequencies than in the lower frequencies, the signal carries a lot less power in the faster oscillations than in the slow brain responses. This 1/f decrease in power is an observation that we already made with the power spectrum density in the filtering tutorial. If we represent the TF maps with a linear color scale, we will always see values close to zero in the higher frequency ranges. Normalizing each frequency separately with respect with a baseline helps obtaining more readable TF maps.

No normalization

The values we were looking at were already normalized (option: "1/f compensation" in the process options), but not with respect to a baseline. We will now compute the non-normalized power obtained with the Morlet wavelets and try the various options available for normalizing them.

In Process1, keep the simulated signal selected.
Run process: Frequency > Time-frequency (Morlet wavelets), "Spectral flattening: None"
Double-click on the file. As expected, we only see the lower frequencies in this representation: the power of the 2Hz oscillation is a lot larger than the power at 20Hz or 60Hz.

Spectral flattening

In Process1: Select this new non-normalized file "Power,1-60Hz".
Run process: Standardize > Spectral flattening, Method=1/f compensation.
This produces exactly the same results as previously (option "1/f compensation" in the time-frequency process). It multiplies the power at each frequency bin with the frequency value (eg. multiplies by 20 the power at 20Hz), in order to correct for the 1/f shape we observe in the power spectrum. This works well for the lower part of the spectrum and up to 60-80Hz, but past this range it tends to overcompensate the higher frequencies.
Note that it does not do any form of baseline correction: the 50Hz oscillation is visible everywhere.

Baseline normalization

The second way to proceed is to normalize the power with respect to its average level during a reference time period. We can consider the oscillations at 2Hz and 20Hz as our events of interest, and the 50Hz as noise we want to get rid of. The segment from 0 to 2 seconds does not contain any of the signals of interest, therefore we can consider it as a baseline.
However, we will not be able to use the full segment [0,2]s because of the edge effects we described at the beginning of this tutorial. The time-frequency map at 2Hz with the display option "Hide edge effects" (left figure below) shows that the power could not be estimated correctly before 0.75s, therefore we shouldn't use it as a baseline. The power time series at 2Hz (right) shows that the power related with the 2Hz oscillation starts to increase significantly after 1.25s, therefore it's not really a "baseline" anymore. This leaves only the segment [0.75,1.25]s available.
At 20Hz, the expected transient effects are only 75ms long, therefore we could use a much longer baseline if we were not interested by the lower frequencies: [0.075, 1.925]s.
In this case, we have a very long "resting" time segment (2s), therefore the edge effects are not a problem for picking a time window for the baseline normalization. We will use the first time window mentioned, [0.75,1.25]s as it is long enough (500ms) to estimate the baseline mean power. In real-life cases, with shorter epochs, it is sometimes difficult to find an acceptable trade-off between the baseline duration and the exclusion of the edge effects, especially for the lower frequencies.
Run process: Standardize > Baseline normalization, Baseline=[0.75, 1.25]s.
Method=Event-related perturbation: ERS/ERD stands for "event-related synchronization / desynchronization", a widely used normalization measure for time-frequency power maps. It evaluates the deviation from the mean over the baseline, in percents: (x-mean)/mean*100.
Double-click on the new "ersd" file. The colormap type changed from "Timefreq" to "Stat2", which uses by default the "rwb" color set and shows relative values. Indeed, the ERS/ERD values can be positive or negative: the power at a given time sample can be higher or lower than the average level during the baseline. In the simple simulation we used, there is no power decrease at any frequency after 2s, so the strong values are mostly positive. However, if you look in the file, you would see that there are many small negative values (due to the random noise we added).
Note that the 50Hz disappeared because it was present during the baseline, while the 2Hz and 20Hz oscillations show high positive values (between 100% and 2000% increase from baseline).
Remember to select your baseline very carefully according to the frequency range you are interested in. See below examples obtained with different baselines: 0.075-1.925s, 0-2s, 0-6s.
Change the colormap to "jet" if you prefer, and adjust the colormap contrast as needed.
This video may help you understand better what are the implications of the baseline selection: http://www.mikexcohen.com/lecturelets/whichbaseline/whichbaseline.html

Advanced

Tuning the wavelet parameters

Time resolution

You can adjust the relative time and frequency resolution of your wavelet transformation by adjusting the parameters of the mother wavelet in the options of the process.

Increasing the option "time resolution" will produce longer wavelets at a given frequency, hence increase the frequency accuracy (lower Δf) and decrease the time accuracy (higher Δt). Expect longer edge effects.
Decrease the time resolution will produce shorter wavelets at a given frequency, hence decrease the frequency accuracy (higher Δf) and increase the time accuracy (higher Δt). Expect shorter edge effects.
You can modify one or the other parameter, what is important is the product of the two. All the following combinations fc/FWHM_t produce the same results because their product is constant: (1Hz,3s), (3Hz,1s), (6Hz,0.5s), (60Hz,0.05s)

Examples for a constant central frequency of 1Hz with various time resolutions: 1.5s, 4s, 10s.

Frequency axis

You can also obtain very different representations of the data by changing the list of frequencies for which you estimate the power. You can change this in the options of the process.

Examples: Log 1:20:150, Log 1:300:150, Linear 15:0.1:25

Advanced

Hilbert transform [TODO]

We can repeat the same analysis with the other approach available for exploring the simulated signal in the time-frequency plane. The process "Frequency > Hilbert transform" first filters the signals in various frequency bands with a band-pass filter, then computes the Hilbert transform of the filtered signal. The magnitude of the Hilbert transform of a narrow-band signal is a measure of the envelope of this signal, and therefore gives an indication of the activity in this frequency band.

No normalization

Let's compute the same three results as before: non-normalized, 1/f compensation, baseline normalization.

In Process1, select the simulated signal.
Run process: Frequency > Hilbert transform, Spectral flattening=None, Mirror signal before.
The mirroring option helps with controlling the amplitude of the edge effects.
In the advanced options panel, keep the default options: Default frequency bands and Power.
Double-click on the new file. The figure now has only 6 rows, one for each frequency band.
Note that the values saved in the file and displayed are now amplitude and not power anymore.
The non-normalized results are already easy to interpret:
- delta (2-4Hz): Includes the 2Hz oscillation, contribution starts at 2s
- beta (15-29Hz): Includes the 20Hz oscillation, contribution starts at 2s
- gamma1(30-59Hz): Includes the 50Hz oscillation, contribution starts at the beginning (0s)
Right-click on the file or figure > Time series. Example for delta and beta.

Normalization

In Process1, select the non-normalized Hilbert-based decomposition.
Run process: Standardize > Spectral flattening, Method=1/f compensation.
Run process: Standardize > Baseline normalization, Baseline=[0.75, 1.25]s, Method=ERS/ERD.
Display the two normalized files side by side.

Method specifications

Filters: [TODO]
Hilbert transformation: Using the Matlab's hilbert() function.
Extraction of the envelope: Power of the complex Hilbert transform, abs(hilbert(x))².

MEG recordings: Single trials

Let's go back to our auditory oddball paradigm and apply the concepts to MEG recordings. We will use all the trials available for one condition to estimate the average time-frequency decomposition.

Spectral flattening

In Process1, select all the deviant trials in Run#01.
Run process: Frequency > Time-frequency (Morlet wavelets), Spectral flattening: None.
In the advanced options, select: Log 1:40:150, Power, Save average time-frequency maps.
Save individual TF maps: This option stops the computation here and saves in the database one time-frequency file for each input file (40 files), with one TF map for each scout.
Save average TF maps: Instead of saving the TF for each file separately, it automatically computes the average of the power of all the TF. This is a good choice if you do not plan to use independently all the TF files, because it saves a lot of time and disk space.
Remove evoked response from each trial before computing TF: This option first computes the average of all the trials, then subtracts this average from each trial before computing the time-frequency decomposition. This brings the signals to a slightly more stationary state, which may help for the evaluation of the frequency contents.

Baseline normalization

Double-click on the new file Avg,Power,1-150Hz (MEG). Select "Hide edge effects".
In the drop-down list, select sensor MLP56 (the one with the strongest respons at 90ms).
Right-click on the TF figure > Time series.
Defining a baseline is now a lot trickier than with 6s-long simulated signal. The epochs are only 600ms long, and the power at many frequencies could not be estimated correctly. If we want all the values to be "good" after 0s, we cannot use anything below 15Hz. If we want to normalize the values, we have to go even higher: 30Hz if we want a baseline of 50ms before 0.
The epochs we use in this tutorial are too short to perform a correct time-frequency analysis. We should have imported at least 200ms more on each side, just for controlling the edge effects. You should always think carefully about the length of the epochs you import in your database if you are planning to run any form of frequency or time-frequency analysis.
For the purpose of illustrating the tools available in Brainstorm, we will keep on working with these short epochs. Let's try to do the best we can with what we have here. We could use a baseline of 50ms to get a correct estimation above 30Hz, but this is probably a bit too short. We propose to include some more of the baseline (75ms), hoping there are not major edge effects in this segment.
In Process1, select the average time-frequency file.
Run process: Standardize > Baseline normalization, Baseline=[-75, 0]ms, Method=ERS/ERD.
The new menus available to display this file are described in the next section.

Things to avoid

Avoid computing the time-frequency decomposition of the average of the trials, you would miss some of the induced response, the brain activity in higher frequencies that is not strictly time-locked to the stimulus, and not aligned in phase across trials. Always prefer the computation of the average of the time-frequency power maps of each trials, as we did here.
This is well documented in: Bertrand O, Tallon-Baudry C (2000).
Avoid using the "spectral flattening: 1/f compensation" with frequencies above 60 or 80Hz. As mentioned before, it may overcompensate the higher frequencies. If you apply this normalization on the file we just computed, you get figures that give you a false idea of dominant high frequencies:

Display: All channels

Three menus display the time-frequency of all the sensors with different spatial organizations. All the figures below represent the ERS/ERD-normalized average. They use the "jet" colormap, which is not the default configuration for these files. To get the same displays, change the colormap configuration:
right-click on the figure > Colormap: Stat2 > Colormap > jet.

All channels: All the maps are displayed one after the other, in the order they are saved in the file.
2D Layout (maps): Show each TF map where the sensor is located on a flattened 2D map. Most display options are available, such as the colormap management and the option "Hide edge effects".
2D Layout (no overlap): Similar to the the previous display, but the positions of the images are reorganized so that they do not overlap.

Useful shortcuts for these figures:

Click: Clicking on any small TF image opens a new figure with only the selected sensor.
Shift + click: Opens the original recordings time series of the selected sensor, when available. Here, we display an average of time-frequency maps, so this menu has no effect.
Mouse wheel: Zoom in/out.
Right click + move: Move in a zoomed figure.

Display: Topography

The menus below show the distribution of TF power over the sensors, for one time point and one frequency bin, very similarly to what was introduced in tutorial Visual exploration.

2D Sensor cap / 2D Disc / 3D Sensor cap: 175ms, 8Hz
2D Layout: 8Hz (black), 35Hz (white)

Useful shortcuts for these figures:

Left/right arrows: Change the current time.
Up/down arrows: Change the current frequency.
Control + E: Display the sensors markers/names.
Shift + click on a sensor: Displays the time-frequency decomposition for that specific sensors.
Right click + move: Select a group of sensors.
Shift + scroll: Change the gain of the time series (2D Layout).
Control + scroll: Change the length of the window displayed around the current time (2D Layout).

Advanced

Scouts

Similar calculations can be done at the level of the sources, either on the full cortex surface or on a limited number of regions of interests. We will start with the latter as it is usually an easier approach.

Drag and drop all the deviant trials from both runs, select [Process sources].
Run process "Frequency > Time-frequency (Morlet wavelets)".
Select the option "Use scouts" and select all the scouts defined in the previous tutorial.
In the advanced options, select "Scout function: After" and "Output: Save average".
Run the process (it may take a while).
The scout function was introduced in the previous tutorial. It is the method we use to group the time series for the 20 dipoles we have in each scout into one unique signal. When computing the TF of one scout, we have the choice between applying this function before or after the time-frequency decomposition itself.
- Before: Extract the 20 source signals, apply the scout function to get one signal, run the TF decomposition of this signal. This is faster but may lead to information loss.
- After: Extract the 20 source signals, run the TF decomposition of the 20 signals, apply the scout function on the power of the TF maps. Always prefer this option when possible.
Rename the new file to add a tag "Deviant" in it. Then right-click > Time-freq: All scouts.
In Process1, select the new average TF file.
Run process: Standardize > Baseline normalization, Baseline=[-75, 0]ms, Method=ERS/ERD.

Advanced

Full cortical maps

Computing the time-frequency decomposition for all the sources of the cortex surface is possible but complicated because it can easily generate gigantic files, completely out of the reach of most computers. For instance the full TF matrix for each trial we have here would be [Nsources x Ntimes x Nfrequencies] = [15000 x 361 x 40] double-complex = 3.2 Gb!

We have two ways of going around this issue: computing the TF decomposition for fewer frequency bins or frequency bands at a time, or as we did previously, use only limited number of regions of interest.

In Process1, keep all the deviant trials both conditions selected, select [Process sources].
Run process "Frequency > Hilbert transform", Spectral flattening=None, Mirror signal before.
To process the entire brain, do not select the option "Use scouts".
In the advanced options, select "Optimize storage: No", this option is not available when computing on the fly the average of multiple trials. Save the power, Save the average Hilbert maps.
Optimize the storage of the time frequency file: Let's describe this option in more details.
- When computing the TF decomposition of a source file, we are actually applying sequentially two linear transformations to the original recordings: the TF analysis and the source inversion. These two processes can be permuted: TF(Inverse(Recordings)) = Inverse(TF(Recordings)).
- Therefore we can optimize the TF computation time by applying the wavelet transformation only to the sensor recordings, and then multiply the wavelet complex coefficients by the inverse operator (ImagingKernel). This trick is always used in the computation of the Hilbert and Morlet transforms.
- When we have the option to the save the complex values (constrained sources and no averaging), this can also be used to optimize the storage of the files. In these cases, we save only the wavelet transformation of the sensor data. Later, when the file is loaded for display, the imaging kernel is applied on the fly. This can be disabled explicitely with this option.
Rename the new Hilbert file to include the tag "Deviant", and select it in Process1.
Run process: Standardize > Baseline normalization, Baseline=[-75, 0]ms, Method=ERS/ERD.
Right-click on the Hilbert file > Display on cortex.
The frequency slider now shows frequency bands ("alpha:8-12Hz") instead of frequencies ("12Hz"). You can explore the source activity in time and frequency dimensions. The screen capture below shows the activity at 175ms: a 60% increase in the alpha band around the auditory cortex and a 20% decrease in the beta band around the motor cortex.
Shift + click on the cortex surface: Displays the TF decomposition of the selected source.
Right-click on the brain: Selects the closest vertex and displays the popup menu at the same time. The first three menus are relative to the source that was just clicked.

Advanced

Unconstrained sources

In the current example, we are working with the simple case: sources with constrained orientations. The unconstrained case is more difficult to deal with, because we have to handle correctly the three orientations we have at each vertex.

Full cortex: Computes the TF decompositions for all the sources (3*15000=45000), then sum at each location the power for the three orientations.
Scouts: Option "Scout function" in the process.
- Before: Extract the 20*3=60 source signals, apply the scout function to get three signals (one per orientation), run the TF decomposition of the three signals, and finally sum the power of the three TF maps. This is faster but may lose some frequency resolution (especially for constrained sources).
- After: Extract the 20*3=60 source signals, run the TF decomposition of the 60 signals, apply the scout function on the power of the TF maps for each orientation separately, and finally sum the power obtained for the three orientations.
The storage optimization option is not available with unconstrained sources.

Advanced

Getting rid of the edge effects

To avoid making mistakes in the manipulation of the data and producing more readable figures, we encourage you to cut out the edge effects from your time frequency maps after computation.

In Process1, select the very first computed in this tutorial: Test/Simulation/Power,1-60Hz | multiply
Run the process: "Extract > Extract time", Time window = [0.75, 5.25]s
Open the new file, select the option "Hide edge effects": Almost everything left in this new file is correctly estimated. Brainstorm keeps track of the edge effects in the TFmask field of the file.
We recommend you do the same when epoching your recordings: import trials that are longer than necessary, and after the time-frequency estimation, remove the unnecessary segments.

Advanced

On the hard drive

Right click one of the first TF file we computed > File > View file contents.

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

TF: [Nsignals x Ntime x Nfreq] matrix containing all the values of the time-frequency decomposition (complex wavelet coefficients, or double values for power/magnitude/Z-score).
TFmask: [Nfreq x Ntime] logical mask indicating the edge effects (0=edge, 1=valid value).
Std: [Nsignals x Ntime x Nfreq] standard deviation if this file is an average.
Comment: String displayed in the database explorer to represent the file.
DataType: From what kind of data this file was computed: {'data', 'results', 'scout, 'matrix'}
Time: [1 x Ntime] Time vector used to estimate this file.
TimeBands: [Ntimebands x 3] Cell array where each line represents a time band:
{'band_name', 'time definition', 'function'}
Freqs: For regular frequency binning: vector containing all the frequencies.
If using frequency bands: [Nfreqbands x 3] cell array, where each line represents a frequency band {'band_name', 'frequency definition', 'function'}
RefRowNames: Used only for connectivity matrices.
RowNames: [1 x Nsignals] Cell array of strings that describes each row of the TF matrix. In this specific case, it would be the list of all the MEG sensor names. But it could also be a list of names of scouts or clusters.
Measure: Contains the name of the function that was applied right after the computation of the wavelet coefficients. So it represents the type of data contained in the TF matrix. Possible values:
- none: No measure applied, TF contains the complex wavelet coefficients.
- power: Power for each frequency, ie. the square of the amplitude: abs(coefficients)²
- magnitude: abs(coefficients)
- log: 10 * log10( abs(coefficients)²)
- phase: angle(coefficients)
Method: String that identifies the process that generated the file:
{'morlet', 'fft', 'psd', 'hilbert', 'corr', 'cohere', 'granger', 'plv'}
DataFile: Initial file from which this file was computed. In the database explorer, the TF file will be shown as a child of this DataFile file.
SurfaceFile / GridLoc / GridAtlas: Source space that was used, only for source files.
nAvg: Number of trials that were averaged to obtain this file.
ColormapType: String, force a specific colormap type to be used when displaying this file.
DisplayUnits: String, force to use specific units when displaying this file.
Options: Options that were selected in the time-frequency options window.
History: List of operations performed on this file (menu File > View file history).

Useful functions

in_bst_timefreq(TimefreqFile): Read a time-frequency file.
in_bst(FileName, TimeWindow): Read any Brainstorm file with the possibility to load only a specific part of the file. "TimeWindow" is an range of time values in seconds: [tStart, tStop].
bst_process('LoadInputFile', FileName, Target, TimeWindow): The most high-level function for reading data files. "Target" is a string with the list of sensor names or types to load.
morlet_transform(): Applies complex Morlet wavelet transform to the time series in input.

Advanced

References

Bertrand O, Tallon-Baudry C (2000)
Oscillatory gamma activity in humans: a possible role for object representation
Int J Psychophysiol, 38(3):211-23
Bruns A (2004)
Fourier-, Hilbert- and wavelet-based signal analysis: are they really different approaches?
J Neurosci Methods, 137(2):321-32
Le Van Quyen M, Foucher J, Lachaux J, Rodriguez E, Lutz A, Martinerie J, Varela FJ (2001)
Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony
J Neurosci Methods, 111(2):83-98
Pfurtscheller G (1992)
Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest
Electroencephalogr Clin Neurophysiol, 83(1):62-9
Series of lectures by Mike X Cohen:
http://www.mikexcohen.com/lectures.html

Advanced

Additional documentation

Forum: Time and frequency resolution: http://neuroimage.usc.edu/forums/showthread.php?1848

-  ⇤ ← Revision 23 as of 2015-08-19 23:41:55 → 
  Size: 32091
  Editor: FrancoisTadel
  Comment:
+   ← Revision 127 as of 2016-05-04 22:57:55 → ⇥
  Size: 46200
  Editor: FrancoisTadel
  Comment:
-Deletions are marked like this.
+Additions are marked like this.
 Line 2:
-''Authors: Francois Tadel, Dimitrios Pantazis, Sylvain Baillet''

This tutorial introduces how to compute time-frequency decompositions of MEG/EEG recordings and cortical currents using complex Morlet wavelets and Hilbert transforms.
+''Authors: Francois Tadel, Dimitrios Pantazis, Elizabeth Bock, Sylvain Baillet''

This tutorial introduces how to compute time-frequency decomposition of MEG/EEG recordings and cortical currents using complex Morlet wavelets and Hilbert transforms.
 Line 11:
-The averaging in time domain may also lead to a cancellation of these oscillations when they are not strictly locked in phase across trials. Averaging trials in time-frequency domain allows to extract the power of the oscillation regardless of the phase shifts. For a better understanding of this topic, we recommend the lecture of the following article: Bertrand O, Tallon-Baudry C (2000), [[http://www.sciencedirect.com/science/article/pii/S0167876000001665|Oscillatory gamma activity in humans: a possible role for object representation]]

In Brainstorm we offer two approaches for computing time-frquency decompositions (TF): the first one is based on the convolution of the signal with series of complex Morlet wavelets, the second filters the signal in different frequency bands and extract the envelope of the filtered signals using the Hilbert transform.

 . {{attachment:slide.gif}}
+The averaging in time domain may also lead to a cancellation of these oscillations when they are not strictly locked in phase across trials. Averaging trials in time-frequency domain allows to extract the power of the oscillation regardless of the phase shifts. For a better understanding of this topic, we recommend the lecture of the following article: Bertrand O, Tallon-Baudry C (2000), [[http://www.iec-lnc.ens.fr/IMG/Files/Visual/Publications/2000_Review_IntJPsychol.pdf|Oscillatory gamma activity in humans: a possible role for object representation]]

In Brainstorm we offer two approaches for computing time-frequency decomposition (TF): the first one is based on the convolution of the signal with series of complex Morlet wavelets, the second filters the signal in different frequency bands and extracts the envelope of the filtered signals using the Hilbert transform.

{{attachment:slide.gif}}
 Line 20:
-Contrary to the standard short-time Fourier transform, wavelets have variable resolution in time and frequency. When designing the wavelet, we basically decide a trade off between temporal and spectral resolution.

To design the wavelet, we first need to choose a ''central frequency'', ie the frequency where we will define the mother wavelet. All other wavelets will be scaled and shifted versions of the mother wavelet. Unless interested in designing the wavelet at a particular frequency band, the default 1Hz should be fine.

Then, the desirable'' time resolution'' for the central frequency should be defined. For example, we may wish to have a temporal resolution of 3 seconds at frequency 1 Hz (default parameters). These two parameters, uniquely define the temporal and spectral resolution of the wavelet for all other frequencies, as shown in the plots below.

Resolution is given in units of Full Width Half Maximum of the Gaussian kernel, both in time and frequency. The relevant plots are given below.

{{attachment:waveletOptions.gif||height="291px",width="637px"}}
+Contrary to the standard short-time Fourier transform, wavelets have variable resolution in time and frequency. For low frequencies, the frequency resolution is high but the time resolution is low. For high frenqucies, it's the opposite. When designing the wavelet, we basically decide a trade-off between temporal and spectral resolution.

 . {{attachment:tf_resolution.gif||height="144",width="336"}}

To design the wavelet, we first need to choose a ''central frequency'', ie. the frequency where we will define the mother wavelet. All other wavelets will be scaled and shifted versions of the mother wavelet. Unless interested in designing the wavelet at a particular frequency band, the default 1Hz should be fine.

Then, the desirable'' time resolution'' for the central frequency should be defined. For example, we may wish to have a temporal resolution of 3 seconds at frequency 1 Hz (default parameters). These two parameters, uniquely define the temporal and spectral resolution of the wavelet for all other frequencies, as shown in the plots below. Resolution is given in units of Full Width Half Maximum of the Gaussian kernel, both in time and frequency.

 . {{attachment:waveletOptions.gif||width="650px"}}
 Line 33:
-From the figure above, which designs the Morlet wavelet, we can see that the default wavelet (central frequency 1Hz, FWHM=3sec) has temporal resolution 0.25sec at 5Hz and 0.1sec at 10Hz. In such case, the edge effects are roughly half these times: 125ms in 5Hz and 50ms in 10Hz. Examples of such edge effects are given in the figures below.
+From the figure above, which designs the Morlet wavelet, we can see that the default wavelet (central frequency '''fc'''=1Hz, '''FWHM_tc'''=3sec) has temporal resolution 0.6sec at 5Hz and 0.3sec at 10Hz. In such case, the edge effects are roughly half these times: 300ms in 5Hz and 150ms in 10Hz.

More precisely, if '''f''' is your frequency of interest, you can expect the edge effects to span over FWHM_t seconds: '''FWHM_t = FWHM_tc * fc / f / 2'''. Examples of such transients are given in the figures below.
-Line 37:
+Line 39:
-== Simple example ==
Let's start with a simple example to explore the options of the time-frequency decomposition process.

 * In Process1, select the average recordings for the '''deviant '''condition for Run#01.
 * Select process: '''Frequency > Time-frequency (Morlet wavelets)'''.<<BR>><<BR>> {{attachment:example_process.gif||height="288",width="507"}}
 * Enter "'''MEG'''" for the sensor types, do '''not '''select the "Spectral flattening" option.<<BR>>Click on the button ['''Edit'''] to see all the process options.<<BR>>Time definition: '''Same as input''', Frequency definition: '''Log 1:40:150''', Compute measure: '''Power'''. <<BR>><<BR>> {{attachment:example_options.gif||height="528",width="392"}}
+== Simulation ==
We will illustrate the time-frequency decomposition process with a simulated signal.

 * The following code generates a sum of three sinusoids (2Hz, 20Hz, 50Hz) with random white noise. The 50Hz and noise are present everywhere, the 2Hz and 20Hz start only after two seconds.<<BR>>
 {{{
f1 = 2; f2 = 20; f3 = 50;
i =2000:6000;
Data(1,i) = sin(f1*2*pi*t(i)) + 0.4 * cos(f2*2*pi*t(i));
Data = Data + 0.2 * sin(f3*2*pi*t) + 0.4 * rand(1,6000);
}}}

 * Empty the Process1 list (right-click > Clear list) then click on [Run].
 * Run process: '''Simulate > Simulate generic signals'''. <<BR>>Ntime='''6000''', Sampling frequency='''1000Hz''' (signal duration = 6000/1000 = 6 seconds).<<BR>>Copy-paste the few lines of code above to generate the sum of three sinusoids.<<BR>><<BR>> {{attachment:simulate_process.gif||height="412",width="489"}}
 * Double-click on the new file to look at the simulated signal.<<BR>><<BR>> {{attachment:simulate_display.gif||height="159",width="667"}}
 * In Process1, select the simulated signal.
 * Run process: '''Frequency > Time-frequency (Morlet wavelets)'''. <<BR>>Select the option '''1/f compensation''': The normalization will be discussed later.<<BR>><<BR>> {{attachment:example_process.gif||height="271",width="557"}}
 * Click on the button ['''Edit'''] to see all the process options.<<BR>>Time definition: '''Same as input''', Frequency definition: '''Linear 1:1:60''', Compute measure: '''Power'''. <<BR>><<BR>> {{attachment:example_options.gif||height="506",width="391"}}
-Line 49:
+Line 62:
- * '''Same as input file''': The output file has the same time definition as the input file. <<BR>>In this example, it means: 361 samples between -100ms and 500ms.
+ * '''Same as input file''': The output file has the same time definition as the input file. <<BR>>In this example, it means: 6000 samples between 0 and 6s.
-Line 51:
+Line 64:
-  * Enter your own time bands in the text area, one line per time band, with following format: '''"name / time definition / function" '''
  * Click on the "Generate" button to automatically  create a list of time bands with the same length. You will be asked the maximal length of each time band.
+  * Enter your own time bands in the text area, one line per time band, with the following format: '''"name / time definition / function" '''
  * Click on the button [Generate] to automatically  create a list of time bands with the same length. You will be asked the maximal length of each time band.
-Line 55:
+Line 68:
-'''Frequency definition''': Frequencies for which the power will be estimated at each time instant
+'''Frequency definition''': Frequencies for which the power will be estimated at each time instant.
-Line 58:
+Line 71:
+ * '''Log''': With the option start:N:stop, produces a list of N frequencies logarithmically scaled between "start" and "stop". For example "1:40:80" is converted to [1, 1.5, 2.1, 2.7, ..., 61.5, 65.8, 75, 80]
-Line 60:
+Line 74:
-  * The frequency definition is a Matlab expression evaluated with an eval() call. If the frequency definition contains only two values, Brainstorm adds two extra values in the middle so that the final averaged value is a bit more robust. Example of valid expressions:<<BR>>"'''2,4'''": Evaluates to [2,4], and then expands to the frequency vector [2, 2.66, 3.33, 4]<<BR>>"'''2:0.5:4'''": evalues to [2 2.5 3 3.5 4]<<BR>>"'''2, 2.5, 3, 3.5, 4'''": Evaluates to [2 2.5 3 3.5 4]
+  * The frequency definition is a Matlab expression evaluated with an eval() call. If the frequency definition contains only two values, Brainstorm adds two extra values in the middle so that the final averaged value is a bit more robust. Example of valid expressions:<<BR>>"'''2,4'''": Evaluates to [2,4], and then expands to the frequency vector [2, 2.66, 3.33, 4]<<BR>>"'''2:0.5:4'''": Evaluates to [2 2.5 3 3.5 4]<<BR>>"'''2, 2.5, 3, 3.5, 4'''": Evaluates to [2 2.5 3 3.5 4]
-Line 70:
+Line 84:
- * The convolution of the signal with complex Morlet wavelets returns the complex coefficients for each frequency/time/sensor. Typically, what we display is the power of the coefficients (square of the amplitude: abs(TF).^2). You can choose if you want to apply this transformation or not.
+ * The convolution of the signal with complex Morlet wavelets returns the complex coefficients for each frequency/time/sensor. Typically, what we display is the power of the coefficients (square of the amplitude: abs(TF)^2^). You can choose if you want to apply this transformation or not.
-Line 72:
+Line 86:
- * '''None''': Save the TF coefficients as they are computed. This can be useful if you plan to use these decompositions for other purposes that require the phase.
+ * '''Magnitude''': Save the magnitude of the complex values instead of the power: abs(TF).
 * '''None''': Save the TF coefficients as they are computed (complex values). This can be useful if you plan to use these decompositions for other purposes that require the phase.
-Line 75:
+Line 90:
-== Display: One channel ==
 * Right on the new time-frequency file > '''One channel''' (same as double-clicking on the file). <<BR>><<BR>> {{attachment:example_file.gif||height="250",width="341"}}

 * This menu displays the time-frequency decomposition of the first sensor. The Brainstorm window shows two new elements: the tab "Display" and the  "current frequency" slider.<<BR>><<BR>> {{attachment:example_display.gif||height="238",width="569"}}
+== Display: Time-frequency map ==
 * Right on the new time-frequency file > '''Time-freq: One matrix''' (same as double-clicking). <<BR>>This menu displays the time-frequency decomposition of the first (and only) signal.  The Brainstorm window shows two new elements: the tab "Display" and the frequency slider.<<BR>><<BR>> {{attachment:example_file.gif||height="210",width="660"}}
 * We can easily identify the three horizontal bars as the three sinusoids in the simulated signals, and the trade-off between accuracy in time and accuracy in frequency. Click on the figure to move the time-frequency cursor and explore the two axes of the plane.
  * '''2Hz''': High frequency resolution but poor time resolution (supposed to start sharply at 2s)
  * '''20Hz''': Better time resolution but poorer frequency resolution (17-24Hz)
  * '''50Hz''': Continuous over the 6s - Frequency resolution gets even worse (40-60Hz). It looks discontinuous because this oscillation has the same amplitude as the white noise we added in the signal (weight 0.4, relatively to the 2Hz oscillation).
-Line 81:
+Line 99:
- * '''List of sensors''': This drop-down list shows the current entry displayed in the selected figure. In this case, it is by default the MEG sensor "MLC11". Try to change it, it will update the display of the figure, showing the frequency power for another sensor.
 * '''Hide edge effects''': When this option is selected, the time-frequency coefficients that could not be properly estimated because of a lack of time samples are hidden. It allows you to see only the information that is really reliable. In this example, you can see that it masks a large portion of the window: 600ms is not enough to perform a correct time-frequency decomposition. Try to work on larger time windows when using these tools.<<BR>><<BR>> {{attachment:example_hide.gif||height="181",width="498"}}
+ * '''List of signals''': This drop-down list shows the signal currently displayed in the selected figure. In this case, there is only one channel of data called "s1". It will be more useful later.
 * '''Hide edge effects''': When this option is selected, the time-frequency coefficients that could not be properly estimated because of a lack of time samples are hidden. It allows you to see only the information that is really reliable. The lower the frequency, the longer the edge effects.<<BR>><<BR>> {{attachment:example_hide.gif||height="207",width="555"}}
-Line 85:
+Line 103:
+ * '''Colormap''': As explained in the previous tutorials, you can change  the colormap by  clicking+moving on the colorbar on the right of the  figure. Double-click on the colorbar to restore the defaults.
-Line 89:
+Line 108:
- * '''Left-click + move''': Select a time/frequency range. The legends of the X/Y axis show the selection. <<BR>><<BR>> {{attachment:example_select.gif||height="134",width="252"}}
+ * '''Left-click + move''': Select a time/frequency range. The legends of the X/Y axis show the selection. <<BR>><<BR>> {{attachment:example_select.gif||height="159",width="362"}}
-Line 95:
+Line 114:
- * '''Colormap''': As explained in the previous tutorials, you can change the colormap by  clicking+moving on the colorbar on the right of the figure. Double-click on the colorbar to restore the defaults.
-Line 109:
+Line 127:
- . {{attachment:example_popup.gif}}
+ . {{attachment:example_popup.gif||height="191",width="480"}}
-Line 116:
+Line 134:
-== Display: Multiple channels ==
Three menus display the same information (the TF maps of all the channels) with different spatial organizations.

 * '''Time-freq: All channels''': All the maps are displayed one after the other, in the order they are saved in the file. <<BR>><<BR>> {{attachment:timeFreqAllSensors.gif}}
 * '''Time-freq: 2D Layout''': Similarly to the "2D Layout" for the recordings, it shows each the position of each TF map where the sensor is located on a flattened 2D map. <<BR>><<BR>> {{attachment:timeFreqAllSensors2.gif}}
 * '''Time-freq: 2D Layout (no overlap)''': Similar to the "2D Layout", but reorganizing the positions of the small images so that they do not overlap. Overlapping in the "2D Layout" views is not a problem with this old CTF system (151 sensors); it gets much worse with the newer CTF systems (275 sensors), or the Elekta-Neuromag systems (102 chips with 3 sensors at each position). This display is a good intermediate between seeing all the information, and having it spatially organized. <<BR>><<BR>> {{attachment:timeFreqAllSensors3.gif}}

What to do with these figures:

 * Maximize the window to be able to see something
 * '''Click''': clicking on any small TF image opens a new figure to display only the selected sensor
 * '''Shift + click''': Same, but opens the original recordings time series of the sensor instead of the TF
 * '''Mouse wheel''': Zoom / unzoom
 * '''Left+right mouse click + move''': Move in a zoomed figure.

== Display: Sensor topography ==
=== 2D topography ===
Right-click on the TF file in condition ''Left ''and  select successively the first three topography menus. All  these three  windows represent the same information, in a  slightly  different way: a  spatial map of the power of the current  frequency, for  all the  sensors at the current time.

 * 3D Sensor cap
 * 2D Sensor  cap
 * 2D Disc

Try  to  move the time and frequency sliders and see what happens. You  can  open at the same time a "One channel" view, to keep track  easily of the  current time and frequency.

 . {{attachment:timeFreqTopo.gif}}

Keyboard shortcuts:

 * '''Left  and right arrows''': Change the current time
 * '''Page-up and page-down''': Also change the current time, but  faster  (10 time samples instead of one)
 * '''Up and down arrows''':  Change the current frequency
 * '''Control  + E''': Display the  sensors markers
 * '''Control  + E again''': Display  the sensors names

 * '''Shift  + Click on a sensor''':  Displays the time-frequency decomposition for that specific  sensors. {{attachment:topoSelectSensor.gif}} {{attachment:topoSelectSensorTf.gif}}

=== 2D Layout ===
The  last display mode available for these TF decomposition of recordings is  this "2D Layout" topography menu. Right-click on the TF file for the  Left condition and select "2D Layout".

This  represents spatially the power of the current frequency for all the   sensors and all the time points. Try to move the current frequency   slider to see how the display changes when increasing the frequency. It   is a good example to show that the time resolution increase with the  frequency. Below: the "2D Layout" for f=8Hz and f=60Hz.

 . {{attachment:topo2DLayout1.gif}} {{attachment:topo2DLayout2.gif}}

Useful operations for this window:

 * '''Mouse wheel''': Zoom / unzoom
 * '''Left+right click + move''': (or middle click+move) Move into the zoomed window
 * '''Shift+mouse wheel''': Increase/decrease the sharpness of the time series
 * '''Control + mouse wheel''': Increase/decrease the length of the time window displayed around the current time.
 * '''Click on a line''': Select a sensor
 * '''Right click + move''': Select a group of sensors
 * '''Shift + click on a line''': Display the time-frequency decomposition for the selected sensor
 * '''Right-click''': Display popup menu. More options in the "2D Layout options" sub-menu.
 * '''Control + E''': Display the channels names
 * All the other shortcuts and the popup  menu are already described  in [[Tutorials/TutExploreRecodings|tutorial   #4  "Exploring the recordings"]].
-Line 174:
+Line 135:
-These two menus generate similar figures: they represent as a line the evolution of each sensor for one parameter (time or frequency) when the other parameter is fixed.

 * Power spectrum: For the given time, show the evolution of all the sensors across the frequencies
 * Time series: For the given frequency, show the evolution of all the sensors in time

Open at the same time three figures for the Left/ERF file:

 * One channel
 * Power spectrum
 * Time series.

Then navigate in time and frequencies, and observe how each figure gets updated at each change.

 . {{attachment:timeFreqSpectrum.gif}}



== Scouts time series ==
 * The  previous section introduces an interesting tool to display the TF for a  single source, but it is hard to extract meaningful information that  can be analyzed statistically or compared between subjects this way. A  better option is to define cortical regions of interest (scouts, see [[Tutorials/TutScouts|tutorial #8]]), and compute the TF for the average time series of these scouts.
 * Define the scouts:<<BR>>
  * Double click for the source file in the Left condition.
  * Make  sure that your LeftSS and RightSS scouts are still available in your  cortex surface.If not, go back to tutorial 8, and recreate them <<BR>><<BR>> {{attachment:createScouts.gif}}

 * Compute the TF decomposition for these scouts:
  * Drag and drop the Left source file (''MN: MEG(Constr)'') into ''Process1''.
  * Select the "Process sources" button.
  * Click on Run. Select "Frequency > Time-frequency (Morlet wavelets)".
  * As some scouts are available for the input, you will see this new box:<<BR>><<BR>> {{attachment:selectScouts.gif}}
  * Select "Use scouts time series", and select both scouts in list.
  * Click  on Edit. This time, most of the options are available again. We are  going to compute the full TF matrix at once for the two scouts, so we  can compute the power right away, together with time and/or frequency  bands.
  * Leave all the options to the default (linear time, linear frequencies, scout function: After), then click on Ok, then Run.
  * A new file appears in the tree, as a child of the sources file: "Scouts,Power,1:1:60Hz".
 * To compute the TF of each scout, the  operations performed by the program are the following:
  * Get time series for each source in the scout
  * If the scout function is applied __before__ the TF:
   * Apply the "scout function", that combines the different sources to get only one time series for each scout.
   * Compute the power of the TF for each scout.
  * If the scout function is applied __after__ the TF:
   * Compute the power of the TF for each source
   * Apply the "scout function", that combines the different sources TF to  get one TF only  for each scout.
  * The selection of the "scout function" is done in the "Scout"  tab.

 * '''Output''':
  * '''Save individual TF maps''': This option stops the computation here, and saves in the database one time-frequency file for each input file.
   * Be very careful using this option because if you have hundreds of trials, it will fill your hard drive very quickly.
   * One full TF file in this case contains [nSensors x nTime x nFreq] = 3.397.500 values = 25 Mb.
   * You  can however limit the amount of data produced by using time and/or  frequency bands. For 15 time bands and 6 frequency bands, the size of  each file drops to: 110 Kb.
  * '''Save average TF maps''': This option is available only when there are more than one input files (number of files you dropped in the Processes tab).
   * Instead of saving the TF for each file separately, it automatically computes the average of the power of all the TF.
   * This is a good choice if you do not plan to use independently all the TF files, because it saves a lot of time and disk space.

 * Right-click on this file and try all the visualization menus available (including "Time series" and "Power spectrum").<<BR>><<BR>> {{attachment:scoutsTfFile.gif}}
+Right-click on the file in the database or directly on the time-frequency figure to access these menus.

 * '''Power spectrum''': For the current time, shows the power for all the frequencies.
 * '''Time series''': For the current frequency, shows the power for all the time samples.
 * Example: Power spectrum density at '''0.5s''' and power time series at '''2Hz'''. <<BR>>We see the oscillation at the 50Hz in the PSD plot, and the oscillation at 2Hz in the TS plot.<<BR>><<BR>> {{attachment:example_psd1.gif||height="148",width="660"}}
 * Example: Power spectrum  density at '''4s''' and power time series at '''20Hz'''. <<BR>>We see all three oscillations in the PSD plot, and the oscillation at 20Hz in the TS plot.<<BR>><<BR>> {{attachment:example_psd2.gif||height="148",width="661"}}
 * Note that if you right-click on the file in the '''database explorer''' and then select one of these menus, it will show '''all the signals'''. If you right-click on an existing '''time-frequency figure''', it will show only the '''selected signal'''. It doesn't make any difference here because there is only one signal, but it will with the MEG recordings.
-Line 228:
+Line 144:
- * Double-click on the new file to display the time-frequency decomposition of the first sensor. <<BR>><<BR>> {{attachment:example_display.gif||height="216",width="623"}}
 * Time-frequency files, once computed, can be processed with ''Process1'' and ''Process2'', exactly as recordings or sources file, as introduced in [[Tutorials/TutProcesses|Tutorial #9]].
 * Only the TF of sensors, clusters and scouts can be processed this way.
  * The   TF of full source files, as explained above, are saved in the database   in an incomplete form, that requires a few more on-the-fly computation   steps. The processing engine requires an access to the full TF matrix,   which is usually too big to be handled in Matlab for these TF/sources   files.
  * If you need to average or do some group-analysis on   time-frequency data at the source level, the problem has to be   simplified a few steps before: for instance by doing your analysis only   on a few scouts instead of the full cortex.
 * Drag and drop one of the TF files in Process1, select the "Time-frequency" button, click on Run.<<BR>><<BR>> {{attachment:panelProcess.gif}}
 * One   of the process that is commonly applied on the time-frequency maps is   the Z-score, to normalize the values in each frequency band with  respect  to a baseline. This way, all the frequency bands can be  compared  directly in terms of amplitude. If we look directly at the  power of the  Morlet coefficients, the values we observe are always  lower in the high  frequencies than in the low frequencies. The Z-score  fixes this and  helps getting figures that are easier to read.
 * Select   the process "Standardize > Compute Z-score (static)". Leave the   baseline to the default value (the pre-stimulus interval). Click on run.
 * Double-click on the new file.<<BR>><<BR>> {{attachment:timeFreqZscore.gif}}
 * Red   values are positive, representing an increase in power with respect   with the baseline, blue values are negative, representing a decrease.
 * The   time-frequency map may appear all white or all blue depending on the   way the colormap is configured. To edit the colormap, right-click on the   figure > Colormap.

== Cortical sources ==
 * How to display the TF decompositions for the source time series?
  * It is possible to estimate the TF for each source of the brain, but it would be unrealistic to save this information to a file. The size of the TF matrix would be [nVertices x nTimes x nFrequencies] = [15010 x 375 x 60] complex-double = 5.2 Gb!
  * We need then to simplify this problem. Two methods: compute the TF only for a few scouts (next section), or use the linear property of the TF decomposition.
  * As both the source reconstruction process (ImagingKernel * recordings) and the TF are linear operators, it is possible to exchange them:<<BR>>TF(Inverse(Recordings)) = Inverse(TF(Recordings))<<BR>>=> Power(TF(Inverse(Recordings))) =  Power(Inverse(TF(Recordings)))
  * So the solution is to estimate the TF of the recordings, and then multiply it on the fly by the ImagingKernel only for the required sources and/or frequency and/or time instant. This is done in a completely transparent way from the user point a view.
 * Start the computation of the TF for a source file:
  * Drag and drop the source file computed for the ''Left / ERF'' file in ''Process1''.
  * Select the "Process sources" button.
  * Click on Run. Select "Time-frequency (Morlet wavelets)". Do NOT check the box "Use scouts time series".
  * Click on Edit. There is a new option available, that offers between saving only the kernel and the TF of sensors, or the full TF matrix. Select the optimized solution.
  * Note that this solution is not compatible with any of the operations that requires to compute explicitly the power of all the TF decomposition: time bands and frequency bands.<<BR>><<BR>> {{attachment:timeFreqKernel.gif}}

 * Click on Ok, then Run. A new file appears in the database explorer, as a child of the source file. Right-click on it, try all the visualization menus.<<BR>><<BR>> {{attachment:computeForSources.gif}}
 * If you understood well how to use the visualization of both the sources and the time-frequency maps, the manipulation of these figures should be rather intuitive. Change several times the current time and the current frequency in different ways. Change current time and frequency so that you can observe the main  response in the somatosensory cortex (eg. t = 46.40ms, f = 45Hz). <<BR>><<BR>> {{attachment:timeFreqSources.gif}}
 * '''Shift + click '''on the cortex surface: This is a useful shortcut that you should remember. It displays the TF decomposition of the selected source.<<BR>><<BR>> {{attachment:selectionSource.gif}}
 * '''Right-click on the brain''': This selects the closest vertex and displays the popup menu at the same time. The first three menus are relative to the source that was just clicked. <<BR>><<BR>> {{attachment:popupSource.gif}}



== Hilbert transform ==
<<TAG(Advanced)>>

== Time and frequency bands ==
=== Frequency bands ===
 * Drag and drop the averaged file ''Left'' /'' ERF ''in  the Process tab, click on Run, select the process "Frequency >  Time-frequency (Morlet wavelets)". Click on "Edit" to get the options.
 * Just  change the frequency definition: select "Group in frequency bands",  leave the default frequency bands, and click on Ok. Click on Run.
 * Right-click on the file you've just computed "Power,FreqBands" > "One channel".<<BR>><<BR>> {{attachment:freqBands.gif}}
 * You  can notice that now the frequency selection is discreet, both in the   figure and in the frequency slider, it is not possible to select a  frequency in between. At each time and for each sensor, there are now  just 6 values, one per frequency band.
 * Now if you try to open  at the same time the first TF that you have computed "Power,1:1:60Hz",  you would get an error message: This frequency definition is not  compatible for display with the linear [1...60] scale.
 * Display  the other possible views (All sensors, 2D Sensor Cap, 2D Layout), and  change the current time and the current frequency in different ways  (sliders, keyboard, clicks on figures). Play with the open figures until  you are completely comfortable with these representations.

=== Time bands ===
 * Compute TF again for file ''Left / ERF'',  but this time select the option "Group in time bands" (with "linear"  frequency definition). To generate a list of regular time intervals,  click on "Generate" and enter "30" ms as the duration of each time band.
 * You get 11 time bands: t-2 and t-1 before 0ms, and t1 to t9 after 0ms. Click on Ok to start the computation.
 * Right-click on this new file "Power,1:1:60Hz" > Time-frequency maps.<<BR>><<BR>> {{attachment:timeBands.gif}}
 * Observe  that the behavior is not the same as with the frequency bands: it is  still possible to select a specific time in a time band (using the  slider or clicking on the figure). This is because the time definition  in the interface is still based on the ERF averaged file time  definition, which is "continuous". Only 11 values in time are computed,  and when requesting values a specific time point, the interface gets the  values associated with the time band that time point belongs to.
 * As  a consequence, it is possible to open at the same time the first TF  file that you have computed: "Power,1:1:60Hz". You would not get an  error message like with the frequencies mismatching.
 * Display  also other types of figures (All sensors, 2D topographies), change  current time and frequency several times in several ways (sliders,  keyboard, clicks on figures).

=== Time bands and frequency bands ===
 * Just for curiosity, compute a file with both time and frequency bands.<<BR>><<BR>> {{attachment:timeFreqTimeBand.gif}}
+The  brain is always active, the MEG/EEG recordings are never flat, some  oscillations are always present in the signals. Therefore we are often  more interested in the transient changes in the power at certain  frequencies than at the actual power. A good way to observe these changes  is to compute a '''deviation of the power with respect with a baseline'''.

There  is another reason for which we are usually interested in standardizing  the TF values. The power of the time-frequency coefficients are always  lower in the higher frequencies than in the lower frequencies, the  signal carries a lot less power in the faster oscillations  than in the slow brain responses. This '''1/f decrease in power''' is an observation that we already made with the power spectrum density in the [[Tutorials/ArtifactsFilter|filtering tutorial]]. If  we represent the TF maps with a linear color scale, we will always see  values close to zero in the higher frequency ranges. Normalizing each  frequency separately with respect with a baseline helps obtaining more  readable TF maps.

==== No normalization ====
The values we were looking at were already normalized (option: "1/f compensation" in the process options), but not with respect to a baseline. We will now compute the non-normalized power obtained with the Morlet wavelets and try the various options available for normalizing them.

 * In Process1, keep the simulated signal selected.
 * Run process: '''Frequency > Time-frequency (Morlet wavelets)''', "Spectral flattening: '''None'''"<<BR>><<BR>> {{attachment:example_process_none.gif||height="223",width="458"}}

 * Double-click on the file. As expected, we only see the lower frequencies in this representation: the power of the 2Hz oscillation is a lot larger than the power at 20Hz or 60Hz. <<BR>><<BR>> {{attachment:example_display_none.gif||height="131",width="433"}}

==== Spectral flattening ====
 * In Process1: Select this new non-normalized file "Power,1-60Hz".

 * Run process: '''Standardize > Spectral flattening''', Method='''1/f compensation'''. <<BR>><<BR>> {{attachment:norm_multiply.gif||height="224",width="454"}}

 * This produces exactly the same results as previously (option "1/f compensation" in the time-frequency process). It '''multiplies the power''' at each frequency bin '''with the frequency''' value (eg. multiplies by 20 the power at 20Hz), in order to correct for the 1/f shape we observe in the power spectrum. This works well for the lower part of the spectrum and up to 60-80Hz, but past this range it tends to '''overcompensate the higher frequencies'''.<<BR>>Note that it does not do any form of baseline correction: the 50Hz oscillation is visible everywhere.

==== Baseline normalization ====
 * The second way to proceed is to normalize the power with respect to its average level during a reference time period. We can consider the oscillations at 2Hz and 20Hz as our events of interest, and the 50Hz as noise we want to get rid of. The segment from '''0 to 2 seconds''' does not contain any of the signals of interest, therefore we can consider it as a '''baseline'''.

 * However, we will not be able to use the full segment [0,2]s because of the edge effects we described at the beginning of this tutorial. The time-frequency map at '''2Hz''' with the display option "Hide edge effects" (left figure below) shows that the power could not be estimated correctly before 0.75s, therefore we shouldn't use it as a baseline. The power time series at 2Hz (right) shows that the power related with the 2Hz oscillation starts to increase significantly after 1.25s, therefore it's not really a "baseline" anymore. This leaves only the segment '''[0.75,1.25]s''' available. <<BR>><<BR>> {{attachment:norm_baseline_2hz.gif||height="195",width="551"}}

 * At '''20Hz''', the expected transient effects are only 75ms long, therefore we could use a much longer baseline if we were not interested by the lower frequencies: '''[0.075, 1.925]s'''. <<BR>><<BR>> {{attachment:norm_baseline_20hz.gif||height="195",width="551"}}

 * In this case, we have a very long "resting" time segment (2s), therefore the edge effects are not a problem for picking a time window for the baseline normalization. We will use the first time window mentioned, [0.75,1.25]s as it is long enough (500ms) to estimate the baseline mean power. In real-life cases, with shorter epochs, it is sometimes difficult to find an acceptable '''trade-off between the baseline duration and the exclusion of the edge effects''', especially for the lower frequencies.

 * Run process: '''Standardize > Baseline normalization''', Baseline='''[0.75, 1.25]s'''.<<BR>>Method='''Event-related perturbation''': ERS/ERD stands for "event-related synchronization / desynchronization", a widely used normalization measure for time-frequency power maps. It evaluates the deviation from the mean over the baseline, in percents: '''(x-mean)/mean*100'''. <<BR>><<BR>> {{attachment:norm_ersd_process.gif||height="341",width="477"}}
 * Double-click on the new "ersd" file. The colormap type changed from "'''Timefreq'''" to "'''Stat2'''", which uses by default the "'''rwb'''" color set and shows relative values. Indeed, the ERS/ERD values can be '''positive or negative''': the power at a given time sample can be higher or lower than the average level during the baseline. In the simple simulation we used, there is no power decrease at any frequency after 2s, so the strong values are mostly positive. However, if you look in the file, you would see that there are many small negative values (due to the random noise we added). <<BR>><<BR>> {{attachment:norm_ersd_file.gif||height="187",width="499"}}
 * Note that the '''50Hz disappeared''' because it was present during the baseline, while the 2Hz and 20Hz oscillations show high positive values (between 100% and 2000% increase from baseline).
 * Remember to select your baseline very carefully according to the frequency range you are interested in. See below examples obtained with different baselines: '''0.075-1.925s''', '''0-2s''', '''0-6s'''. <<BR>>Change the colormap to "jet" if you prefer, and adjust the colormap contrast as needed. <<BR>><<BR>> {{attachment:norm_ersd_compare.gif||height="202",width="661"}}
 * This video may help you understand better what are the implications of the baseline selection: http://www.mikexcohen.com/lecturelets/whichbaseline/whichbaseline.html

<<TAG(Advanced)>>

== Tuning the wavelet parameters ==
==== Time resolution ====
You can adjust the relative time and frequency resolution of your wavelet transformation by adjusting the parameters of the mother wavelet in the options of the process.

 * Increasing the option "time resolution" will produce longer wavelets at a given frequency, hence increase the frequency accuracy (lower <<HTML(&Delta;)>>f) and decrease the time accuracy (higher <<HTML(&Delta;)>>t). Expect longer edge effects.
 * Decrease the time resolution will produce shorter wavelets at a given  frequency, hence decrease the frequency accuracy (higher <<HTML(&Delta;)>>f) and increase the time accuracy (higher <<HTML(&Delta;)>>t). Expect shorter edge effects.
 * You can modify one or the other parameter, what is important is the '''product of the two'''. All the following combinations fc/FWHM_t produce the same results because their product is constant: (1Hz,3s), (3Hz,1s), (6Hz,0.5s), (60Hz,0.05s) <<BR>><<BR>> {{attachment:tf_resolution_process.gif}}

Examples for a constant central frequency of '''1Hz''' with various time resolutions: '''1.5s''', '''4s''', '''10s'''.

 . {{attachment:tf_resolution_compare.gif||height="183",width="599"}}

==== Frequency axis ====
You can also obtain very different representations of the data by changing the '''list of frequencies''' for which you estimate the power. You can change this in the options of the process.

 . {{attachment:tf_freqlist_process.gif}}

Examples: Log 1:20:150, Log 1:300:150, Linear 15:0.1:25

 . {{attachment:tf_freqlist_compare.gif||height="183",width="598"}}

<<TAG(Advanced)>>

== Hilbert transform [TODO] ==
We can repeat the same analysis with the other approach available for exploring the simulated signal in the time-frequency plane. The process "'''Frequency > Hilbert transform'''" first filters the signals in various frequency bands with a band-pass filter, then computes the Hilbert transform of the filtered signal. The magnitude of the Hilbert transform of a narrow-band signal is a measure of the envelope of this signal, and therefore gives an indication of the activity in this frequency band.

 . {{attachment:slide_hilbert.gif}}

==== No normalization ====
Let's compute the same three results as before: non-normalized, 1/f compensation, baseline normalization.

 * In Process1, select the simulated signal.
 * Run process: '''Frequency > Hilbert transform''', Spectral flattening='''None''', '''Mirror signal''' before. <<BR>> The mirroring option helps with controlling the amplitude of the edge effects.<<BR>><<BR>> {{attachment:hilbert_process.gif||height="242",width="460"}}
 * In the advanced options panel, keep the default options: Default frequency bands and '''Power'''.<<BR>><<BR>>  {{attachment:hilbert_options.gif||height="353",width="362"}}
 * Double-click on the new file. The figure now has only 6 rows, one for each frequency band. <<BR>>Note that the values saved in the file and displayed are now '''amplitude''' and not power anymore. <<BR>><<BR>> {{attachment:hilbert_file.gif||height="197",width="606"}}
 * The non-normalized results are already easy to interpret:
  * '''delta (2-4Hz)''': Includes the 2Hz oscillation, contribution starts at 2s
  * '''beta (15-29Hz)''': Includes the 20Hz oscillation, contribution starts at 2s
  * '''gamma1(30-59Hz)''': Includes the 50Hz oscillation, contribution starts at the beginning (0s)
 * Right-click on the file or figure > Time series. Example for delta and beta.<<BR>><<BR>> {{attachment:hilbert_ts.gif||height="136",width="429"}}

==== Normalization ====
 * In Process1, select the non-normalized Hilbert-based decomposition.
 * Run process: '''Standardize > Spectral flattening''', Method='''1/f compensation'''.
 * Run process: '''Standardize > Baseline normalization''', Baseline='''[0.75, 1.25]s''', Method='''ERS/ERD'''.
 * Display the two normalized files side by side. <<BR>><<BR>> {{attachment:hilbert_normalized.gif||height="155",width="609"}}

==== Method specifications ====
 * Filters: '''[TODO]'''
 * Hilbert transformation: Using the [[http://www.mathworks.com/help/signal/ref/hilbert.html|Matlab's hilbert()]] function.
 * Extraction of the envelope: Power of the complex Hilbert transform, abs(hilbert(x))^2^.

== MEG recordings: Single trials ==
Let's go back to our auditory oddball paradigm and apply the concepts to MEG recordings. We will use all the trials available for one condition to estimate the average time-frequency decomposition.

==== Spectral flattening ====
 * In Process1, select all the '''deviant trials''' in '''Run#01'''.
 * Run process: '''Frequency > Time-frequency (Morlet wavelets)''', Spectral flattening: '''None'''. <<BR>>In the advanced options, select: '''Log 1:40:150''', '''Power''', Save '''average''' time-frequency maps. <<BR>><<BR>> {{attachment:tf_trials_process.gif||height="214",width="352"}} {{attachment:tf_trials_options.gif||height="318",width="300"}}
 * '''Save individual TF maps''': This   option stops the computation here and saves in the database one   time-frequency file for each input file (40 files), with one TF map for   each scout.
 * '''Save average TF maps''':   Instead of saving the TF for each file separately, it automatically   computes the average of the power of all the TF. This is a good choice   if you do not plan to use independently all the TF files, because it   saves a lot of time and disk space.
 * '''Remove evoked response from each trial before computing TF''': This option first computes the average of all the trials, then subtracts this average from each trial before computing the time-frequency decomposition. This brings the signals to a slightly more stationary state, which may help for the evaluation of the frequency contents.

==== Baseline normalization ====
 * Double-click on the new file '''Avg,Power,1-150Hz (MEG)'''. Select "'''Hide edge effects'''". <<BR>>In the drop-down list, select sensor '''MLP56''' (the one with the strongest respons at 90ms). <<BR>>Right-click on the TF figure > '''Time series'''. <<BR>><<BR>> {{attachment:tf_trials_evaluate.gif||height="163",width="559"}}
 * Defining a baseline is now a lot trickier than with 6s-long simulated signal. The epochs are only 600ms long, and the power at many frequencies could not be estimated correctly. If we want all the values to be "good" after 0s, '''we cannot use anything below 15Hz'''. If we want to normalize the values, we have to go even higher: '''30Hz''' '''if we want a baseline of 50ms''' before 0.
 * The epochs we use in this tutorial are '''too short to perform a correct time-frequency analysis'''. We should have imported at least 200ms more on each side, just for controlling the edge effects. You should always think carefully about the length of the epochs you import in your database if you are planning to run any form of frequency or time-frequency analysis.
 * For the purpose of illustrating the tools available in Brainstorm, we will keep on working with these short epochs. Let's try to do the best we can with what we have here. We could use a baseline of 50ms to get a correct estimation above 30Hz, but this is probably a bit too short. We propose to include some more of the baseline (75ms), hoping there are not major edge effects in this  segment.
 * In Process1, select the average time-frequency file.
 * Run process: '''Standardize > Baseline normalization''', Baseline='''[-75, 0]ms''', Method='''ERS/ERD'''. <<BR>><<BR>> {{attachment:tf_trials_normalize.gif||height="285",width="517"}}
 * The new menus available to display this file are described in the next section. <<BR>><<BR>> {{attachment:tf_trials_popup.gif||height="194",width="625"}}

==== Things to avoid ====
 * Avoid computing the time-frequency decomposition of the average of the trials, you would miss some of the '''induced response''', the brain activity in higher frequencies that is not strictly time-locked to the stimulus, and not aligned in phase across trials. Always prefer the computation of the average of the time-frequency power maps of each trials, as we did here. <<BR>>This is well documented in: [[http://www.iec-lnc.ens.fr/IMG/Files/Visual/Publications/2000_Review_IntJPsychol.pdf|Bertrand O, Tallon-Baudry C (2000)]].
 * Avoid using the "spectral flattening: '''1/f compensation'''" with frequencies '''above 60 or 80Hz'''. As mentioned before, it may overcompensate the higher frequencies. If you apply this normalization on the file we just computed, you get figures that give you a false idea of dominant high frequencies: <<BR>><<BR>> {{attachment:tf_trials_multiply.gif||height="129",width="216"}}

== Display: All channels ==
Three menus display the time-frequency of all the sensors with different spatial organizations. All the figures below represent the '''ERS/ERD'''-normalized average. They use the "jet" colormap, which is not the default configuration for these files. To get the same displays, change the colormap configuration: <<BR>>right-click on the figure > Colormap: Stat2 > Colormap > '''jet'''.

 * '''All channels''': All the maps are displayed one after the other, in the order they are saved in the file. <<BR>><<BR>> {{attachment:tf_allchannels.gif||height="295",width="345"}}
 * '''2D Layout (maps)''': Show each TF map where the sensor is located on a flattened 2D map. Most display options are available, such as the colormap management and the option "Hide edge effects".<<BR>><<BR>> {{attachment:tf_2dlayout_maps.gif||height="279",width="329"}} {{attachment:tf_2dlayout_hide.gif||height="279",width="328"}}
 * '''2D Layout (no overlap)''': Similar to the the previous display, but the positions of the images are reorganized so that they do not overlap. <<BR>><<BR>> {{attachment:tf_2dlayout_nooverlap.gif||height="294",width="345"}}

Useful shortcuts for these figures:

 * '''Click''': Clicking on any small TF image opens a new figure with only the selected sensor.
 * '''Shift + click''': Opens the original recordings time series of the selected sensor, when available. Here, we display an average of time-frequency maps, so this menu has no effect.
 * '''Mouse wheel''': Zoom in/out.
 * '''Right click + move''': Move in a zoomed figure.

== Display: Topography ==
The menus below show the distribution of TF power over the sensors, for one time point and one frequency bin, very similarly to what was introduced in tutorial [[Tutorials/ExploreRecordings|Visual exploration]].

 * '''2D Sensor  cap''' / '''2D Disc''' / '''3D Sensor cap''': 175ms, 8Hz<<BR>><<BR>> {{attachment:tf_topo.gif||height="137",width="474"}}
 * '''2D Layout''': 8Hz (black), 35Hz (white)<<BR>><<BR>> {{attachment:tf_2dlayout8.gif||height="241",width="260"}} {{attachment:tf_2dlayout35.gif||height="240",width="259"}}

Useful shortcuts for these figures:

 * '''Left/right arrows''': Change the current time.
 * '''Up/down arrows''':  Change the current frequency.
 * '''Control  + E''': Display the  sensors markers/names.

 * '''Shift  + click on a sensor''':  Displays the time-frequency decomposition for that specific  sensors.

 * '''Right click + move''': Select a group of sensors.

 * '''Shift +'''''' scroll''': Change the gain of the time series (2D Layout).
 * '''Control +'''''' scroll''': Change the length of the window displayed around the current time (2D Layout).

<<TAG(Advanced)>>

== Scouts ==
Similar calculations can be done at the level of the sources, either on the full cortex surface or on a limited number of regions of interests. We will start with the latter as it is usually an easier approach.

 * Drag and drop all the '''deviant trials''' from '''both runs''', select ['''Process sources'''].
 * Run process "'''Frequency > Time-frequency (Morlet wavelets)'''". <<BR>>Select the option "Use scouts" and '''select all the scouts''' defined in the previous tutorial. <<BR>><<BR>> {{attachment:scouts_process.gif||height="381",width="420"}}
 * In the advanced options, select "Scout function: '''After'''" and "Output: '''Save average'''". <<BR>>Run the process (it may take a while).<<BR>><<BR>> {{attachment:scouts_options.gif||height="598",width="361"}}
 * The '''scout function''' was introduced in the previous tutorial. It is the method we use to group the time series for the 20 dipoles we have in each scout into one unique signal. When computing the TF of one scout, we have the choice between applying this function before or after the time-frequency decomposition itself.
  * '''Before''': Extract the 20 source signals, apply the scout function to get one signal, run the TF decomposition of this signal. This is faster but may lead to information loss.
  * '''After''': Extract the 20 source signals, run the TF decomposition of the 20 signals, apply the scout function on the power of the TF maps. Always prefer this option when possible.
 * Rename the new file to add a tag "Deviant" in it. Then right-click > '''Time-freq: All scouts'''. <<BR>><<BR>> {{attachment:scouts_file.gif||height="193",width="505"}}
 * In Process1, select the new average TF file.
 * Run process: '''Standardize > Baseline normalization''', Baseline='''[-75, 0]ms''', Method='''ERS/ERD'''. <<BR>><<BR>> {{attachment:scouts_ersd.gif||height="186",width="562"}}

<<TAG(Advanced)>>

== Full cortical maps ==
Computing the time-frequency decomposition for all the sources of the cortex surface is possible but complicated because it can easily generate gigantic files, completely out of the reach of most computers. For instance the full TF matrix for each trial we have here would be [Nsources x Ntimes x Nfrequencies] = [15000 x 361 x 40] double-complex = '''3.2 Gb'''!

We have two ways of going around this issue: computing the TF decomposition for '''fewer frequency bins''' or frequency bands at a time, or as we did previously, use only limited number of '''regions of interest'''.

 * In Process1, keep all the '''deviant''' trials both conditions selected, select '''[Process sources]'''.
 * Run process "'''Frequency > Hilbert transform'''", Spectral flattening='''None''', '''Mirror signal''' before.<<BR>>To process the '''entire brain''', do '''not '''select the option "Use scouts".<<BR>><<BR>> {{attachment:sources_process.gif||height="398",width="432"}}
 * In the advanced options, select "Optimize storage: '''No'''", this option is not available when computing on the fly the average of multiple trials. Save the '''power''', Save the '''average '''Hilbert maps. <<BR>><<BR>> {{attachment:sources_options.gif||height="477",width="355"}}
 * '''Optimize the storage of the time frequency file''': Let's describe this option in more details.
  * When computing the TF decomposition of a source file, we are actually applying sequentially two linear transformations to the original recordings: the TF analysis and the source inversion. These two processes can be permuted: TF(Inverse(Recordings)) = Inverse(TF(Recordings)).
  * Therefore we can optimize the TF computation time by applying the wavelet transformation only to the sensor recordings, and then multiply the wavelet complex coefficients by the inverse operator (ImagingKernel). This  trick is always used in the computation of the Hilbert and Morlet  transforms.
  * When we have the option to the save the complex values (constrained sources and no averaging), this can also be used to optimize the storage  of the files. In these cases, we save only the wavelet transformation of the sensor data. Later, when the file is loaded for display, the imaging kernel is applied on the fly. This can be disabled explicitely with this option.

 * Rename the new Hilbert file to include the tag "Deviant", and select it in Process1.
 * Run process: '''Standardize > Baseline normalization''', Baseline=[-75, 0]ms, Method='''ERS/ERD'''.

 * Right-click on the Hilbert file > '''Display on cortex'''. <<BR>>The frequency slider now shows frequency bands ("alpha:8-12Hz") instead of frequencies ("12Hz"). You can explore the source activity in time and frequency dimensions. The screen capture below shows the activity at 175ms: a 60% increase in the alpha band around the auditory cortex and a 20% decrease in the beta band around the motor cortex.<<BR>><<BR>> {{attachment:sources_display.gif||height="161",width="540"}}

 * '''Shift + click '''on the cortex surface: Displays the TF decomposition of the selected source.
 * '''Right-click on the brain''': Selects the closest vertex and displays the popup menu at the same time. The first three menus are relative to the source that was just clicked. <<BR>><<BR>> {{attachment:sources_vertex.gif||height="163",width="480"}}

<<TAG(Advanced)>>

== Unconstrained sources ==
In the current example, we are working with the simple case: sources with constrained orientations. The unconstrained case is more difficult to deal with, because we have to handle correctly the three orientations we have at each vertex.

 * '''Full cortex''': Computes the TF decompositions for all the sources (3*15000=45000), then sum at each location the power for the three orientations.
 * '''Scouts''': Option "Scout function" in the process.
  * '''Before''': Extract the 20*3=60 source signals,  apply the scout function to get three signals (one per orientation), run  the TF decomposition of the three signals, and finally sum the power of  the three TF maps. This is faster but may lose some frequency  resolution (especially for constrained sources).
  * '''After''':  Extract the 20*3=60 source signals, run the TF decomposition of the 60  signals, apply the scout function on the power of the TF maps for each  orientation separately, and finally sum the power obtained for the three orientations.
 * The storage optimization option is not available with unconstrained sources.

<<TAG(Advanced)>>

== Getting rid of the edge effects ==
To avoid making mistakes in the manipulation of the data and producing more readable figures, we encourage you to cut out the edge effects from your time frequency maps after computation.

 * In Process1, select the very first computed in this tutorial: Test/Simulation/Power,1-60Hz | multiply
 * Run the process: "'''Extract > Extract time'''", Time window = '''[0.75, 5.25]s'''
 * Open the new file, select the option "'''Hide edge effects'''": Almost everything left in this new file is correctly estimated. Brainstorm keeps track of the edge effects in the TFmask field of the file. <<BR>><<BR>> {{attachment:tf_cut.gif||height="166",width="463"}}
 * We recommend you do the same when epoching your recordings: import trials that are longer than necessary, and after the time-frequency estimation, remove the unnecessary segments.
-Line 286:
+Line 350:
-Right  click one of the TF files, and select the menu File > View file  contents, to have a look to what is the actual contents of these  structures.

 . {{attachment:viewMatFile.gif}}

 * '''TF''':  [nSensors x nTime x nFrequency] matrix containing all the values of the  time-frequency decomposition (complex wavelet coefficients, or power of  amplitudes)
  * nSensors: Number of sensors for which the TF has been estimated
  * nTime:
   * If there are time bands defined, nTime = size(TimeBands,1) = number of time bands
   * If there are no time bands (linear time sampling): nTime = length(Time)
  * nFrequency:
   * If there are frequency bands defined, nFrequency = size(Freqs,1) = number of frequency bands
   * If there are no frequency bands (linear frequency sampling): nFrequency = length(Freqs)
+Right  click one of the first TF file we computed > File > '''View file  contents'''.

 . {{attachment:tf_contents.gif||height="450",width="563"}}

==== Structure of the time-frequency files: timefreq_*.mat ====
 * '''TF''':  [Nsignals x Ntime x Nfreq] matrix containing all the values of the  time-frequency decomposition (complex wavelet coefficients, or double values for power/magnitude/Z-score).
 * '''TFmask''': [Nfreq x Ntime] logical mask indicating the edge effects (0=edge, 1=valid value).
 * '''Std''': [Nsignals x Ntime x Nfreq] standard deviation if this file is an average.
-Line 299:
+Line 359:
- * '''DataType''': Explains from what kind of data this file was computed. Possible values are: 'data', 'results', 'scout, 'matrix'...
 * '''Time''': Time vector used to estimate this file.
 * '''TimeBands''':  Description of the time bands. Cell array {nTimeBands x 3}, where each  line represents a time band {'band_name', 'time definition', 'function'}
 * '''Freqs''':
  * If frequency  bands: Cell array {nFreqBands x 3}, where each line represents a  frequency band  {'band_name', 'frequency definition', 'function'}
  * In linear frequency sampling: vector containing all the frequencies
 * '''RowNames''':  Cell array of strings that describes each row of the TF matrix. In this  specific case, it would be the list of all the MEG sensor names. But it  could also be a list of names of scouts or clusters.
+ * '''DataType''': From what kind of data this file was computed: {'data', 'results', 'scout, 'matrix'}
 * '''Time''': [1 x Ntime] Time vector used to estimate this file.
 * '''TimeBands''':  [Ntimebands x 3] Cell array where each  line represents a time band: <<BR>>{'band_name', 'time definition', 'function'}
 * '''Freqs''': For regular frequency binning: vector containing all the frequencies.<<BR>>If using frequency  bands: [Nfreqbands x 3] cell array, where each line represents a  frequency band  {'band_name', 'frequency definition', 'function'}

 * '''RefRowNames''': Used only for connectivity matrices.
 * '''RowNames''': [1 x Nsignals]  Cell array of strings that describes each row of the TF matrix. In this  specific case, it would be the list of all the MEG sensor names. But it  could also be a list of names of scouts or clusters.
-Line 307:
+Line 367:
-  * 'none': No measure applied, TF contains the complex wavelet coefficients.
  * 'power': TF contains the power for each frequency, ie. the square of the amplitude: abs(coefficients).^2
  * 'magnitude':  abs(coefficients)
  * 'log': 10 * log10( abs(coefficients).^2)
  * 'phase': angle(coefficients)
 * '''Method''':  String that identifies the process that generated the file; possible  values: {'morlet', 'fft', 'psd', 'hilbert', 'corr', 'cohere', 'granger',  'plv'}
+  * '''none''': No measure applied, TF contains the complex wavelet coefficients.
  * '''power''': Power for each frequency, ie. the square of the amplitude: abs(coefficients)^2^
  * '''magnitude''':  abs(coefficients)
  * '''log''': 10 * log10( abs(coefficients)^2^)
  * '''phase''': angle(coefficients)
 * '''Method''':  String that identifies the process that generated the file: <<BR>>{'morlet', 'fft', 'psd', 'hilbert', 'corr', 'cohere', 'granger',  'plv'}
-Line 314:
+Line 374:
- * '''nAvg''': Number of trials that were averaged to obtain this file (copied from the recordings file at computation)
+ * '''SurfaceFile''' / '''GridLoc''' / '''GridAtlas''': Source space that was used, only for source files.
 * '''nAvg''': Number of trials that were averaged to obtain this file.
 * '''ColormapType''': String, force a specific colormap type to be used when displaying this file.
 * '''DisplayUnits''': String, force to use specific units when displaying this file.
-Line 316:
+Line 379:
- * '''History''': List of operations performed on this file (better visualized with popup menu File > View file history)

'''Document file tags '''

== Additional discussions on the forum ==
 * Time and frequency resolution: http://neuroimage.usc.edu/forums/showthread.php?1848
+ * '''History''': List of operations performed on this file (menu File > View file history).

==== Useful functions ====
 * '''in_bst_timefreq'''(TimefreqFile): Read a time-frequency file.
 * '''in_bst'''(FileName, TimeWindow): Read any Brainstorm file with the possibility to load  only a specific part of the file. "TimeWindow" is an range of time  values in seconds: [tStart, tStop].
 * '''bst_process'''(''''LoadInputFile'''',   FileName, Target, TimeWindow): The most high-level function  for  reading data files. "Target" is a string with the list of sensor names  or types to load.
 * '''morlet_transform'''(): Applies complex Morlet wavelet transform to the time series in input.

<<TAG(Advanced)>>

== References ==
 * Bertrand O, Tallon-Baudry C (2000)<<BR>> [[http://www.iec-lnc.ens.fr/IMG/Files/Visual/Publications/2000_Review_IntJPsychol.pdf|Oscillatory gamma activity in humans: a possible role for object representation]]<<BR>>Int J Psychophysiol, 38(3):211-23
 * Bruns A (2004)<<BR>>[[http://www.sciencedirect.com/science/article/pii/S0165027004001098|Fourier-, Hilbert- and wavelet-based signal analysis: are they really different approaches?]]<<BR>>J Neurosci Methods, 137(2):321-32
 * Le Van Quyen M, Foucher J, Lachaux J, Rodriguez E, Lutz A, Martinerie J, Varela FJ (2001) <<BR>>[[http://brainimaging.waisman.wisc.edu/~lutz/LeVanQuyen_et_all_JNM_2001.pdf|Comparison of Hilbert transform and wavelet methods for the analysis of neuronal synchrony]]<<BR>>J Neurosci Methods, 111(2):83-98
 * Pfurtscheller G (1992)<<BR>>[[http://www.sciencedirect.com/science/article/pii/0013469492901333|Event-related synchronization (ERS): an electrophysiological correlate of cortical areas at rest]]<<BR>>Electroencephalogr Clin Neurophysiol, 83(1):62-9
 * Series of lectures by Mike X Cohen: <<BR>>http://www.mikexcohen.com/lectures.html

<<TAG(Advanced)>>

== Additional documentation ==
 * Forum: Time and frequency resolution: http://neuroimage.usc.edu/forums/showthread.php?1848
-Line 325:
+Line 403:
-<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Scouts&next=Tutorials/Statistics")>>
+<<EmbedContent("http://neuroimage.usc.edu/bst/get_prevnext.php?prev=Tutorials/Scouts&next=Tutorials/Difference")>>

Feedback: Comments, bug reports, suggestions, questions

Email address (if you expect an answer):