Time-frequency analysis using wavelets, Hanning windows and multitapers

Contents

Introduction

In this tutorial we will demonstrate time-frequency analysis of a single subject's MEG data using a Hanning window, multitapers and wavelets. This tutorial also shows how to visualize time-frequency results (TFRs).

We assume that the reader already knows how to do basic preprocessing in FieldTrip. Since there are multiple events (left/right cue, fixation, stimulus onset, speed change, left/right response), the sequence of triggers in the dataset is considerably more complex than in the somatosensory experiment on the same subject. We will show how to parse the sequence of triggers, select the data segments of interest and preprocess the data.

There is no information in this tutorial about how to compare conditions, how to average over subjects, how to analyze at the source level or how to do statistical analysis on the time-frequency data. Some of these issues are covered in other tutorials (see Summary and suggested further reading).

Background

Oscillatory components contained in the ongoing EEG or MEG signal often show power changes relative to experimental events. These signals are not necessarily phase-locked to the event and will not be represented in event related fields and potentials. The goal of this analysis is to compute and visualize event related changes in power by calculating time-frequency representations (TFRs). This will be done using analysis based on Fourier analysis and wavelets. The Fourier analysis will include the application of multitapers, which allow a better control of time and frequency smoothing.

Calculating time-frequency representations of power is done using a sliding time window. This can be done according to two principles: either the time window has a fixed length independent of frequency, or the time window decreases in length with increased frequency. For each time window the power is calculated. Prior to calculating the power one or more tapers are multiplied with the data. The aim of the tapers is to reduce spectral leakage and control the frequency smoothing.

We will work on a dataset of an visual change-detection experiment. Magneto-encephalography (MEG) data was collected using the 306-channel Neuromag Triux MEG system at NatMEG. The subject was instructed to fixate on the screen and press a left- or right-hand button upon a speed change of an inward moving stimulus Each trial started with the presentation of a cue pointing either rightward or leftward, indicating the response hand. Subsequently a fixation cross appeared. After a baseline interval of 1s, an inward moving circular sine-wave grating was presented. After an unpredictable delay, the grating changed speed, after which the subjects had to press a button with the cued hand. In a small fraction of the trials there was no speed change and subjects should not respond.

The dataset allows for a number of different questions, e.g.

  1. What is the effect of cueueing?
  2. What is the effect of stimulus onset?
  3. What is the effect of the ongoing visual stimulation and the maintained attention?
  4. What is the effect of the response target, i.e. the speed change?

These questions require different analysis strategies on different sections of the data. Hence we will spend some time on segmenting and preparing the data for the different analyses, although we will focus on question 3.

Procedure

To calculate the time-frequency analysis for the example dataset we will perform the following steps:

In this tutorial four strategies for time-frequency analysis will be demonstrated, numbered 1 to 4.

Defining data segments of interest

The first step is to read the data using the function ft_preprocessing, which requires the specification of the trials by ft_definetrial.

Let us start by looking at the sequence of triggers that is present in the data:

clear all
close all

filename = 'workshop_visual_sss.fif';

hdr   = ft_read_header(filename);
event = ft_read_event(filename);

sample = [event.sample];
value  = [event.value];

figure
plot(sample/hdr.Fs, value, '.');
xlabel('time (s)')
	306 MEG channel locations transformed
Opening raw data file workshop_visual_sss.fif...
	Range : 7000 ... 994999  =      7.000 ...   994.999 secs
Ready.
	306 MEG channel locations transformed
Opening raw data file workshop_visual_sss.fif...
	Range : 7000 ... 994999  =      7.000 ...   994.999 secs
Ready.
Reading 7000 ... 994999  =      7.000 ...   994.999 secs... [done]

Along the horizontal axis is time in the experiment, along the vertical axis is the trigger value. Each dot represents one trigger.

There are in total 14 different triggers present in the data:

unique([event.value])
ans =

  Columns 1 through 13

    10    11    12    13    14    16    17    20    21    22    23    24    26

  Column 14

    27

You can zoom in in the figure to show the trigger sequence of a few trials.

axis([0.1168   41.4003    9.7295   17.3887])

Looking at the triggers, you can recognize a certain diagonal 'stripe' structure that repeats itself. Fully deciphering these triggers requires more information from the experimenter or from the stimulus presentation script (which is shared along with the data). The Presentation PCL file is an ascii file with the code that controls the stimulation computer and sends out triggers to the MEG system. You can open it in an editor. It contains the following description of the triggers, whose values can be recognized in the MEG dtaa.

int code_cue_left       = 10;
int code_cue_right      = 20;
int code_bason_L        = 11;   #LEFT baseline, onset of baseline stimulus
int code_stimon_L       = 12;   #LEFT stimulus, onset of stimulus
int code_speedup_L      = 13;   #LEFT stimulus, speed change of stimulus
int code_rtype_rok_L    = 14;   #LEFT response-type, correct response, speed detected.
int code_rtype_nrok_L   = 15;   #LEFT response-type, correct response, there was no speedchange.
int code_rtype_early_L  = 16;   #LEFT response-type, too early.
int code_rtype_late_L   = 17;   #LEFT response-type, too late response, as a matter of fact: no response.
int code_bason_R        = 21;   #RIGHT baseline, onset of baseline stimulus
int code_stimon_R       = 22;   #RIGHT stimulus, onset of stimulus
int code_speedup_R      = 23;   #RIGHT stimulus, speed change of stimulus
int code_rtype_rok_R    = 24;   #RIGHT response-type, correct response, speed detected.
int code_rtype_nrok_R   = 25;   #RIGHT response-type, correct response, there was no speedchange.
int code_rtype_early_R  = 26;   #RIGHT response-type, too early.
int code_rtype_late_R   = 27;   #RIGHT response-type, too late response, as a matter of fact: no response.

You can see that not the only responses themselves are coded, but also whether they are correct. The interesting segments of data are not simply a piece of data around a single trigger, but around a specific sequence. E.g. the trigger sequence [10 11 12 13 14] codes for a complete trial in which a correct response was given with the left hand.

The ft_definetrial function, which you were introduced to in the preprocessing hands-on, has built-in support for the segmentation on the basis of single triggers, not for complex trigger sequences. To support complex trigger sequences, ft_definetrial allows you to provide your own piece of MATLAB code. That piece of code should be written as a function. In FieldTrip we refer to such a function as a 'trial function' or trialfun.

INSTRUCTION: Open the provided trial function

edit trialfun_visual_cue

This trialfun_visual_cue function is defined such that it takes the cfg structure as input, and produces a 'trl' matrix as output. Furthermore, it returns the event structure.

function [trl, event] = trialfun_visual_cue(cfg)

The code of this function makes use of the cfg.dataset field and of cfg.trialdef.pre = time before the first trigger cfg.trialdef.post = time after the last trigger

which you will have to specify. At line 47 it starts looping over all cues, and for each checks whether the subsequent trigger sequence corresponds to a correctly executed left- or right-hand trial.

On line 61-63 the most important characteristics of the 'trl' matrix are specified: the begin and endsample, and the offset. These code where you want the data segment to start and end, and which sample of the data segment you want to consider as time t=0, i.e. as the time at which you want to align all trials in the analysis.

On line 65-73 some additional characteristics of each trial are determined, which are all concatenated to the trl matrix on line 76. The first three columns of 'trl' are required, all other columns are optional and can for example be used to perform analysis on subsets or for statistical analyses.

Also provided are the trialfun for selecting a data segment around the response (trialfun_visual_response) and around the stimulus (trialfun_visual_stimulus). Since we are for now interested in the visual response, we will continue with trialfun_visual_stimulus.

edit trialfun_visual_stimulus

You pass the name of your trial function (or a function handle to it) to ft_definetrial:

cfg = [];
cfg.dataset   = filename;
cfg.trialfun  = 'trialfun_visual_stimulus';
cfg = ft_definetrial(cfg);
evaluating trialfunction 'trialfun_visual_stimulus'
	306 MEG channel locations transformed
Opening raw data file workshop_visual_sss.fif...
	Range : 7000 ... 994999  =      7.000 ...   994.999 secs
Ready.
	306 MEG channel locations transformed
Opening raw data file workshop_visual_sss.fif...
	Range : 7000 ... 994999  =      7.000 ...   994.999 secs
Ready.
Reading 7000 ... 994999  =      7.000 ...   994.999 secs... [done]
found 1473 events
created 125 trials
the call to "ft_definetrial" took 14 seconds and required the additional allocation of an estimated 0 MB

The ft_definetrial function will do some general stuff and subsequently call your function, passing the cfg along. If you compare the cfg structure before and after the call to ft_definetrial, you see that cfg.trl has been added

disp(cfg)
           dataset: 'workshop_visual_sss.fif'
          trialfun: @trialfun_visual_stimulus
          callinfo: [1x1 struct]
           version: [1x1 struct]
       trackconfig: 'off'
       checkconfig: 'loose'
         checksize: 100000
      showcallinfo: 'yes'
             debug: 'no'
     trackcallinfo: 'yes'
     trackdatainfo: 'no'
    trackparaminfo: 'no'
          datafile: 'workshop_visual_sss.fif'
        headerfile: 'workshop_visual_sss.fif'
        dataformat: 'neuromag_fif'
...
disp(cfg.trl)
   1.0e+05 *

    0.0514    0.0720   -0.0075    0.0000    0.0000
    0.0872    0.1039   -0.0075    0.0000    0.0000
    0.1188    0.1326   -0.0075    0.0000    0.0000
    0.2041    0.2173   -0.0075    0.0000    0.0000
    0.2324    0.2472   -0.0075    0.0000    0.0000
    0.2898    0.3055   -0.0075    0.0000    0.0000
    0.3203    0.3367   -0.0075    0.0000    0.0000
    0.4082    0.4221   -0.0075    0.0000    0.0000
    0.4371    0.4535   -0.0075    0.0000    0.0000
    0.4985    0.5133   -0.0075    0.0000    0.0000
    0.5284    0.5476   -0.0075    0.0000    0.0000
    0.5626    0.5803   -0.0075    0.0000    0.0000
    0.5952    0.6092   -0.0075    0.0000    0.0000
...

If you were to look at the code of ft_definetrial, you would see that by default it uses the (separate) function ft_trialfun_general.

which ft_trialfun_general
/Users/robert/matlab/fieldtrip/trialfun/ft_trialfun_general.m
help ft_trialfun_general
  FT_TRIALFUN_GENERAL determines trials/segments in the data that are
  interesting for analysis, using the general event structure returned
  by read_event. This function is independent of the dataformat
 
  The trialdef structure can contain the following specifications
    cfg.trialdef.eventtype  = 'string'
    cfg.trialdef.eventvalue = number, string or list with numbers or strings
    cfg.trialdef.prestim    = latency in seconds (optional)
    cfg.trialdef.poststim   = latency in seconds (optional)
 
  If you want to read all data from a continous file in segments, you can specify
     cfg.trialdef.triallength = duration in seconds (can be Inf)
     cfg.trialdef.ntrials     = number of trials
 
  If you specify
...

There are some other example trial functions in the fieldtrip/trialfun directory, such as ft_trialfun_example1, which detects the occurence of two subsequent triggers, and ft_trialfun_example2, which detects the onset of muscle activity in an EMG channel. The FieldTrip mechanism of "trialfuns" allows you to segment your data, regardless of the complexity, using straight-forward MATLAB code. More information can be found in the online documentation.

Reading and preprocessing the raw data

Having defined the segments of interest, we can specify additional preprocessing options in the cfg-structure and call ft_preprocessing to read and baseline-correct the data.

cfg.channel         = 'MEGGRAD';
cfg.baselinewindow  = [-inf 0];
cfg.demean          = 'yes';
data = ft_preprocessing(cfg);
	306 MEG channel locations transformed
Opening raw data file workshop_visual_sss.fif...
	Range : 7000 ... 994999  =      7.000 ...   994.999 secs
Ready.
processing channel { 'MEG0112' 'MEG0113' 'MEG0122' 'MEG0123' 'MEG0132' 'MEG0133' 'MEG0142' 'MEG0143' 'MEG0212' 'MEG0213' 'MEG0222' 'MEG0223' 'MEG0232' 'MEG0233' 'MEG0242' 'MEG0243' 'MEG0312' 'MEG0313' 'MEG0322' 'MEG0323' 'MEG0332' 'MEG0333' 'MEG0342' 'MEG0343' 'MEG0412' 'MEG0413' 'MEG0422' 'MEG0423' 'MEG0432' 'MEG0433' 'MEG0442' 'MEG0443' 'MEG0512' 'MEG0513' 'MEG0522' 'MEG0523' 'MEG0532' 'MEG0533' 'MEG0542' 'MEG0543' 'MEG0612' 'MEG0613' 'MEG0622' 'MEG0623' 'MEG0632' 'MEG0633' 'MEG0642' 'MEG0643' 'MEG0712' 'MEG0713' 'MEG0722' 'MEG0723' 'MEG0732' 'MEG0733' 'MEG0742' 'MEG0743' 'MEG0812' 'MEG0813' 'MEG0822' 'MEG0823' 'MEG0912' 'MEG0913' 'MEG0922' 'MEG0923' 'MEG0932' 'MEG0933' 'MEG0942' 'MEG0943' 'MEG1012' 'MEG1013' 'MEG1022' 'MEG1023' 'MEG1032' 'MEG1033' 'MEG1042' 'MEG1043' 'MEG1112' 'MEG1113' 'MEG1122' 'MEG1123' 'MEG1132' 'MEG1133' 'MEG1142' 'MEG1143' 'MEG1212' 'MEG1213' 'MEG1222' 'MEG1223' 'MEG1232' 'MEG1233' 'MEG1242' 'MEG1243' 'MEG1312' 'MEG1313' 'MEG1322' 'MEG1323' 'MEG1332' 'MEG1333' 'MEG1342' 'MEG1343' 'MEG1412' 'MEG1413' 'MEG1422' 'MEG1423' 'MEG1432' 'MEG1433' 'MEG1442' 'MEG1443' 'MEG1512' 'MEG1513' 'MEG1522' 'MEG1523' 'MEG1532' 'MEG1533' 'MEG1542' 'MEG1543' 'MEG1612' 'MEG1613' 'MEG1622' 'MEG1623' 'MEG1632' 'MEG1633' 'MEG1642' 'MEG1643' 'MEG1712' 'MEG1713' 'MEG1722' 'MEG1723' 'MEG1732' 'MEG1733' 'MEG1742' 'MEG1743' 'MEG1812' 'MEG1813' 'MEG1822' 'MEG1823' 'MEG1832' 'MEG1833' 'MEG1842' 'MEG1843' 'MEG1912' 'MEG1913' 'MEG1922' 'MEG1923' 'MEG1932' 'MEG1933' 'MEG1942' 'MEG1943' 'MEG2012' 'MEG2013' 'MEG2022' 'MEG2023' 'MEG2032' 'MEG2033' 'MEG2042' 'MEG2043' 'MEG2112' 'MEG2113' 'MEG2122' 'MEG2123' 'MEG2132' 'MEG2133' 'MEG2142' 'MEG2143' 'MEG2212' 'MEG2213' 'MEG2222' 'MEG2223' 'MEG2232' 'MEG2233' 'MEG2242' 'MEG2243' 'MEG2312' 'MEG2313' 'MEG2322' 'MEG2323' 'MEG2332' 'MEG2333' 'MEG2342' 'MEG2343' 'MEG2412' 'MEG2413' 'MEG2422' 'MEG2423' 'MEG2432' 'MEG2433' 'MEG2442' 'MEG2443' 'MEG2512' 'MEG2513' 'MEG2522' 'MEG2523' 'MEG2532' 'MEG2533' 'MEG2542' 'MEG2543' 'MEG2612' 'MEG2613' 'MEG2622' 'MEG2623' 'MEG2632' 'MEG2633' 'MEG2642' 'MEG2643' }
reading and preprocessing
Reading 12142 ... 14201  =     12.142 ...    14.201 secs... [done]
reading and preprocessing trial 2 from 125
Reading 15719 ... 17386 reading and preprocessing trial 3 from 125
Reading 18878 ... 20262 reading and preprocessing trial 4 from 125
Reading 27408 ... 28725  =     27.408 ...    28.725 secs... [done]
Reading 30242 ... 31718 reading and preprocessing trial 6 from 125
Reading 35978 ... 37546 reading and preprocessing trial 7 from 125
Reading 39030 ... 40664  =     39.030 ...    40.664 secs... [done]
Reading 47817 ... 49210 reading and preprocessing trial 9 from 125
...

The raw data is represented in a MATLAB structure and you can use standard MATLAB code to access and visualize it:

disp(data);

figure
plot(data.time{1}, data.trial{1});
           hdr: [1x1 struct]
         label: {204x1 cell}
          time: {1x125 cell}
         trial: {1x125 cell}
       fsample: 1000
    sampleinfo: [125x2 double]
     trialinfo: [125x2 double]
          grad: [1x1 struct]
           cfg: [1x1 struct]

In the data.trialinfo field you can recognize the additional columns that were added to the 'trl' matrix by trialfun_visual_stimulus. They represent the cued hand (1=left, 2=right) and the response latency, i.e. time between the speed change and the response.

We can have a closer look at the behavioral data, i.e, the response time.

figure
plot(data.trialinfo(:,2), '.');
xlabel('trial number')
ylabel('response latency (s)')

% There is an unexpected structure in this scatterplot, whish shows up as
% stripes. Sorting the trials according to response latency reveals this:

figure
plot(sort(data.trialinfo(:,2)), '.');
xlabel('sorted trial number')
ylabel('response latency (s)')

You can now clearly see that the response latencies are quantized at approximately 8 ms, plus or minus one sample. This is a consequence of the processing of the responses in Presentation, which synchronizes with the screen refresh rate of 120 Hz. It means that the behavioral data has not been recorded with the highest possible accuracy (1 ms, since the sampling frequency is 1000 Hz).

We continue processing the data. Using ft_rejectvisual we can perform a quick scan of the overall data quality and exclude noisy trials:

cfg             = [];
cfg.method      = 'summary';
cfg.channel     = 'MEG';
cfg.keepchannel = 'yes';
dataclean = ft_rejectvisual(cfg, data);
the input is raw data with 204 channels and 125 trials
Warning: correcting numerical inaccuracy in the time axes 
showing a summary of the data for all channels and trials
computing metric [--------------------------------------------------------/]
125 trials marked as GOOD, 0 trials marked as BAD
204 channels marked as GOOD, 0 channels marked as BAD
no trials were removed
the call to "ft_rejectvisual" took 8855 seconds and required the additional allocation of an estimated 12 MB

Let us save the cleaned data, to allow resuming the analysis at a later point in time.

save visual_data.mat      data
save visual_dataclean.mat dataclean

clear data % this is not needed any more

Inspect the ERFs

Although we are not interested in the ERFs in this tuturial, it is good practice to check at various stages in your analysis whether the expected features in the data are present. Here we know that a strong visual stimulus switches on at t=0, which means that there should be a Visual Evoked Field (VEF).

load visual_dataclean.mat

cfg = [];
cfg.vartrllength  = 2;  % the default is NOT to allow variable length trials
timelock_planar   = ft_timelockanalysis(cfg, dataclean);

figure
plot(timelock_planar.time, timelock_planar.avg);
grid on
the input is raw data with 204 channels and 125 trials
Warning: correcting numerical inaccuracy in the time axes 
averaging trials
averaging trial 125 of 125

the call to "ft_timelockanalysis" took 3 seconds and required the additional allocation of an estimated 19 MB

The figure above shows the VEF. You should note that towards the end, from about 1 second post-stimulus, the ERF traces start to diverge. The very last piece of the data at t=1.3 is really far off from the rest.

QUESTION: Can you explain the data getting more noisy towards the end?

Since the speed change happened at an unpredictable point in time, the subject's response can be anywhere from 0.5 seconds onwards. We did not want to include the resonses in the analysis, hence we only selected (in trialfun_visual_stimulus) the data segments up to the speed change, excluding the response. We can visualize this for a few trials with

channel = 131;
scale = 4e-11;
for k = 1:10
  plot(dataclean.time{k}, dataclean.trial{k}(1,:)+k*scale);
  hold on;
end
plot([0 0], [0 1], 'k');
ylim([0 11*scale]);
set(gca, 'ytick', (1:10).*scale);
set(gca, 'yticklabel', 1:10);
ylabel('trial number');
xlabel('time (s)');
title(dataclean.label{channel});

We can look at the full distribution of the trial durations using some standard MATLAB code:

figure
for i=1:length(dataclean.time)
  starttime = dataclean.time{i}(1);
  endtime   = dataclean.time{i}(end);
  x = [starttime endtime];
  y = [i i];
  line(x, y);
end

You see that the data segments end due to the speed-change anywhere from 0.5 seconds post-stimulus up to 1.3 seconds post-stimulus. Up to 0.5 seconds there are many trials that can be averaged. Towards the end, the number of trials that contributes to the average decreases to one, causing the last section of the average to be very noisy.

The number of trials that went into the average at every sample is also represented in the average timelock structure:

figure
plot(timelock_planar.time, timelock_planar.dof(1,:)); % it is the same for all channels
grid on

Let us focus on the data again. Having computed the ERF, we can look at its spatial distribution:

cfg = [];
cfg.layout = 'neuromag306planar.lay';
figure
ft_multiplotER(cfg, timelock_planar);
selection avg along dimension 1
selection dof along dimension 1
selection var along dimension 1
reading layout from file neuromag306planar.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the call to "ft_multiplotER" took 4 seconds and required the additional allocation of an estimated 13 MB

For a correct interpretation of the planar gradient topography, you should combine the two planar gradients at each location into a single (absolute) magnitude:

cfg = [];
timelock_cmb = ft_combineplanar(cfg, timelock_planar);

cfg = [];
cfg.layout = 'neuromag306cmb.lay';
ft_multiplotER(cfg, timelock_cmb);
the input is timelock data with 204 channels and 2094 timebins
the input is timelock data with 204 channels and 2094 timebins
Warning: the trial definition in the configuration is inconsistent with the
actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive
segments of a continuous recording 
the input is raw data with 102 channels and 1 trials
the call to "ft_combineplanar" took 1 seconds and required the additional allocation of an estimated 0 MB
selection avg along dimension 1
reading layout from file neuromag306cmb.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the call to "ft_multiplotER" took 1 seconds and required the additional allocation of an estimated 0 MB

INSTRUCTION: Use the interactive mode of the figure to zoom in on the occipital channels. Can you see the occipital topography of the VEF peak around 100ms?

If you look at the single timecourse (i.e. the average ERF of the occipital channels), you see that the planar magnitude increases towards the end. You also see that the baseline is not zero. This is a consequence of the non-linear operation of taking the planar-gradient magnitude. The non-linear sensitivity of the planar gradient to the noise has consequences for statistical comparisons, for example in a group study: in conditions with fewer trials, you always expect more noise. After taking the magnitute (absolute value) this noise is always positive and does not average out over subjects. So when comparing a condition with on average 50 trials with a condition with 150 trials, you expect the 150-trial condition to be less biassed in each of the subjects, and hence have a smaller magnitude of the planar gradient. A statistical comparison of unbalanced data that has been processed like this is not meaningful, since you already expect the difference to be visible in the baseline of your exeriment.

INSTRUCTION: Explore the activity at the start of the window of interest. What seems to be going on in the latency range from -0.8 to -0.4 seconds?

Time-frequency analysis

To speed up the subsequent computations and reduce memory requirements, we will downsample the data from 1000 to 250 Hz. This is not recommended for normal analysis, where you would better wait a bit longer. However, in this tutorial we want to compare a number of analyses without having the time to wait.

load visual_dataclean

cfg = [];
cfg.resamplefs  = 250;
cfg.demean      = 'no';
cfg.detrend     = 'no';
resampled = ft_resampledata(cfg, dataclean);

clear dataclean
the input is raw data with 204 channels and 125 trials
Warning: correcting numerical inaccuracy in the time axes 
resampling data
Warning: not all of the trials have the same original time axis: to avoid
rounding issues in the resampled time axes, data will be zero-padded to the left
prior to resampling 
resampling data in trial 125 from 125

original sampling rate = 1000 Hz
new sampling rate = 250 Hz
the call to "ft_resampledata" took 4 seconds and required the additional allocation of an estimated 84 MB

We use the ft_freqanalysis function for the time-frequency analysis, which takes a raw data structure as input. It implements various algorithms, which can be specified using cfg.method. Each of the methods has additional configuration options.

INSTRUCTION: look at the help of ft_freqanalysis and see whether you can find the options for the 'wavelet' method.

(1) Time-frequency analysis with Hanning taper, fixed window length

We will describe how to calculate time frequency representations using Hanning tapers. When choosing for a fixed window length, the frequency resolution is defined according to the length of the time window (delta T). The frequency resolution (delta F in figure 1) is 1/length of time window in sec (delta T in figure 1). Thus a 500 ms time window results in a 2 Hz frequency resolution (1/0.5 sec = 2 Hz), meaning that power can be calculated for 2 Hz, 4 Hz, 6 Hz etc. An integer number of cycles must fit in the time window.

cfg              = [];
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.output       = 'pow';
cfg.pad          = 2.5;
cfg.foi          = 2:4:120;                        % estimate power at 2 to 120 Hz in steps of 4 Hz
cfg.t_ftimwin    = ones(length(cfg.foi),1).*0.5;   % length of time window = 0.5 sec
cfg.toi          = -0.750:0.05:1.350;              % time window "slides" in steps of 50 ms
TFRhann = ft_freqanalysis(cfg, resampled);
the input is raw data with 204 channels and 125 trials
Warning: the trial definition in the configuration is inconsistent with the
actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive
segments of a continuous recording 
Warning: correcting numerical inaccuracy in the time axes 
Warning: correcting numerical inaccuracy in the time axes 
processing trials
Warning: output time-bins are different from input time-bins 
trial 125, frequency 30 (118.00 Hz), 1 tapers

the call to "ft_freqanalysis" took 55 seconds and required the additional allocation of an estimated 0 MB

Since this step takes quite some time, you might want to save the data afterwards

save visual_TFRhann.mat TFRhann

Regardless of the method used for calculating time-frequency power, the output format is an identical MATLAB structure with the following elements:

disp(TFRhann)
        label: {204x1 cell}
       dimord: 'chan_freq_time'
         freq: [1x30 double]
         time: [1x43 double]
    powspctrm: [204x30x43 double]
         grad: [1x1 struct]
          cfg: [1x1 struct]

The element TFRhann.powspctrm contains the temporal evolution of the raw power values for each specified frequency. The element TFRhann.dimord describes the three dimensions: channels by frequency by time. Each of the dimensions is additionally described by TFRhann.label, TFRhann.freq and TFRhann.time.

To visualize the data, we will combine the two planar channels. Whereas for ERFs it requires application of Pythagoras' rule, in the case of power it simply consists of adding the power of the two respective planar channels.

cfg = [];
TFRhann_cmb = ft_combineplanar(cfg, TFRhann);
the input is freq data with 204 channels, 30 frequencybins and 43 timebins
the call to "ft_combineplanar" took 0 seconds and required the additional allocation of an estimated 2 MB

Visualization

This part of the tutorial shows how visualize the results of any type of time-frequency analysis.

To visualize the event-related power changes, a normalization with respect to a baseline interval is desired. EEG and MEG data usually has much larger overall power in the lower frequencies than the higher frequencies. The consequence in a color-coded visual representation would be that the power at low frequencies would dominate the figure.

There are three possibilities for normalizing: # subtracting, for each frequency, the average power in a baseline interval from all power values. This gives, for each frequency, the absolute change in power with respect to the baseline interval. # expressing, for each frequency, the raw power values as the relative increase or decrease with respect to the power in the baseline interval. This means active period/baseline. Note that the relative baseline is expressed as a ratio; i.e. no change is represented by 1. # take the previous relative increase and subtract 1, i.e. express it as relative change. A value of 0.1 means 10% increase, a value of -0.1 means 10% decrease.

cfg = [];
cfg.layout        = 'neuromag306cmb.lay';
cfg.baselinetype  = 'relchange';
cfg.baseline      = [-inf 0];
cfg.renderer      = 'painters'; % the default is OpenGL, which is faster, but it gives visual artifacts on some computers
cfg.colorbar      = 'yes';
figure
ft_multiplotTFR(cfg, TFRhann_cmb);
reading layout from file neuromag306cmb.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the input is freq data with 102 channels, 30 frequencybins and 43 timebins
the call to "ft_freqbaseline" took 0 seconds and required the additional allocation of an estimated 0 MB
Warning: (one of the) axis is/are not evenly spaced, but plots are made as if
axis are linear 
the call to "ft_multiplotTFR" took 3 seconds and required the additional allocation of an estimated 0 MB

Note that by using the options cfg.baseline and cfg.baselinetype when calling plotting functions, baseline correction can be applied to the data. Baseline correction can also be applied directly by calling ft_freqbaseline. See the plotting tutorial for more details.

INSTRUCTION: Explore the TFR in detail, especially over the occipital cortex. Do you see gamma-band activity?

The figures use a relative baseline (cfg.basline='relative'), which means that power is expressed relative to the baseline. A value of 0.5 means that the power has halved, a value of 2 means that the power has doubled.

NOTE: The figures by default use an automatic scaling of the colorbar along the min-max range. The colors in one figure can therefore be different than in the other figure. You can specify cfg.zlim to control the colorbar. Please type 'colorbar'

There are three ways of graphically representing the data: # time-frequency plots of all channels, in a quasi-topographical layout # time-frequency plot of an individual channel, # topographical 2-D map of the power changes in a specified time-frequency interval. Starting from the first (ft_multiplotTFR) which shows all data, you can use the interactive feature of the figure to get the second (ft_singleplotTFR) and third (ft_topoplotTFR). For fine-grained control of the figures, e.g. for publications, you would not use the interactive features of the figure, but rather call ft_topoplotTFR itself.

(2) Time-frequency analysis with Hanning taper, frequency-dependent window length

It is also possible to calculate the TFRs with a sliding time-window that varies with frequency. For low frequencies the oscillations are long in time, which requires a long window. For high frequencies you can use a shorter window. The main advantage of this approach is that the temporal smoothing decreases with higher frequencies; however, this is at the expense of frequency smoothing. We will here show how to perform this analysis with a Hanning window. The approach is very similar to wavelet analysis, but differs in the choise of the tapering function. A Morlet wavelet analysis uses a Gaussian-shaped taper, which in principle is infinitely wide, whereas a Hanning taper has a clearly limited width.

The analysis is best done by first selecting the numbers of cycles per time window, which will be the same for all frequencies. For instance if the number of cycles per window is 7, the time window is 1000 ms for 7 Hz (1/7 x 7 cycles); 700 ms for 10 Hz (1/10 x 7 cycles) and 350 ms for 20 Hz (1/20 x 7 cycles). The frequency can be chosen arbitrarily - however; too fine a frequency resolution is just going to increase the redundancy rather than providing new information.

Below is the cfg for a 7-cycle time window.

cfg              = [];
cfg.method       = 'mtmconvol';
cfg.taper        = 'hanning';
cfg.output       = 'pow';
cfg.pad          = 2.5;
cfg.foi          = 2:4:120;
cfg.t_ftimwin    = 7./cfg.foi;  % 7 cycles per time window
cfg.toi          = -0.750:0.050:1.350;

Let us look at the relation between length of the time window versus frequency before computing the TFR:

figure
plot(cfg.foi, cfg.t_ftimwin)
xlabel('frequency (Hz)');
ylabel('window length (s)')
TFRhann7 = ft_freqanalysis(cfg, resampled);
the input is raw data with 204 channels and 125 trials
Warning: the trial definition in the configuration is inconsistent with the
actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive
segments of a continuous recording 
Warning: correcting numerical inaccuracy in the time axes 
Warning: correcting numerical inaccuracy in the time axes 
processing trials
Warning: output time-bins are different from input time-bins 
trial 125, frequency 30 (118.00 Hz), 1 tapers

the call to "ft_freqanalysis" took 52 seconds and required the additional allocation of an estimated 0 MB

Since this step takes quite some time, you might want to save the data again.

save visual_TFRhann7.mat TFRhann7

cfg = [];
TFRhann7_cmb = ft_combineplanar(cfg, TFRhann7);
the input is freq data with 204 channels, 30 frequencybins and 43 timebins
the call to "ft_combineplanar" took 0 seconds and required the additional allocation of an estimated 0 MB
cfg = [];
cfg.layout        = 'neuromag306cmb.lay';
cfg.baselinetype  = 'relchange';
cfg.baseline      = [-inf 0];
cfg.renderer      = 'painters'; % the default is OpenGL, which is faster, but it gives visual artifacts on some computers
cfg.colorbar      = 'yes';
figure
ft_multiplotTFR(cfg, TFRhann7_cmb);
reading layout from file neuromag306cmb.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the input is freq data with 102 channels, 30 frequencybins and 43 timebins
the call to "ft_freqbaseline" took 0 seconds and required the additional allocation of an estimated 0 MB
Warning: (one of the) axis is/are not evenly spaced, but plots are made as if
axis are linear 
the call to "ft_multiplotTFR" took 2 seconds and required the additional allocation of an estimated 1 MB

Note the boundary effects for lower frequencies (the even blue time frequency points in the plot): there is no power value calculated for these time frequency points. The power value is assigned to the middle time point in the time window and the time window extends in both directions. For example for 2 Hz the time window has a length of 3.5 sec (1/2 * 7 cycles = 3.5 sec), which means 1.75 seconds in both directions from the centre of the 'tile'. There is simply not enough data to fill that time window and the power will be set to NaN (not-a-number). The higher the frequency, the shorter the time window, the closer you can get to the edges of the data window. It is important to consider these boundary effects in your experimental design and during preprocessing to get all desired time-frequency points for your time window of interest.

(3) Time-frequency analysis using Morlet wavelets

An alternative to calculating TFRs with the mtmconvol and frequency-dependent time windows is to use Morlet wavelets. The approach is nearly similar to calculating TFRs with time windows that depend on frequency using a taper with a Gaussian shape, except that the Gausian function is infinitely wide. The commands below illustrate how to do this.

One crucial parameter to set is cfg.gwidth. It determines the width of the wavelets in number of cycles. Making the value smaller will increase the temporal resolution at the expense of frequency resolution and vice versa. The spectral bandwidth at a given frequency F is equal to F/width*2 (so, at 30 Hz and a width of 7, the spectral bandwidth is 30/7*2 = 8.6 Hz) while the wavelet duration is equal to width/F/pi (in this case, 7/30/pi = 0.074s = 74ms).

The parameter cfg.width specifies where the Gausian taper will be truncated (cut of and set to zero). So with cfg.gwidth=3 and cfg.width=5, there are three cycles under the bell-shaped curve from -1 to +1 times the standard-deviation, and the bell-shaped curve is truncated at 5 cycles. Ideally you would set cfg.width to infinite to get a proper Morlet wavelet decomposition, but since the data is not infinitely long, this cannot be achieved. Furthermore, cfg.width puts a hard constraint on the temporal characteristics and you can use it to be sure that a stimulus artefact (or response) does not inadvertedly leak into power estimates at other latencies.

Calculate TFRs using Morlet wavelets:

cfg = [];
cfg.method  = 'wavelet';
cfg.pad     = 4;
cfg.toi     = -0.750:0.050:1.500;
cfg.foi     = 5:2.5:120;
cfg.gwidth  = 3;
cfg.width   = 5;
TFRwavelet  = ft_freqanalysis(cfg, resampled);
the input is raw data with 204 channels and 125 trials
Warning: the trial definition in the configuration is inconsistent with the
actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive
segments of a continuous recording 
Warning: correcting numerical inaccuracy in the time axes 
Warning: correcting numerical inaccuracy in the time axes 
processing trials
Warning: output time-bins are different from input time-bins 
trial 125, frequency 47 (120.00 Hz)

the call to "ft_freqanalysis" took 164 seconds and required the additional allocation of an estimated 14 MB

Since this step takes quite some time, you might want to save the data afterwards

save visual_TFRwavelet.mat TFRwavelet

Again we will visualize the data after combining the planar channels:

cfg = [];
TFRwavelet_cmb = ft_combineplanar(cfg, TFRwavelet);

cfg = [];
cfg.layout        = 'neuromag306cmb.lay';
cfg.baselinetype  = 'relchange';
cfg.baseline      = [-inf 0];
cfg.renderer      = 'painters'; % the default is OpenGL, which is faster, but it gives visual artifacts on some computers
cfg.colorbar      = 'yes';
figure
ft_multiplotTFR(cfg, TFRwavelet_cmb);
the input is freq data with 204 channels, 47 frequencybins and 46 timebins
the call to "ft_combineplanar" took 0 seconds and required the additional allocation of an estimated 0 MB
reading layout from file neuromag306cmb.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the input is freq data with 102 channels, 47 frequencybins and 46 timebins
the call to "ft_freqbaseline" took 0 seconds and required the additional allocation of an estimated 0 MB
Warning: (one of the) axis is/are not evenly spaced, but plots are made as if
axis are linear 
the call to "ft_multiplotTFR" took 2 seconds and required the additional allocation of an estimated 0 MB

(4) Time-frequency analysis using multi-tapering

Finally we will demonstrate how to use the mtmconvol method with multitapering. Multitapering allows us to fully control the division of the the time-frequency space into 'tiles'. Here we will use a constant-length time window, with an increasing smoothing for the higher frequencies.

cfg = [];
cfg.output     = 'pow';
cfg.method     = 'mtmconvol';
cfg.taper      = 'dpss';
cfg.pad        = 2.5;
cfg.foi        = 10:5:120;
cfg.toi        = -0.750:0.05:1.350;
cfg.t_ftimwin  = 0.2 * ones(size(cfg.foi));  % use a fixed window length of 0.2 seconds
cfg.tapsmofrq  = nan(size(cfg.foi));  % start with NaN for all frequencies
cfg.tapsmofrq(:)          = 20;       % use 20Hz smoothing for the high frequencies
cfg.tapsmofrq(cfg.foi<50) = 10;       % use 10Hz smoothing for the intermediate frequencies
cfg.tapsmofrq(cfg.foi<30) = 5;        % use 5Hz smoothing (minimum) for the low frequencies

TFRmult = ft_freqanalysis(cfg, resampled);
the input is raw data with 204 channels and 125 trials
Warning: the trial definition in the configuration is inconsistent with the
actual data 
Warning: reconstructing sampleinfo by assuming that the trials are consecutive
segments of a continuous recording 
Warning: correcting numerical inaccuracy in the time axes 
Warning: correcting numerical inaccuracy in the time axes 
processing trials
Warning: output time-bins are different from input time-bins 
10 Hz: WARNING: using only one taper for specified smoothing
15.2 Hz: WARNING: using only one taper for specified smoothing
20 Hz: WARNING: using only one taper for specified smoothing
25.2 Hz: WARNING: using only one taper for specified smoothing
trial 1, frequency 22 (115.20 Hz), 7 tapers
10 Hz: WARNING: using only one taper for specified smoothing
...

Since this step takes quite some time, you might want to save the data afterwards.

save visual_TFRmult.mat TFRmult
cfg = [];
TFRmult_cmb = ft_combineplanar(cfg, TFRmult);

cfg = [];
cfg.layout        = 'neuromag306cmb.lay';
cfg.baselinetype  = 'relchange';
cfg.baseline      = [-inf 0];
cfg.renderer      = 'painters'; % the default is OpenGL, which is faster, but it gives visual artifacts on some computers
cfg.colorbar      = 'yes';
figure
ft_multiplotTFR(cfg, TFRmult_cmb);
the input is freq data with 204 channels, 23 frequencybins and 43 timebins
the call to "ft_combineplanar" took 0 seconds and required the additional allocation of an estimated 1 MB
reading layout from file neuromag306cmb.lay
the call to "ft_prepare_layout" took 0 seconds and required the additional allocation of an estimated 0 MB
the input is freq data with 102 channels, 23 frequencybins and 43 timebins
the call to "ft_freqbaseline" took 0 seconds and required the additional allocation of an estimated 0 MB
Warning: (one of the) axis is/are not evenly spaced, but plots are made as if
axis are linear 
the call to "ft_multiplotTFR" took 2 seconds and required the additional allocation of an estimated 1 MB

The time available for the first analysis of this tutorial dataset has been limited. None of the time-frequency analyses demonstrated above will therefore be fully optimal for this experimental design and for the features present in the data of this subject. However, we do hope that this tutorial demonstrates the choices that can be made in optimizing the analysis.

Optimizing the analysis involves tweaking the parameters to best capture the features of the data. If the data has different features at lower than at higher frequencies, you should not hesitate splitting the analysis over the two ranges. For this particular dataset, the most optimal analysis probably will consist of a wavelet or Hanning-tapered analysis for the lower frequencies (say, up to 30 Hz) and a multi-taper analysis for the higher frequencies.

Summary and suggested further reading

This tutorial showed how to do time-frequency analysis on a single subject MEG dataset and how to plot the time-frequency representations. There were 4 methods shown for calculating time- frequency representations and three functions for plotting the results.

After having finished this tutorial on time-frequency analysis, you can continue with the beamformer source reconstruction tutorial if you are interested in the source-localization of the power changes or the cluster-based permutation tests on time-frequency data tutorial if you are interested how to do statistics on the time-frequency representations.