Simplified introduction to time-frequency analysis in EEGLAB

  -->

Version 1.1, updated 14th February 2013

This simplified account of time-frequency analysis was written by a non-expert who was learning to use the newtimef() command of EEGLAB. I suspect that there are others, like me, who come to EEGLAB with a background in analysis of averaged ERPs and who find the account of time-frequency analysis in the EEGLAB manual assumes more background knowledge than they have. As I gradually battled throught the confusion, I felt it might help others if I laid out what I understood in simple terms.

If you find errors in this account, please let me know. However, please be aware that, because of my lack of expertise, I am unlikely to be able to answer technical questions.

To show you how time-frequency analysis works in EEGLAB, we will first create some simulated data, and then run the analysis with different parameters. The best way to learn about the parameters of time-frequency analysis is to modify the Matlab scripts below and experiment with settings.

Creating a simple sine wave

To create a sine wave, you need four parameters:

D = duration of the signal, here specified in seconds

S = sampling rate, samples per second

F = frequency, in Hz

P = phase, as an angle

Here's a little program to illustrate the method

%-------------------------------------------------------------------------------------

% make_sig.m

%-------------------------------------------------------------------------------------

D = 1; % signal duration, 1 second

S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave

F = 50; % frequency in Hz

P = .25 % NB. 2 pi radians = 360 degrees

T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals

t = [T:T:D]; % time vector %NB corrected from previous version

myphi=2*pi*P;

mysig = sin(2*F*t*pi+ myphi); % sinusoidal signal; 2*F*t*pi is freq in radians

figure; plot(t,mysig);

xlabel('Time (seconds)');

ylabel('Amplitude');

%sound(mysig); %You can play the sound if you like!

%-------------------------------------------------------------------------------------

If you don't understand about sampling rate, frequency or phase, experiment with altering values of S, F and P to see the effect on the wave. Hint: it is easiest to see how changing phase alters the signal if you select a frequency less than 10.

Note that if you select a value of S that is less than twice F, the waveform changes, because the sampling is too slow to capture the fluctuations. 2*F is known as the Nyquist frequency, which is the lowest sampling frequency that can be used to correctly identify a sine wave of frequency F.

Next step is to add some sine waves together, with different frequencies and phases.

Creating a complex wave

The next program is based on make_sin, but adds together waves of different frequency and phase.

%-----------------------------------------------------------------------

% create_complex.m

%-----------------------------------------------------------------------

D = 1.2; %signal duration

S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave

F = [30 40 100 200]; % 4 frequencies in Hz

w = 2*pi*F; % convert frequencies to radians

P = [0 .5 .25 .3]; % 4 corresponding phases

A = [1 .5 .3 .2]; % corresponding amplitudes

T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals

t = [T:T:D]; % time vector %NB this has been corrected from previous version

mysig=zeros(1,length(t)); %initialise mysig

myphi=2*pi*P; % compute phase angle

nfreqs=length(F); % N frequencies in the complex

% Add all sine waves together to give composite

for thisfreq=1:nfreqs

mysig = mysig+A(thisfreq)*(sin(w(thisfreq)*t + myphi(thisfreq)));

end

figure; plot(t,mysig);

xlabel('Time (seconds)');

ylabel('Amplitude');

%-----------------------------------------------------------------------

Fourier transform

The Discrete Fourier Transform (FFT) is a method that takes a complex waveform and decomposes it into a set of component sine waves. There's a clear explanation of Matlab's FFT transform here:

The FFT requires a dataset that is 2^N elements long, i.e. 2, 4, 8, 16, 32, etc., so the transform uses the value of N that is closest to the signal length you have (NFFT). This is computed below for illustration, but not used.

%-----------------------------------------------------------------------

% do_FFT.m

%-----------------------------------------------------------------------

% assumes you have already created mysig.

L = length(mysig);

myFFT = fft(mysig,S);

myFFT=myFFT/L; %scale the output to 1

freq = S/2*linspace(0,1,S/2);%create a range with all frequency values

% Plot single-sided amplitude spectrum.

figure; stem(freq,abs(myFFT(1:length(freq))));

xlabel('Frequency (Hz)');

ylabel('Amplitude ');

%-----------------------------------------------------------------------

If you look at the myFFT variable, you will see it has a real and imaginary part. These correspond to the amplitude and phase of the signal. We select just the real part by using the abs() command.

Event-related spectral perturbation (ERSP)

The time-frequency analysis options of EEGLAB allow you to calculate event-related spectral perturbation (ERSP), which reflects the extent to which the power at different frequencies in a signal is altered in relation to a specific time point, such as signal onset.

It is easiest to gain an understanding of ERSP using constructed stimuli where the spectral content is known. Modify the create_complex.m script above to generate a signal, mysig, that is 2s long, with sampling rate 1000 Hz, composed of frequencies 10, 20 and 45 Hz. Have phase settings all at zero, and amplitude settings all at one.

Now modify the signal so that the power of the 20 Hz component is increased after the first 500 ms, by the following commands:

%-----------------------------------------------------------------------

% assumes you have already created mysig

A(2)=2; %amplitude of 2nd frequency increased from 1 to 2

mysig2=zeros(1,length(t)); %initialise mysig2

for thisfreq=1:nfreqs

mysig2 = mysig2+A(thisfreq)*(sin(w(thisfreq)*t + myphi(thisfreq)));

end

mysig(500:1000)=mysig2(500:1000);

t=t-.5; %subtract 500 ms to indicate first 500 ms are baseline, i.e. -ve time

figure; plot(t,mysig);

xlabel('Time (seconds)');

ylabel('Amplitude');

%-----------------------------------------------------------------------

In effect, we are simulating the situation where power is constant during a 500 ms baseline, but then increases for 500 ms when a stimulus is presented.

Now we are ready to see how newtimef displays ERSP for this stimulus, using the newtimef() command.

%------------------------------------------------------------------------------------

eeglab %load EEGLAB commands if you have not already done so

figure;[ersp,itc,powbase,times,freqs]=...

newtimef( mysig,2000,[-500 1500],1000, 0,'plotitc','off');

%------------------------------------------------------------------------------------

This gives a plot of the ERSP, which is shown as a colormap against time (x-axis) and frequency (y-axis). The default is to treat time values less than zero as the baseline. You will see a yellow patch in the region 0-500 ms at a point corresponding to 20 Hz. This indicates that the power in the signal went above the baseline value in this interval. The colorbar on the right hand side indicates the ERSP scale, which includes negative as well as positive values. A negative value denotes that there is a decrease in power relative to the baseline.

Notice that although we entered data for the time range from -500 to 1500 ms, the actual time range displayed is less than this. This is because the frequency analysis works by considering how a signal changes over time, and cannot give instantaneous values. You can see the time scale that is output by newtimef() by inspecting the 'times' term of newtimef() output. You can also see the centres of the frequency bands that the analysis generates by inspecting 'freqs'. You will see there are 24 frequencies ranging from 3.09 Hz to 50 Hz. Now look at the variable 'ersp'; you will see it is a 24 x 200 matrix of numbers, i.e. there is a value of ersp for each frequency bin and each time point. You can generate a plot resembling the ERSP plot using the imagesc command, which translates numbers into colors, with the command below. The [-30 30] sets the scale for the colors; you may want to see what happens if you omit or modify this.

figure;imagesc(ersp,[-30 30])

set(gca,'YDir','normal') % y axis is inverted if you don't specify this

h = colorbar;

Note that the ERSP plot is not quite what we might expect: there is both temporal and spectral splattering: the yellow region extends around from 20 Hz, and a yellow region at 30 Hz starts before zero. There is also a blue band suggesting a drop in power at 30-35 Hz It seems that these settings of newtimef() introduce some artefact. This is because the FFT is useful for computing a frequency spectrum over a long signal, but less good at identifying local frequency content.

The plot on the left-hand of the y-axis shows the overall signal power at different frequencies: as we would expect, this peaks at 10, 20 and 45 Hz.

The plot beneath the x-axis shows the averaged signal in the time domain.

Save mysig as mysig1, and ersp as ersp1, as you will want to use these again.

Re-run the program setting A(2) to .2, i.e. make the 20 Hz component substantially smaller during the 0-500 ms interval. Instead of yellow, we now see a blue region in this interval and frequency, denoting reduction of power.

The ERSP scale is in dB, because it is a ratio of the power in the epoch from 0-1500 ms to the average power in the baseline. Now reset mysig1 to mysig. To see the raw power, you can rerun the command as:

%------------------------------------------------------------------------------------

figure;[ersp,itc,powbase,times,freqs]=...

newtimef( mysig,2000,[-500 1500],1000, 0,'baseline',NaN,'plotitc','off');

%------------------------------------------------------------------------------------

The colors now depict raw log power, rather than ERSP (although the label of the colorbar still states ERSP). The ERSP in the initial analysis is found by finding the average log power in the baseline for each frequency, and subtracting these values from the low power for each timepoint in the post-zero period. (Footnote: I think this is not strictly true – my attempts to recreate the ERSP by subtraction gives slightly different results, and I'd welcome clarification on why this is so. However, note that the subtraction gives results that are sufficiently close to computed ERSP that they look very similar when plotted).

Let's now take a look at the arguments used with the newtimef() command.

The first argument is the dataset to be analysed. In the context of EEGLAB this will be an EEG dataset. Using the command pop_newtimef(), you will have the opportunity to specify the dataset, whether the analysis is done on channels or components, and which channel is analysed. An example is given below. For the time being, though, we are keeping things simple, and just using our toy simulation to explore the effect of altering other parameters of the newtimef() command.

The second argument is 'frames', i.e. number of points per epoch. We treated our simulated data as representing a single ERP epoch, but if we altered the 'frames' parameter, we could represent it as two consecutive epochs, i.e.

%-----------------------------------------------------------------------------------------

figure;[ersp,itc,powbase,times,freqs]=...

newtimef( mysig,1000,[-500 500],1000, 0,'plotitc','off');

%-----------------------------------------------------------------------------------------

Note that we then have to also change the next argument, which is the interval over which the power is computed, as our epoch length is now only 1000 points.

Try running this command and see if you can understand why the output looks different.

The fourth argument in the command is the sampling rate. Try and anticipate what will happen if you changed this to 500. Then try it out.

%-----------------------------------------------------------------------------------------

figure;[ersp,itc,powbase,times,freqs]=...

newtimef( mysig,2000,[-500 1500],500, 0,'plotitc','off');

%-----------------------------------------------------------------------------------------

You will see a change in the frequency location of the yellow area. This is because at the lower sampling rate, the distance between the points in the epoch is doubled.

The fifth argument, currently set to zero, specifies how the frequency analysis is done. When zero is used, the analysis uses the FFT transform, with Hanning window tapering. The effect of a Hanning window can be seen in red in the following plot of a Fourier output of a signal composed of five sine waves:

To get an analysis of frequency at different time-points, the FFT was applied to the signal over different time windows. Note that with this approach there is an inevitable trade-off between time resolution and frequency resolution. At one extreme we could take the whole epoch as a window, which would give a good representation of the frequencies in the signal, but no time resolution. Or we could apply the FFT to very short windows of a few ms, but this would not be useful for detecting the range of frequencies in the signal, because it would not give a big enough sample to detect oscillations.

Let us now try a different method of frequency analysis, using Morlet wavelets. This method is adopted if we use a positive number for the 5^th argument of newtimef(). Try the following script, which will analyse the same data with different values specified for the number of cycles in each Morlet wavelet. The first loop, with cycles set at zero, just repeats the FFT analysis.

%-----------------------------------------------------------------------------------------

for i=0:4

figure;[ersp,itc,powbase,times,freqs]=...

newtimef( mysig,2000,[-500 1500],1000,i,'plotitc','off','erspmax',30);

end

% erspmax specifies +/- range for color bar

% Important to set this when comparing settings, to over-ride automatic

% default

%-----------------------------------------------------------------------------------------

This plot shows the output when cycles = 1.

Note that this analysis removes some of the spurious bands of low and high ERSP that were seen in the FFT analysis. As the number of cycles increases, note that the analysis of lower frequencies drops out. If you try the analysis with cycles set to 10, you will see that the analysis only covers the frequency range from 40 to 50 Hz, and so misses the 20 Hz increase completely.

One advantage of wavelet analysis over a moving Fourier spectrum is that it the window size can be varied according to the frequency being analysed. This is best illustrated if we create a new signal in which there are bursts of increased amplitude at different time points for different frequencies. Here is a little routine that creates a signal with a 500 ms baseline, then bursts of increased power, each lasting 200 ms, at 0 ms, 400 ms and 800 ms post baseline, for frequencies 8 Hz, 25 Hz and 50 Hz.

%-----------------------------------------------------------------------

% Make complex signal with bursts of power at 8, 25 and 40 Hz

%-----------------------------------------------------------------------

D = 2; %signal duration 2s

S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave

F = [8 25 40]; % 4 frequencies in Hz

w = 2*pi*F; % convert frequencies to radians

P = [0 0 0]; % 4 corresponding phases

T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals

t = [T:T:D]; % time vector %corrected from previous version

nfreqs=length(F); % N frequencies in the complex

A=ones(nfreqs,length(t)); %amplitudes initialised at one

A(1,500:700)=3; %boost of increased amplitude at each time pt, nb 500 corresponds to end of baseline period

A(2,900:1100)=3;

A(3,1300:1500)=3;

% Add all sine waves together to give composite

for thispt=1:2000

for thisfreq=1:nfreqs

mysig(thisfreq,thispt) =A(thisfreq,thispt)*sin(w(thisfreq)*t(thispt));

end

end

mysig2=sum(mysig,1);

t=t-.5; %subtract 500 ms to indicate baseline is negative

figure; plot(t,mysig2);

xlabel('Time (seconds)');

ylabel('Amplitude');

for i=0:3

% First plot ERSP with constant values of cycles

scrsz = get(0,'ScreenSize'); %used to adjust figure size; not normally needed

figure('Position',[1 scrsz(4)/4 scrsz(3)/2 scrsz(4)/4])

[ersp,itc,powbase,times,freqs]=...

newtimef( mysig2,2000,[-500 1500],1000,i,'plotitc','off', 'erspmax',30,

);

end

%-----------------------------------------------------------------------

Now try this, which gives the output immediately below:

%-----------------------------------------------------------------------

i = [1 5];

[ersp,itc,powbase,times,freqs]=...

newtimef( mysig2,2000,[-500 1500],1000,i,'plotitc','off', 'erspmax',30);

%-----------------------------------------------------------------------

This time, the number of cycles used varies with the frequency being analysed. The differences between settings with this sample waveform are not huge, but nevertheless, it should be possible to see that there are advantages in specifying a range of cycles if you want to inspect ERSP over a wide frequency range.

Pad ratio

As already mentioned, there is an unavoidable tradeoff between frequency resolution and time resolution in time-frequency analysis. In the analyses done so far, the temporal specificity of the detecting ERSP is good, but the frequency precision is less good. Pad ratio is another parameter that can be modified to select an analysis most suitable for identifying key features in your data. Padding is literally that – it involves adding zeros to your signal to make it longer, which improves frequency resolution.

To demonstrate this point, we will first create a signal similar to mysig2, but with both long and short bursts of power at two frequencies. We will then run newtimef() with different values of padratio. Having read about padratio, I had anticipated that a higher padratio would give an analysis that homed in more precisely on the precise frequency level where the power increased. This example suugests not – the number of frequencies that are distinguished in the output increases (as you can see if you look at length(freqs)), but little else changes.

%-----------------------------------------------------------------------------------

% Make complex signal with long and short bursts of power at 8 and 45 Hz

%-----------------------------------------------------------------------------------

D = 2; %signal duration 2s

S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave

F = [8 8 45 45]; % 4 frequencies in Hz

w = 2*pi*F; % convert frequencies to radians

P = [0 0 0 0]; % 4 corresponding phases

T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals

t = [T:T:D]; % time vector %corrected from previous version

t=t(1:length(t)-1); % takes off extra point arising from starting at zero

nfreqs=length(F); % N frequencies in the complex

A=ones(nfreqs,length(t)); %amplitudes initialised at one

A(1,550:570)=3; %boost of increased amplitude at each time pt

A(2,700:900)=3;

A(3,950:970)=3;

A(4,1100:1300)=3;

% Add all sine waves together to give composite

for thispt=1:2000

for thisfreq=1:nfreqs

mysig(thisfreq,thispt) =A(thisfreq,thispt)*sin(w(thisfreq)*t(thispt));

end

end

mysig2=sum(mysig,1);

t=t-.5; %subtract 500 ms to indicate baseline is negative

figure; plot(t,mysig2);

xlabel('Time (seconds)');

ylabel('Amplitude');

scrsz = get(0,'ScreenSize');

for pr=[.5, 1,2,4,8,16] %pad ratio

cyc=[1 ]; % cycles (can experiment with different values)

figure('Position',[1 scrsz(4)/4 scrsz(3)/2 scrsz(4)/3])

[ersp,itc,powbase,times,freqs]=...

newtimef( mysig2,2000,[-500 1500],1000, cyc, 'erspmax', 20, 'plotitc', 'off', 'padratio',pr);

text(52,20,strcat('Pad ratio= ',num2str(pr),': N freqs=',num2str(length(freqs))));

end

%-----------------------------------------------------------------------------------------

Note that in each case, the ability of the analysis to detect the short burst of frequency at 8 Hz is poor, whereas the short 45 Hz burst is clearer.

ERSP with multiple epochs

We have illustrated methods for looking at ERSP with a single simulated trial, but in practice this method will usually be applied to a set of ERP epochs.

The simplest way to achieve this with real data saved as an epoched .set file is with the pop_newtimef() command, which will automatically determine factors such as sampling rate and epoch length from an EEG .set file.

pop_newtimef(INEEG, typeproc, num, tlimits,cycles)

Inputs:

INEEG - input EEG dataset

typeproc - type of processing.

1 process the raw data

0 process ICA components

num component or channel number

tlimits [mintime maxtime] (ms) sub-epoch time limits

cycles as for newtimef()

There are numerous optional arguments that can be applied, as for newtimef().

ERSP with subtracted datasets

It is common for ERP paradigms to use subtraction methodology, whereby the averaged waveform from one condition is subtracted from another condition. For instance, the mismatch negativity (MMN) is measured in auditory experiments where the waveform for a common standard repeating stimulus is subtracted from the waveform for a rarer deviant stimulus. Traditionally, analyses are done on subtracted averages. In the context of ERSP, this raises the question of whether one should compare the ERSP associated with standards and that with deviants, or first subtract the waveforms and then consider the ERSP. Does it make a difference, and if so, which is preferable?

To answer that question, I generated two waveforms, simulating standard and deviant waves, each composed of summed sinosoids at 10, 25 and 34 Hz. These frequency values were arbitrarily selected. The power at all frequencies was equal for the standard, but doubled at specific time windows for the deviant. Random noise was then added to each waveform, and the mismatch wave computed. The script is given below, followed by ERSP plots for (a) standard (b) deviant waveform (c) mismatch, created by subtracting the standard from the deviant, and (d) the ERSP for standard subtracted from ERSP from deviant.

%-----------------------------------------------------------------------

% Simulated mismatch

%-----------------------------------------------------------------------

D = 2; %signal duration 2s

S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave

F = [10 25 34]; % 4 frequencies in Hz

w = 2*pi*F; % convert frequencies to radians

P = [0 0 0 ]; % 4 corresponding phases

T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals

t = [T:T:D]; % time vector %corrected from previous version on 12 Jul 2010

t=t(1:length(t)-1); % takes off extra point arising from starting at zero

nfreqs=length(F); % N frequencies in the complex

A=ones(nfreqs,length(t)); %amplitudes initialised at one

A(1,500:570)=2; %boost of increased amplitude at each time pt

A(2,700:900)=2;

A(3,950:970)=2;

% Add all sine waves together to give composite

for thispt=1:2000

for thisfreq=1:nfreqs

mysig(thisfreq,thispt) =A(thisfreq,thispt)*sin(w(thisfreq)*t(thispt));

end

end

mystan=sum(mysig,1);

A(1,850:1050)=3; %boost of increased amplitude for freq 1

A(2,1100:1200)=3;

A(3,1350:1400)=3;

for thispt=1:2000

for thisfreq=1:nfreqs

mysig(thisfreq,thispt) =A(thisfreq,thispt)*sin(w(thisfreq)*t(thispt));

end

end

mydev=sum(mysig,1);

b=normrnd(0,.5,2000);

mystan=mystan+b(1,:);

mydev=mydev+b(2,:);

myMM=mydev-mystan;

t=t-.5; %subtract 500 ms to indicate baseline is negative

figure; plot(t,mystan);hold

plot(t,mydev,'color','r');

legend({'standard','deviant'})

xlabel('Time (seconds)');

ylabel('Amplitude');

scrsz = get(0,'ScreenSize');

cyc=[1 ]; % cycles

figure('Position',[1 scrsz(4)/4 scrsz(3)/2 scrsz(4)/3])

[ersp,itc,powbase,times,freqs]=...

newtimef( mystan,2000,[-500 1500],1000,cyc,'erspmax',20,'plotitc','off');

text(52,20,'Standard');

erspstan=ersp;

figure('Position',[500 scrsz(4)/4 scrsz(3)/2 scrsz(4)/3])

[ersp,itc,powbase,times,freqs]=...

newtimef( mydev,2000,[-500 1500],1000,cyc,'erspmax',20,'plotitc','off');

text(52,20,'Deviant');

erspdiff=ersp-erspstan;

figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/3])

[ersp,itc,powbase,times,freqs]=...

newtimef( myMM,2000,[-500 1500],1000,cyc,'erspmax',20,'plotitc','off');

text(52,20,'Difference');

figure('Position',[500 scrsz(4)/2 scrsz(3)/2 scrsz(4)/3])

imagesc(erspdiff,[-20 20])

set(gca,'YDir','normal')

h = colorbar;

text(52,20,'Subtracted ERSP');

This simulation suggests that it may be preferable when considering difference waveforms to first subtract the waveforms and then conduct time-frequency analysis, as opposed to conducting time-frequency analysis on each waveform and then subtracting the resulting ERSP.