Pageviews from the past week

Saturday 19 November 2011

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 2N 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 5th 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.

Friday 18 November 2011

Simplified Introduction, part 2



-->

Creating multiple trials with varying phase from trial to trial

So far we have used newtimef() to look at a single simulated waveform. To generate more lifelike data, and to explore another main output of newtimef(), inter-trial coherence or ITC, we need to work with multiple trials (or epochs). We need to introduce some random variation from trial to trial in order to produce meaningful results, and so I will first demonstrate a short bit of code to do that.

Creating a random, autocorrelated waveform

The functions rand() or randn() can be used in Matlab to generate random numbers. Here we use randn(), which has the virtue that it is equally likely to generate positive or negative numbers. We could just add a new random number to each point of a sin wave function, but for a more realistic signal, we generated autocorrelated noise, where each point in the function is correlated with the prior point.
%-------------------------------------------------------------------------
% Make_autocorr-noise
% D.V.M. Bishop, 10th July 2010
% Demonstration of how to make an autocorrelated random signal
% N.B. if you have an old version of Matlab you may not have
% autocorr function; you can just delete lines containing that function
% and ACF
%-------------------------------------------------------------------------
for thistrial=1:10;
est_ac=thistrial/10;
% est_ac determines the proportion of the signal that
% is random, and is directly related to the first order autocorrelation
% So this program demonstrates what random signal looks like when
% first-order autocorrelation is approx .1, .2, .3 etc
myrand(1)=.3*randn();%1st point in series is a normal random number from
% distribution with mean 0 and SD = .3;
% Each point on random wave based on previous value +/- random
for myi=2:2000;
myrand(myi)=(1-est_ac)*.3*randn()+est_ac*myrand(myi-1);
end
[ACF,Lags,Bounds]=autocorr(myrand);
% ACF is autocorrelation function; gives correlations of point N with
% points N+1, N+2 etc. The correlation of N and N+1 is first order
% autocorrelation.
% You can see correlation of points N and N+2 by looking at ACF(3),
% and so on.
figure;plot(myrand(:));
axis([0,2000,-3,3])
mytitle=['1st order autocorrelation = ' num2str(ACF(2))]
title(mytitle);
end
%-------------------------------------------------------------------------

Creating multiple trials with noise, and periods of increased power and phase coherence

We will now use a loop to generate complex sine waves as before, adding some random noise to each iteration of the loop.
%-----------------------------------------------------------------------
% make_mult_trial.m
%-----------------------------------------------------------------------
% generate simulated multi-trial epochs with 2 frequencies
close all;
scrsz = get(0,'ScreenSize');
D = 2; %signal duration 2 s
S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave
F = [8 35]; % 2 frequencies in Hz
w = 2*pi*F; % convert frequencies to radians
P = [.15 .55 ]; % 2 corresponding phases
A = [1 .5 ]; % corresponding amplitudes
T = 1/S; % sampling period, i.e. for this e.g. points at 1 ms intervals
t = [T:T:D]; % time vector %corrected on 12th Jul 2010
myphi=2*pi*P; % compute phase angles
nfreqs=length(F); % N frequencies in the complex
ntrial=50;
mysin=zeros(nfreqs+1,ntrial,length(t)); %initialise mysin
% last value for first dimension will hold noise
boostrange=[500 600;800 900];%ms ranges for each freq where signal amplified
phasereset=[650 850;1050 1150];%ms ranges for each freq where signal phase reset
phasenoise=0;ampnoise=0; %initialise values
boostlevel=2; %N times signal multiplied before adding to unboosted signal
noiselevel=.2; %The lower this is, the greater autocorrelation in wave
for thistrial=1:ntrial
% Make random component
mysin(nfreqs+1,thistrial,1)=randn()/10;%1st point is a random number
% Each point on random wave based on previous +/- random
for myi=2:length(t);
mysin(nfreqs+1,thistrial,myi)=noiselevel*randn()/10 ...
+(1-noiselevel)*mysin(nfreqs+1,thistrial,myi-1);
end
%Make systematic component, combination of sine waves
for thisfreq=1:nfreqs
phasenoise=2*pi+1000*randn(); %jitter phase on each trial
ampnoise=.3*randn(); %vary amplitude on each trial
mysin(thisfreq,thistrial,:) = (A(thisfreq)+ampnoise)*(sin(w(thisfreq)*t ...
+ (myphi(thisfreq)+phasenoise)));
%reset phase for relevant interval for each frequency
phasetimebit=[phasereset(thisfreq,1):phasereset(thisfreq,2)];
mysin(thisfreq,thistrial,phasetimebit) = ...
(A(thisfreq)+ampnoise)*(sin(w(thisfreq)*phasetimebit/S + (myphi(thisfreq))));
%multiply signal just for this time range
ampbit=[boostrange(thisfreq,1):boostrange(thisfreq,2)];
mysin(thisfreq,thistrial,ampbit)=(boostlevel*A(thisfreq)+ampnoise)*...
(sin(w(thisfreq)*ampbit/S + (myphi(thisfreq)+phasenoise)));
% uncomment next 2 lines to make signal much smaller in baseline, multiply x .1
% basebit=1:500;
% mysin(thisfreq,thistrial,basebit)=(.1*A(thisfreq)+ampnoise)*(sin(w(thisfreq)*basebit/S + (myphi(thisfreq)+phasenoise)));
end
mysig(thistrial,:)=sum(mysin(:,thistrial,:),1);
end
t=t-.5; %subtract 500 ms to indicate baseline is negative
figure('Position',[1 scrsz(4)/2 scrsz(3)/2 scrsz(4)/3])
plot(t,mysig);
xlabel('Time (seconds)');
ylabel('Amplitude');
title('Individual trials');
figure('Position',[1 50 scrsz(3)/2 scrsz(4)/3]);
plot(t,mean(mysig,1));
title('Grand Average');
cyc=[1]; % cycles
figure('Position',[scrsz(3)/2 scrsz(4)/2-50 scrsz(3)/2 scrsz(4)/2])
[ersp,itc,powbase,times,freqs]=...
newtimef( mysig',2000,[-500 1500],1000,cyc,...
'plotphase','off');
text(52,20,'Complex wave');
%-----------------------------------------------------------------------
Run this script and you will get as output figures showing all individual trials superimposed, and the grand average, as well as the follow time-frequency plot:

Interpreting ERSP and ITC output from multiple trial simulation

The ERSP shows bursts of power in a lower and upper frequency range. These are specified in the program by the variable:
boostrange=[500 600;800 900];
This gives the time range (start and end points) for the boost in amplitude for the first and second frequencies respectively. The time is measured from the start of the trial, and so, given that the baseline is 500 ms, the boosts start at 0 and 300 ms post trial onset (zero time point). You can re-run the program with different values to check you understand how this works. Note that the frequency specificity of the ERSP output is not precise, given that the signal is based on a combination of 8 Hz and 35 Hz sine waves.
The lower portion of the figure shows the inter-trial coherence or ITC. This is shown by default when you use newtimef() or pop_newtimef(), but previously we suppressed this part of output by specifying 'plotitc','off' in the command.
The regions where ITC is red are those that were specified with:
phasereset=[650 850;1050 1150]
In general, on each trial, the phase for each frequency is set to be random. However, for these regions, the script specified that the same phase is used on each trial.
The ITC may be seen as complementary to the ERSP, in that it is insensitive to the power in the signal, but measures the extent to which different trials are in phase. The clearest technical account for beginners I have found is:
Roach, B. J., & Mathalon, D. H. (2008). Event-related EEG time-frequency analysis: An overview of measures and an analysis of early gamma band phase locking in schizophrenia. Schizophrenia Bulletin, 34(5), 907-926.
Essentially, when ITC is high this means that the signal in a given frequency band is in phase on different trials. If ITC is zero it means there is no relationship between phase from one trial to the next. Note that we included in the newtimef() command the setting:
'plotphase','off'
This means that the plot simply shows the magnitude of phase coherence, and not the sign, which is not meaningful (because phase is angular and so the same phase difference could be represented as either a positive or negative number: e.g. a 15 degree phase advance is equivalent to a 345 degree phase delay.
Now change boostlevel to zero; this means you retain the phase resetting but no longer have the periods of boosted signal power. You will see that the ITC is unchanged, but ERSP disappears. Importantly, though, the averaged signal still shows distinct periods of increased amplitude. This illustrates a key point in the debate about mechanisms of ERP generation, first made in this paper:
Sayers, B. M., Beagley, H. A., & Henshall, W. R. (1974). The mechanism of auditory evoked EEG responses. Nature, 247, 481-483.
Sayers et al noted that an ERP would be generated if a signal onset reset the phase of the background EEG, because signal frequencies in phase will add, whereas frequencies which are out of phase will tend to cancel out.

Converting ITC to real numbers

If you inspect the variable itc, you will find each value has a real and imaginary portion. To convert this to the numbers that are plotted in the figure, you have to square each portion, add them, and then take the square root. You can use the following commands to do this and to plot the result:
nuitc=sqrt(real(itc).^2+imag(itc).^2);
figure;imagesc(nuitc,[0 1]);
set(gca,'YDir','normal') % y axis is inverted if you don't specify this
h = colorbar;

Masking non-significant values

When using time-frequency analysis with real data, it can be useful to mask out nonsignificant values of ERSP or ITC, so you can readily visualise the results in terms of regions of significant changes in power or inter-trial coherence. To do this, you use the alpha parameter; this treats all nonsignificant values as zero, and so they will be plotted as green. For instance, in the penultimate line of the program, add:
'alpha', .05
before the final close bracket, and re-run the program. This won't make a very large difference, because with these simulated data, changes to power and coherence are deliberately made very clearcut and dramatic, but with real data, using alpha often makes the output much easier to read.
How do you think the output should change if you set 'alpha' to .001? Try it and see if you are right.
The default method for determining significance levels uses bootstrapping. For instance, for ERSP, subsamples of trials are repeatedly taken from the baseline and the power at different time-points is calculated; this allows one to obtain a distribution of power levels in the baseline, and so to identify when obtained power in the trial period is statistically abnormal (i.e. such an extreme value would occur in only 1/20 baseline samples, if alpha is set to .05). If you want to see the values determined as cutoffs by the bootstrap procedure, you just modify the newtimef() command as follows:
[ersp,itc,powbase,times,freqs, erspboot,itcboot] = ...
For erspboot there will be two values for each frequency band (because ERSP can be abnormally high or abnormally low), whereas for ITC there is just one (because ITC ranges from 0 to 1).

Playing with the script to learn more

A good way to gain an understanding of time-frequency analysis is to play with this script, altering parameters and seeing how this affects the results. If you do so, be aware that the program tries to be helpful by adjusting the scale for ERSP and ITC according to the range in the data. This can actually be unhelpful if you want to see the effect, for instance, of altering the amount of noise in the signal, or the amount of amplitude boost. You can ensure that the colorbar corresponds to the same values on different runs by setting 'erspmax' and 'itcmax'. For instance, if you include:
'itcmax', 1
in the newtimef command, the ITC colorbar will always extend to one.
There are many optional commands with newtimef, which you can see if you type:
help newtimef
Most will not be of use to you, but it is well worth investing time in trying out different options, as this will greatly increase your understanding of how the function works.