Time-frequency analysis

One of the main limitations of the frequency analysis is that the Fourier transform is not well suited for the processing of non-linear signals. It is also not optimal for characterizing transient responses, as in the case of evoked potentials. To estimate the power spectrum, we assume that the signal is stationary (EEG signals can be regarded as quasi-stationary only within temporal windows of a few seconds). Consequently, we assume that the oscillatory components at different frequencies are constant throughout the whole signal. When we analyze evoked potentials, we would like to track the oscillatory activity varying over time. Unfortunately, the Fourier transform cannot resolve time-varying features of the signal. In other words, the Fourier transform does not give any time information.

A simple, intuitive solution to the lack of time resolution is to chop the data into pieces, then to taper each piece with an appropriate window function, and finally to estimate the power spectrum for each windowed segment. This procedure is called the short-time Fourier transform (STFT). The stationarity assumption is fulfilled in good approximation, because the signals are quasi-stationary within each window. The bad news, however, is that choosing the window size is not so straightforward. If the window is too large, it will give a good frequency resolution (the frequency resolution is inversely proportional to the data length), but at the price of a reduced time localization. On the contrary, if the window is too narrow, we will have a good localization in time, but a poor localization in frequency. The trade-off between time and frequency resolution is called the uncertainty principle of signal analysis, in analogy to Heisenberg’s uncertainty principle, which states that it is not possible to determine the position and velocity of a given particle at the same time, with arbitrary accuracy.

The main limitation of the STFT is that the size of the tapering window is fixed, and the trade-off between time and frequency resolution may be optimal for a certain frequency, but not for others. As in the case of evoked potentials, a single window may not be sufficient, and ideally, we should set different window sizes for different frequency ranges.

In the late 1970s, Jean Morlet proposed a straightforward solution to the main limitation of the STFT. Morlet worked with a cosine function tapered with a Gaussian function (i.e., a Gabor function), and to get the different frequencies, he compressed or stretched the Gabor function in time. In this way, he manipulated the same wave shape at different scales, in fact introducing a variable size. This simple trick laid the foundation of modern wavelets.

The key idea of wavelet analysis is to consider different window sizes for different frequencies by compressing or stretching the same function, which is called the mother wavelet. In principle, the mother wavelet does not necessarily need to be a sinusoid, and, as usual, the choice of the appropriate wavelet function depends on the application in mind and on the type of patterns that we want to localize in time and frequency. Since in the wavelet analysis we adopt functions at different scales, rather than sinusoids of different frequencies, it would be more accurate to say that the result is a time-scale decomposition, rather than a time-frequency decomposition.

The main advantage of using wavelets is that the time-frequency resolution is variable with frequency. For low frequencies, the wavelet has a large time window, thus increasing the frequency resolution at the cost of time resolution. Conversely, for high frequencies, the time window is small, thus giving a better time resolution at the cost of larger frequency uncertainty.

In FieldTrip, the time-frequency decomposition is computed by the function ft_freqanalysis, specifying one of the supported algorithms in cfg.method. The method mtmconvol implements a STFT, whereas the method wavelet is based on Morlet wavelets (a Morlet wavelet is a complex sinusoid tapered with a Gaussian function).

Before illustrating how to perform Morlet wavelet analysis in FieldTrip, it is important to point out that M/EEG researchers distinguish between three types of oscillatory activities: spontaneous, induced and evoked oscillations. Spontaneous oscillations occur independently of the presentation of any external stimulus; induced oscillations are elicited by the onset of an external stimulus, but are not phase-locked to the stimulus onset; evoked oscillations are elicited by the onset of a stimulus and are phase-locked to the stimulus onset. Because the phase (i.e., timing) may vary across trials, induced oscillations can be distorted or even obscured, if they are analyzed on ERPs, namely after the average across trials is taken in the time-domain. For this reason, a common approach to the time-frequency analysis is to estimate the time-frequency representation of each single trial, and then compute an average over trials.

In general, the average of the spectral decomposition over trials contains both evoked and induced oscillations. Like ERPs, spectral decompositions can be measured with respect to some pre-stimulus baseline but, unlike ERPs, spectral decompositions can be also expressed as a proportion with respect to a baseline period. We usually interpret an increase in signal power (amplitude^2), relative to baseline, to reflect an increase in the synchronicity with which a neural population is firing, or to reflect an increase in the size of the neural population firing synchronously. An increase in signal power relative to baseline is also called ‘event-related synchronization’. Likewise, event-related de-synchronizations are often interpreted to reflect decreases in the synchronization or size of a neural population with respect to a baseline.

Return to the wiki
home page

%% Initializations

dataPath = 'C:\Users\CIMeC\Desktop\wikiDatasets\'; % It depends on where the
                                                   % wikiDatasets folder is
                                                   % stored
fileNameDescription = 'Oddball.set';
 
subjects = {'S01' 'S02' 'S03' ...
            'S04' 'S05'}; % the variable subjects contains the list of 
                          % subject-specific directories
nSubjs = numel(subjects);

↑ Go up

%% Time-frequency domain analysis

% Standard tones
for i = 1:nSubjs
    fileName = 'dataStd.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'data'
    
    % Computing time-frequency representation for each trial
    cfg = [];
    cfg.method = 'wavelet';
    cfg.foi = 1:1:40; % 1 = frequency resolution = 1/trial duration
    cfg.toi = -0.05:0.05:0.5;
    cfg.width = 3; % width of the wavelet ( = number of cycles). 
                   % Small values increase the temporal resolution at 
                   % the expense of frequency resolution
    cfg.keeptrials = 'yes';
    cfg.channel = {'all'};
    cfg.outputfile = fullfile(dataPath,subjects{i},'timeFreqStd');
    ft_freqanalysis(cfg,data);
    clear data
end
 
timeFreqAvgStd = cell(1,nSubjs);
for i = 1:nSubjs
    fileName = 'timeFreqStd.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'freq'
    
    % Computing the average over trials
    cfg = [];
    timeFreqAvgStd{i} = ft_freqdescriptives(cfg,freq);
    clear freq
end
 
bslnTimeFreqAvgStd = cell(1,nSubjs);
for i = 1:nSubjs
    fileName = 'timeFreqStd.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'freq'
    
    % Computing the change relative to a baseline
    cfg = [];
    cfg.baselinetype = 'relchange'; % 'absolute','relative','relchange' or 'db'
    cfg.baseline = [-0.05 0];
    bslnTimeFreqAvgStd{i} = ft_freqbaseline(cfg,freq);
    
    % Computing the average over trials
    cfg = [];
    bslnTimeFreqAvgStd{i} = ft_freqdescriptives(cfg,bslnTimeFreqAvgStd{i});
    clear freq
end
 
% Deviant tones
for i = 1:nSubjs
    fileName = 'dataDev.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'data'
    
    % Computing time-frequency representation for each trial
    cfg = [];
    cfg.method = 'wavelet';
    cfg.foi = 1:1:40; % 1 = frequency resolution = 1/trial duration
    cfg.toi = -0.05:0.05:0.5;
    cfg.width = 3; % width of the wavelet ( = number of cycles). 
                   % Small values increase the temporal resolution at 
                   % the expense of frequency resolution
    cfg.keeptrials = 'yes';
    cfg.channel = {'all'};
    cfg.outputfile = fullfile(dataPath,subjects{i},'timeFreqDev');
    ft_freqanalysis(cfg,data);
    clear data
end
 
timeFreqAvgDev = cell(1,nSubjs);
for i = 1:nSubjs
    fileName = 'timeFreqDev.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'freq'
    
    % Computing the average over trials
    cfg = [];
    timeFreqAvgDev{i} = ft_freqdescriptives(cfg,freq);
    clear freq
end
 
bslnTimeFreqAvgDev = cell(1,nSubjs);
for i = 1:nSubjs
    fileName = 'timeFreqDev.mat';
    dataSet = fullfile(dataPath,subjects{i},fileName);
    load(dataSet); % loading the variable 'freq'
    
    % Computing the change relative to a baseline
    cfg = [];
    cfg.baselinetype = 'relchange'; % 'absolute','relative','relchange' or 'db'
    cfg.baseline = [-0.05 0];
    bslnTimeFreqAvgDev{i} = ft_freqbaseline(cfg,freq);
    
    % Computing the average over trials
    cfg = [];
    bslnTimeFreqAvgDev{i} = ft_freqdescriptives(cfg,bslnTimeFreqAvgDev{i});
    clear freq
end

↑ Go up

%% Grand average time-frequency data
% Standard tones
cfg = [];
gAvgTimeFreqStd = ft_freqgrandaverage(cfg,timeFreqAvgStd{:});
gAvgBslnTimeFreqStd = ft_freqgrandaverage(cfg,bslnTimeFreqAvgStd{:});
 
% Deviant tones
cfg = [];
gAvgTimeFreqDev = ft_freqgrandaverage(cfg,timeFreqAvgDev{:});
gAvgBslnTimeFreqDev = ft_freqgrandaverage(cfg,bslnTimeFreqAvgDev{:});

↑ Go up

%% Plotting time-frequency analysis
% Defining channel layout
templateLayout = 'EEG1005.lay'; % one of the template layouts included
                                % in FieldTrip
cfg = [];
cfg.layout = which(templateLayout); 
layout = ft_prepare_layout(cfg);
 
% Increasing layout width and height
if strcmp(templateLayout,'EEG1005.lay')
    layout.width = 0.07 * ones(length(layout.width),1);
    layout.height = 0.04 * ones(length(layout.height),1);
end
 
% Multi plot
cfg = [];
cfg.layout = layout; %'easycapM1.mat';
cfg.interactive = 'yes';
cfg.showoutline = 'yes';
cfg.showlabels = 'yes';
cfg.fontsize = 10;
cfg.comment = sprintf('\n');
cfg.xlim = [-0.05 0.5];
cfg.ylim = [8 40]; % we cannot analyze delta and theta band because 
                   % the trial lenght is too short
cfg.zlim = [-6 6]; 
cfg.colorbar = 'yes'; % or 'southoutside'
cfg.colormap = jet; % 'parula' or 'jet'
figure
ft_multiplotTFR(cfg,gAvgBslnTimeFreqStd);
title('Grand average - Standard (relchange)');
figure
cfg.zlim = [-6 6];
ft_multiplotTFR(cfg,gAvgBslnTimeFreqDev);
title('Grand average - Deviant (relchange)');
 
% Plotting single average across electrodes of interest
cfg = [];
cfg.ylim = [8 40];
cfg.zlim = [-6 6];
cfg.channel = {'Fz','FCz','Cz'};
cfg.colormap = jet;
figure;
ft_singleplotTFR(cfg,gAvgBslnTimeFreqStd);
set(gca,'Fontsize',20);
title('Mean over fronto-central electrodes');
set(gca,'box','on');
xlabel('time (s)');
ylabel('frequency (Hz)');
c = colorbar;
c.LineWidth = 1;
c.FontSize  = 18;
title(c,'Rel. change');

↑ Go up

%% Statistical analysis of time-frequency representations at the group level
% Creating a neighbourhood structure to define how spatial clusters are formed
fileName = 'dataStd.mat';
dataSet = fullfile(dataPath,subjects{1},fileName);
load(dataSet); % loading the variable 'data'
 
templateElec = 'standard_1005.elc';
elec = ft_read_sens(which(templateElec));
 
idx = ismember(elec.label,data.label);
elec.chanpos = elec.chanpos(idx,:);
elec.chantype = elec.chantype(idx);
elec.chanunit = elec.chanunit(idx);
elec.elecpos = elec.elecpos(idx,:);
elec.label = elec.label(idx);
 
cfg = [];
cfg.method = 'distance'; %'template' or 'distance' or 'triangulation'
cfg.neighbourdist = 80; % maximum distance between neighbouring sensors
                        % (only for 'distance')
cfg.feedback = 'yes'; % show a neighbour plot
neighbours = ft_prepare_neighbours(cfg,elec);
 
% Preparing the design matrix for the statistical evaluation
% For within-subjects analysis, the design matrix contains two rows
design = [1:nSubjs 1:nSubjs; ones(1,nSubjs) ones(1,nSubjs)*2]; 
 
% Test the null-hypothesis of exchangeability between conditions
cfg = [];
cfg.channel = {'all'};
cfg.avgoverfreq = 'no';
cfg.method = 'montecarlo';
cfg.clusterthreshold = 'nonparametric_common';
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan = 2;
cfg.neighbours = neighbours;
cfg.tail = 0;
cfg.clustertail = 0;
cfg.alpha = 0.025;
cfg.numrandomization = 1000;
cfg.statistic = 'ft_statfun_depsamplesT';
cfg.design = design;
cfg.uvar = 1; % row in design indicating subjects, repeated measure
cfg.ivar = 2; % row in design indicating condition for contrast
statTimeFreq = ft_freqstatistics(cfg,bslnTimeFreqAvgStd{:},bslnTimeFreqAvgDev{:});
 
% Plotting stat output over fronto-central electrodes
cfg = [];
cfg.layout = layout;
cfg.channel = {'Fz','FCz','Cz'};
cfg.parameter = 'stat';
cfg.colormap = jet;
cfg.ylim = [8 40];
cfg.marker = 'off';
cfg.style = 'fill';
cfg.comment = 'off';
cfg.maskparameter = 'mask';
cfg.maskstyle = 'outline';
cfg.colorbar = 'yes';
cfg.zlim = [-5 5];
figure;
ft_singleplotTFR(cfg,statTimeFreq);
set(gca, 'Fontsize',20);
title ('Mean over fronto-central electrodes');
set(gca,'box','on');
xlabel('time (s)');
ylabel('frequency (Hz)');
c = colorbar;
c.LineWidth = 1;
c.FontSize  = 18;
title(c,'t-val');

↑ Go up