ERSP

原文

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.

%-----------------------------------------------------------------------
% modified create_complex.m
%-----------------------------------------------------------------------
D = 2; %signal duration
S = 1000; % sampling rate, i.e. N points pt sec used to represent sine wave
F = [10 20 45]; % 4 frequencies in Hz
w = 2*pi*F; % convert frequencies to radians
P = [0 .5 .25]; % 4 corresponding phases
A = [1 .5 .3]; % 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);
ylim([-2,2])
xlim([0,2])
title('a modified complex wave')
xlabel('Time (seconds)');
ylabel('Amplitude');

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.

newtimef()

[ersp,itc,powbase,times,freqs]=newtimef( mysig,2000,[-500 1500],1000, 0,'plotitc','off');

% log info========================================================
% Computing Event-Related Spectral Perturbation (ERSP) and
% Inter-Trial Phase Coherence (ITC) images based on 1 trials
% of 2000 frames sampled at 1000 Hz.
% Each trial contains samples from -500 ms before to
% 1500 ms after the timelocking event.
% Image frequency direction: normal
% Adjust min freq. to 1.95 Hz to match FFT output frequencies
% Adjust max freq. to 50.78 Hz to match FFT output frequencies
% Using hanning FFT tapering
% Generating 200 time points (-371.9 to 1371.9 ms) !!!!!!
% Finding closest points for time variable
% Time values for time/freq decomposition is not perfectly uniformly distributed
% The window size used is 256 samples (256 ms) wide.
% Estimating 26 linear-spaced frequencies from 2.0 Hz to 50.8 Hz.

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.

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.

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.

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 26 frequencies ranging from 3.09 Hz to 50 Hz. Now look at the variable ersp; you will see it is a 26 x 200 matrix of numbers, i.e. there is a value of ersp for each frequency bin and each time point.

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.

newtimef() raw power

figure;[ersp,itc,powbase,times,freqs]=newtimef( mysig,2000,[-500 1500],1000, 0,'baseline',NaN,'plotitc','off');

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.

newtimef() parameters

[ersp,itc,powbase,times,freqs,erspboot,itcboot] = newtimef(data, frames, epochlim, srate, cycles,'key1',value1, 'key2',value2, ... );

The first argument is the dataset to be analysed. In the context of EEGLAB this will be an EEG dataset.

The second argument is 'frames', i.e. number of points per epoch.

The third argument is the interval over which the power is computed.

The fourth argument in the command is the sampling rate. 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.

To get an analysis of frequency at different time-points, the FFT was applied to the signal over different time windows.

Trade-off between time and frequence resolusion

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.