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.