Simplified Introduction, part 3

-->

Understanding more about phase and ITC

Here is the latest addition to the manual for beginners. If you want to begin at the beginning, please scroll down to part 1. I hope to make the whole document available for download in Word format once completed, but if anyone would like a draft copy, please email me: dorothy.bishop@psy.ox.ac.uk.

This week I've learned a lot more about computation of ITC, thanks to James Desjardins from Brock University, who explained how compass plots could be used to see what ITC is measuring. The following incorporates scripts he sent, plus some basic stuff I worked out on how things are computed.

Start with the script for make_mult_trial.m.

Alter the newtimef command at the end of the script so that it gives you one additional piece of output, tfdata.

[ersp,itc,powbase,times,freqs,erspboot,itcboot,tfdata]=newtimef(

mysig',2000,[-500 1500],1000,cyc,'plotphasesign','on');

Now look at the outputs from the function newtimef:

freqs

times

powbase

You will see that freqs gives 24 values ranging from 3.9 to 50 Hz. These are the frequency bins that are used in the time-frequency analysis. You can get a different set of values by specifying a frequency range to analyse (using 'maxfreq' argument) and/or by altering the padratio. The times vector contains 200 time values (in seconds), but the values do not correspond to the values in the input signal: the range is truncated because the moving window used to compute the power at each frequency needs to sample points before and after the time point under consideration, and so must have a smaller range than the original signal.

The powbase variable gives the mean power in the baseline interval: this is subtracted from the power during the post-event period when computing ERSP unless the user specifies 'baseline' NaN.

If you type

size(itc)

size(ersp)

You will see that both variables are matrices of dimension 24 x 200; i.e., they have one value for each frequency and each time point.

Now take a look at tfdata. It is matrix of 24 x 200 x 50 numbers, each with a real and imaginary component. This is the time-frequency decomposition by frequency (24 values) x time points (200 values) x trials (50 values). That's too much data to look at, so we will just select a small portion to explore. The next command selects the values for the 3rd frequency value, the 10th time point and the first five trials:

mytfdata=squeeze(tfdata(3,63,1:5));

(N.B. the squeeze command simply reduces the dimensionality of the matrix by removing any singleton dimensions).

If you looks at freqs you will see that the 3rd frequency is 7.9 Hz, which is nice and close to the 8 Hz frequency that was used to create our sample waveforms. Similarly, inspecting times(63) tells you this corresponds corresponds to a point 177 ms after the baseline, i.e. within the period where phase reset of the signal was specified (150-350 ms).

Now type:

figure;compass(mytfdata);

This gives a figure showing for each trial the amplitude (length of line) and phase (angle of line) at the selected time point. This will look something like this (the precise pattern is not predicatable, because each run of the script has different values for the random component)

As you can see, the phase is similar across all these trials.

Now repeat this exercise using a time period during the baseline:

mytfdata2=squeeze(tfdata(3,10,1:5));

figure;compass(mytfdata2);

You may wonder how the amplitude and phase are extracted from the numbers in tfdata. Let's just work with the first complex number to illustrate:

n=1;

x=real(mytfdata2(n)); %real portion of complex number

y = imag(mytfdata2(n)); % imaginary portion of complex number

z = sqrt(x^2+y^2); % amplitude

p = atan(y/x)*180/pi; % calculated phase in radians, converted to degrees

if p<0 o:p="">

p=p+180;

end;

if y<0>

p=p+180;

end;

z, p

You will see that the length (z) and angle (p) correspond to one of the lines in the figure (precise results vary from run to run). You can repeat for different values of n, and you should see the values of amplitude and phase for each of the five arrows.

ITC is calculated by setting the absolute value (length of the arrow in compass plots) of the complex coefficients to 1 and then calculating the average coefficient. This is illustrated for two different time points in the compass plots generated by the following script:

%Set absolute value of complex coefficients = 1

itcdata=tfdata./abs(tfdata);

figure;compass(itcdata(4,40,:))

hold on;compass(mean(itcdata(4,40,:),3),'r') %red line = average

figure;compass(itcdata(4,64,:))

hold on;compass(mean(itcdata(4,64,:),3),'r')

Notice how the length of the red line corresponds to the values in ITC plots (unambiguously when "plotphasesign" is set to "off"), and that the length of this line is determined only by the degree to which the phase angles cluster together across trials at a give time point and frequency.

Now let's consider how ersp is computed from tfdata. We will just do this for the 3^rd frequency bin, for illlustration

.

mytf3=squeeze(tfdata(3,:,:)); %select just 3^rd frequency bin

mytf3x=real(mytf3).^2+imag(mytf3).^2; %compute power

mytf3m=mean(mytf3x(:,:),2); %mean power across trials per time pt

mye=ersp(3,:); % ersp computed by newtimef for 3^rd freq

myb=powbase(3); % mean power in baseline for 3^rd freq

figure;scatter(10*log10(mytf3m)-myb,mye);

This demonstrates you get the ersp value by taking sum of squares of real and imaginary portions of tfdata for each trial,averaging these across trials, taking 10 x the log value and subtracting the power in the baseline.

If you like, you can compute the mean baseline power directly by adding the 'baseline', NaN input argument to the newtimef command, to generate uncorrected power values. Don't run the whole script (you would generate different values for the random part of the wave) – just run the newtimef command followed by the section of script above to compute mytf3x. You can see from inspecting times that the baseline covers points 1 to 42. Then you get the powbase value with:

mytf3mb=10*log10(mean(mean(mytf3x(1:42,:),2)))