Frequency-domain analysisFollowing Fourier ideas, EEG signals can be represented in the time domain (when recording the brain activity from scalp EEG surface electrodes) or alternatively in the frequency domain (when decomposing time-domain signals into weighted sums of sine and cosine functions). In the frequency-domain analysis (but also in the time-frequency domain analysis), EEG signals are represented in terms of oscillatory components, each of which is believed to isolate the activity of localized neuronal populations. To describe any oscillatory component, we estimate its amplitude and phase. The amplitude refers to the magnitude of an oscillation and can be measured in units of signal amplitude (µV), in power (µV2) or in decibel units (20*log10(µV) or 10*log10(µV2)). The phase refers to the position of a point in time on a waveform cycle and is commonly measured in degrees or radians. The Fourier Transform is usually employed to represent a signal in the frequency-domain. Especially after the introduction of the Fast Fourier Transform (a very fast and efficient algorithm to calculate the Fourier Transform), the Fourier Transform has become one of the most used tools in different scientific disciplines. Let us consider a digital signal x[n], n = 0,…,N-1, which has been obtained from a continuous signal x(t) by sampling at equal time intervals Δt (i.e., with a sampling frequency fs = 1/Δt). The discrete Fourier Transform is given by the following formula:
|
Return to the wiki
|
|
|
|
|
The signal x[n] can be reconstructed with the inverse discrete Fourier Transform:
|
|
|
|
|
|
The Fourier coefficients X[k] are complex numbers that can be represented both in Cartesian and polar forms.
|
|
|
|
|
|
where XR and XI denote the real and imaginary parts in the Cartesian representation, whereas |X[k]| and φ denote the amplitude and phase in the polar representation. If we restrict our attention only to signal power, from the complex Fourier coefficients we can compute the periodogram as: |
|
|
|
|
|
|
|
|
The periodogram is a raw estimation of the power spectrum of the signal, provided that the signal is stationary (i.e., the main characteristics of the signal, such as its mean and variance do not change over time). To illustrate why the periodogram is a raw estimate of the power spectrum, let us consider the following example.
|
|
Figure 4.1: Spectral leakage
|
|
|
The sinusoid on the left panel of Figure 4.1 has an exact number of cycles in the 2 s period of the signal. Its corresponding periodogram (on the left bottom side) shows a single peak at 3 Hz. On the other hand, the sinusoid on the right has a non-integer number of cycles and its periodogram shows an activity that is spread between 2 and 5 Hz. This phenomenon, consisting of smearing the power spectrum estimation, is called leakage (for further details, see Oppenheim and Schaffer, 1999). A simple way to avoid leakage would be to take an integer number of cycles. However, in general it is not possible to define a single periodicity because real signals contain activity at multiple frequencies. A useful approach for reducing the leakage effects is to use an appropriate window function that tapers the borders of the signal. This approach is referred to as a ‘windowing’. Several window functions have been proposed and their advantages and disadvantages depend on the specific application. Among the most popular window functions are the Hanning, Hamming, Barlett and Blackman. Another reason why the periodogram is a raw estimate of the power spectrum is the following. Because the variance of the periodogram does not go to zero, even in the limit of an infinite sample size, the periodogram is not a statistically consistent estimate of the power spectrum. For large sample sizes, moreover, the periodogram tends to vary rapidly with frequency, and the resulting power spectrum tends to look like a random pattern. To obtain a smoother estimate of the power spectrum, Barlett (1953) proposed to divide the signal into segments, to calculate the periodogram for each segment, and then to average the periodograms. Later on, Welch (1967) showed that better estimates can be achieved by using half-overlapping segments. The number of segments to be used for averaging the periodograms depends on the specific application. Ideally, it should be more than 30. As regards the length of the segments, it is noteworthy to point out that the larger the segments, the better the frequency resolution (the frequency resolution Δf is defined by the formula Δf = 1 / (N*Δt), where N is the number of data points and Δt is the time resolution; note that N*Δt is the duration of the signal). However, it is advisable to use segments that are not too long, so that the signal can be considered stationary to a first approximation. Now we will upload the dataset that we employed to illustrate the preprocessing steps, and we will estimate the power spectrum using the Welch method (click here to download all datasets). |
|
%% Initializations
dataPath = 'C:\Users\CIMeC\Desktop\wikiDatasets\'; % It depends on where
% the wikiDatasets
% folder is stored
fileName = 'S01RestPreprocessed.set';
dataSet = fullfile(dataPath,fileName);
|
↑ Go up |
%% Defining trials (or epochs)
% Reading the data as one long continuous segment
cfg = [];
cfg.dataset = dataSet; % file path
data = ft_preprocessing(cfg);
% Segmenting the data into trials
cfg = [];
cfg.length = 2; % seconds
cfg.overlap = 0.5; % value between 0 and 1 (percent)
dataSeg = ft_redefinetrial(cfg,data);
% Reading the events from the data
event = ft_read_event(dataSet);
% Removing trials containing discontinuities
trialsToReject = zeros(length(dataSeg.sampleinfo),1);
for eventId = 1:length(event)
if isempty(event(eventId).type)
tmpSample = event(eventId).sample;
for sampleInfoId = 1:length(dataSeg.sampleinfo)
if (tmpSample - dataSeg.sampleinfo(sampleInfoId,1)) * ...
(tmpSample - dataSeg.sampleinfo(sampleInfoId,2)) < 0
trialsToReject(sampleInfoId) = sampleInfoId;
end
end
end
end
trialsToReject = trialsToReject(trialsToReject ~= 0);
trialsToKeep = 1:length(dataSeg.sampleinfo); % initialize vector with all
% trial indices
trialsToKeep(trialsToReject) = []; % create trial indices to keep
cfg = [];
cfg.trials = trialsToKeep;
dataSeg = ft_selectdata(cfg,dataSeg);
|
↑ Go up |
%% Frequency-domain analysis
% Computing power spectrum for each trial
cfg = [];
cfg.foilim = [0 45];
cfg.method = 'mtmfft';
cfg.taper = 'hanning';
cfg.keeptrials = 'yes';
cfg.channel = {'all' '-AFp9' '-AFp10'}; % all but ocular channels
freq = ft_freqanalysis(cfg,dataSeg);
% Computing the average over trials
cfg = [];
%cfg.keeptrials = 'yes';
freqAvg = ft_freqdescriptives(cfg,freq);
|
↑ Go up |
%% Plot PSD
freqAvg.powSpctrmDB = 10*log10(freqAvg.powspctrm); % convert to
% decibel (dB)
figure; plot(freqAvg.freq,freqAvg.powSpctrmDB')
set(gca,'ylim',[-30 15])
xlabel('Frequency (Hz)')
ylabel('Power spectral density (dB/Hz)')
% 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
% Topographic plot
cfg = [];
cfg.xlim = [8 12]; % alpha band
cfg.zlim = [-13 -1];
cfg.parameter = 'powSpctrmDB';
cfg.marker = 'off';
cfg.highlight = 'off';
cfg.comment = 'auto';
cfg.layout = layout;
cfg.colorbar = 'yes'; % or 'southoutside'
cfg.colormap = jet;
figure; ft_topoplotER(cfg,freqAvg);
c = colorbar;
c.LineWidth = 1;
c.FontSize = 14;
title(c,'dB')
|
↑ Go up |
|
Analogously to the ft_timelockgrandaverage function, the ft_freqgrandaverage allows computing the grand average spectra, including all individual subjects’ spectrum as an input argument. Differences in spectral power between conditions can be statistically tested using the function ft_freqstatistic. The configuration structure remains largely the same, as compared to the statistical evaluation of time-domain data described above. Again, we can use the function ft_topoplotER to show the outcome of the statistical comparisons (see the previous section). |
|
References
|
↑ Go up |