Exercises Part 3 | Analyzing Neural Time Series Data

Exercises 10 | Convolution

10.6.1

1. Create two kernels for convolution: one that looks like a U and one that looks like a decay function. There is no need to be too sophisticated in generating, for example, a Gaussian and an exponential; numerical approximations are fine.

%% Exercises 10
% 1. Create two kernels 
kernel_U = [1 0.8 0.3 0.1 0 0.1 0.3 0.8 1];
kernel_decay = [1 0.9 0.8 0.7 0.5 0.3 0.2 0.1 0];

10.6.2

1. Convolve these two kernels with 50 time points of EEG data from one electrode. Make a plot showing the kernels, the EEG data, and the result of the convolution between the data and each kernel. Use time-domain convolution as explained in this chapter and as illustrated in the online Matlab code. Based on visual inspection, what is the effect of convolving the EEG data with these two kernels?

%% 2
% Snippet of EEG data
load sampleEEGdata.mat
impEEG = EEG.data(47,100:149,10);

% 用自编函数计算impEEG与kernel_U的convolution
% data that we'll use for convolution (must be zero-padded).
dat4conv = [zeros(1,length(kernel_U)-1) impEEG zeros(1,length(kernel_U)-1) ];

% used for cutting the result of convolution
half_of_kernel_size = ceil((length(kernel_U)-1)/2);

% initialize convolution output
convolution_result_U = zeros(1,length(impEEG)+length(kernel_U)-1);

% run convolution (note that kernel is flipped backwards) ti以外的范围会被舍弃，没必要计算
for ti=half_of_kernel_size+1:length(convolution_result_U)-half_of_kernel_size         
    convolution_result_U(ti) = sum(dat4conv(ti:ti+length(kernel_U)-1).*kernel_U(end:-1:1));
end

% cut off edges
convolution_result_U = convolution_result_U(half_of_kernel_size+1:end-half_of_kernel_size);

% 用Matlab自带函数conv()计算impEEG与kernel_decay的convolution
convolution_result_decay = conv(impEEG,kernel_decay,'same');    %% 'same':conv的结果与impfun的长度相同

% plot kernels
figure
subplot(311)
plot(kernel_U,'bo-','LineWidth',1)
hold on
plot(kernel_decay,'ro-','LineWidth',1)
legend('U','decay','Location',[0.6 0.84 0.13512 0.064683])
xlabel('Points')
title('Convolution kernels')

% plot Snippet of EEG data
subplot(312)
plot(impEEG,'b-','LineWidth',1)
xlabel('Time points')
ylabel('Voltage \muV')
title('Snippet of EEG data (EEG.data(47,100:149,10))')

% plot EEG data before and after convolution
subplot(313)
plot(impEEG,'k-','LineWidth',1)
hold on
plot(convolution_result_U./sum(kernel_U),'b-','LineWidth',1)
plot(convolution_result_decay./sum(kernel_decay),'r-','LineWidth',1)
legend('EEG data','U-convolved','decay-convolved','Location',[0.65 0.245 0.2 0.09])
xlabel('Time points')
ylabel('Voltage \muV')
title('EEG data before and after convolution')

---

Exercises 11 | Fourier Transform

11.12.1

Reproduce the top three panels of figure 11.12 three times.

• First, perform time-domain convolution using the “ manual ” convolution method shown in chapter 10.

• Second, perform frequency-domain convolution using the discrete time Fourier transform presented at the beginning of this chapter.

• Finally, perform frequency-domain convolution using the Matlab functions fft and ifft (do not use the function conv). (You can optionally reproduce the bottom panel of figure 11.12 for the frequency domain analyses; keep in mind that the power scaling is for display purposes only.)

%% Exercises 11
% 1. 
load sampleEEGdata.mat
eegdat4convol = squeeze(EEG.data(47,:,1));
% create Gaussian
time = -1:1/EEG.srate:1;
s = 5/(2*pi*30);
gaussian = exp((-time.^2)/(2*s^2))/30;

%% (1) perform time-domain convolution using the "manual" convolution method
longEEG = [zeros(1,length(gaussian)-1) eegdat4convol zeros(1,length(gaussian)-1)];  % 扩展eegdata左右两边的长度
half_of_kernel_size = ceil((length(gaussian)-1)/2);
conv_result_manual = zeros(1,length(eegdat4convol)+length(gaussian)-1);

for ti = half_of_kernel_size+1 : length(conv_result_manual)-half_of_kernel_size
    conv_result_manual(ti) = sum(longEEG(ti:ti+length(gaussian)-1).*gaussian(end:-1:1));
end
conv_result_manual = conv_result_manual(half_of_kernel_size+1:end-half_of_kernel_size);

%% (2) perform frequency-domain convolution using the discrete time Fourier transform
% 对EEG信号和kernal(Gaussian)进行延长，保证Fourier变换时二者长度一致
padding_length = length(eegdat4convol) + length(gaussian) - 1;
EEG_padded = [eegdat4convol, zeros(1,padding_length-length(eegdat4convol))];
Gaussian_padded = [gaussian, zeros(1,padding_length-length(gaussian))];
fourier_eeg = zeros(1,padding_length);
fourier_gaussian = zeros(1,padding_length);
conv_result_Fourier = zeros(1,padding_length);

time_fourier = (0:padding_length-1)/padding_length;
% Fourier transform
for fi = 1:padding_length
    % create sine wave
    sine_wave  = exp(-1i*2*pi*(fi-1).*time_fourier);
    fourier_eeg(fi) = sum(sine_wave.*EEG_padded);
    fourier_gaussian(fi) = sum(sine_wave.*Gaussian_padded);
end
fourier_result = fourier_eeg .* fourier_gaussian;

% inverse Fourier transform
for fi = 1:padding_length
    sine_wave = exp(1i*2*pi*(fi-1).*time_fourier);
    conv_result_Fourier(fi) = sum(fourier_result.*sine_wave);
end
conv_result_Fourier = real(conv_result_Fourier(half_of_kernel_size+1:end-half_of_kernel_size))./padding_length;

%% (3) perform frequency-domain convolution using the Matlab functions fft and ifft
conv_result_fft = ifft(fft(EEG_padded) .* fft(Gaussian_padded));
conv_result_fft = conv_result_fft(half_of_kernel_size+1:end-half_of_kernel_size);

%% plot
figure
plot(EEG.times,eegdat4convol);
hold on

plot(EEG.times,conv_result_manual,'LineWidth',2);
plot(EEG.times,conv_result_Fourier,'o');
plot(EEG.times,conv_result_fft,'*');
xlabel('Time (ms)')
legend('EEG data','"manual" convolution method','discrete time Fourier transform','Matlab functions fft and ifft')

11.12.2

From the three sets of Matlab code you have for reproducing figure 11.12 , run a computation time test. That is, time how long it takes Matlab to perform 1000 repetitions of each of the three methods for computing convolution that you generated in the previous exercise (do not plot the results each time). You can use the Matlab function pairs tic and toc to time a Matlab process. Plot the results in a bar plot, similar to figure 11.8 .

%% Exercises 11
% 2. 
load sampleEEGdata.mat
eegdat4convol = squeeze(EEG.data(47,:,1));
% create Gaussian
time = -1:1/EEG.srate:1;
s = 5/(2*pi*30);
gaussian = exp((-time.^2)/(2*s^2))/30;
t = zeros(1,3);
%% (1) perform time-domain convolution using the "manual" convolution method
longEEG = [zeros(1,length(gaussian)-1) eegdat4convol zeros(1,length(gaussian)-1)];  % 扩展eegdata左右两边的长度
half_of_kernel_size = ceil((length(gaussian)-1)/2);
conv_result_manual = zeros(1,length(eegdat4convol)+length(gaussian)-1);

% repeat 1000 times
tic;
for i = 1:1000
    for ti = half_of_kernel_size+1 : length(conv_result_manual)-half_of_kernel_size
        conv_result_manual(ti) = sum(longEEG(ti:ti+length(gaussian)-1).*gaussian(end:-1:1));
    end
    conv_result_manual = conv_result_manual(half_of_kernel_size+1:end-half_of_kernel_size);
end
t(1) = toc;
%% (2) perform frequency-domain convolution using the discrete time Fourier transform
% 对EEG信号和kernal(Gaussian)进行延长，保证Fourier变换时二者长度一致
padding_length = length(eegdat4convol) + length(gaussian) - 1;
EEG_padded = [eegdat4convol, zeros(1,padding_length-length(eegdat4convol))];
Gaussian_padded = [gaussian, zeros(1,padding_length-length(gaussian))];
fourier_eeg = zeros(1,padding_length);
fourier_gaussian = zeros(1,padding_length);
conv_result_Fourier = zeros(1,padding_length);

time_fourier = (0:padding_length-1)/padding_length;

% repeat 1000 times
tic
for i = 1:1000
    % Fourier transform
    for fi = 1:padding_length
        % create sine wave
        sine_wave  = exp(-1i*2*pi*(fi-1).*time_fourier);
        fourier_eeg(fi) = sum(sine_wave.*EEG_padded);
        fourier_gaussian(fi) = sum(sine_wave.*Gaussian_padded);
    end
    fourier_result = fourier_eeg .* fourier_gaussian;
    
    % inverse Fourier transform
    for fi = 1:padding_length
        sine_wave = exp(1i*2*pi*(fi-1).*time_fourier);
        conv_result_Fourier(fi) = sum(fourier_result.*sine_wave);
    end
    conv_result_Fourier = real(conv_result_Fourier(half_of_kernel_size+1:end-half_of_kernel_size))./padding_length;
end
t(2) = toc;
%% (3) perform frequency-domain convolution using the Matlab functions fft and ifft
% repeat 1000 times
tic
for i = 1:1000
    conv_result_fft = ifft(fft(EEG_padded) .* fft(Gaussian_padded));
    conv_result_fft = conv_result_fft(half_of_kernel_size+1:end-half_of_kernel_size);
end
t(3) = toc;
%% plot results in a bar plot
b = bar(t);
labels = string(b(1).YData);
text(b(1).XEndPoints,b(1).YEndPoints,labels,'HorizontalAlignment','center','VerticalAlignment','top','Color',[1 1 1])
set(gca,'xticklabel',{'manual convolution';'DTFT&IDTFT';'FFT&IFFt'})
set(gca, 'YScale', 'log');
ylabel('time(s)')

11.12.3

Generate a time series by creating and summing sine waves, as in figure 11.2B . Use between two and four sine waves, so that the individual sine waves are still somewhat visible in the sum. Perform a Fourier analysis (you can use the fft function) on the resulting time series and plot the power structure. Confirm that your code is correct by comparing the frequencies with nonzero power to the frequencies of the sine waves that you generated. Now try adding random noise to the signal before computing the Fourier transform. First, add a small amount of noise so that the sine waves are still visually recognizable. Next, add a large amount of noise so that the sine waves are no longer visually recognizable in the time domain data. Perform a Fourier analysis on the two noisy signals and plot the results. What is the effect of a small and a large amount of noise in the power spectrum? Are the sine waves with noise easier to detect in the time domain or in the frequency domain, or is it equally easy/difficult to detect a sine wave in the presence of noise?

%% Exercises 11

% Create sine waves
srate = 400;
time=-1:1/srate:1;
N = length(time);

% create three sine waves
s1 = 20*sin(2*pi*3*time);
s2 = 10*sin(2*pi*8*time);
s3 = 5*sin(2*pi*15*time);
sum_Sine = s1+s2+s3;

% add random noise
Sine_LitteNoise = sum_Sine + 10*randn(1,N);
Sine_LargeNoise = sum_Sine + 40*randn(1,N);

f = linspace(0,srate/2,floor(N/2)+1);

FFT_Sine = fft(sum_Sine)./(N/2);
FFT_LitteNoise = fft(Sine_LitteNoise)./(N/2);
FFT_LargeNoise = fft(Sine_LargeNoise)./(N/2);

FFT_Sine = abs(FFT_Sine(1:ceil(N/2)));
FFT_LitteNoise = abs(FFT_LitteNoise(1:ceil(N/2)));
FFT_LargeNoise = abs(FFT_LargeNoise(1:ceil(N/2)));

%% Plot
figure
subplot(4,1,1);
plot(time,sum_Sine);
xlabel('Time(s)')
ylabel('amplitude')
title('Original time series')

subplot(4,1,2)
bar(f,FFT_Sine);
set(gca,'xlim',[0,20],'ylim',[0,25])
xlabel('Frequency(Hz)')
ylabel('amplitude')
title('Power spectrum')

subplot(4,2,5)
plot(time,Sine_LitteNoise)
xlabel('Time(s)')
ylabel('amplitude')
title('Time series, small noise')

subplot(4,2,7)
bar(f,FFT_LitteNoise)
set(gca,'xlim',[0,20],'ylim',[0,25])
xlabel('Frequency(Hz)')
ylabel('amplitude')
title('Power spectrum, small noise')

subplot(4,2,6)
plot(time,Sine_LargeNoise)
xlabel('Time(s)')
ylabel('amplitude')
title('Time series, large noise')

subplot(4,2,8)
bar(f,FFT_LargeNoise)
set(gca,'xlim',[0,20],'ylim',[0,25])
xlabel('Frequency(Hz)')
ylabel('amplitude')
title('Power spectrum, large noise')

---

Exercises 12 | Morlet wavelets

1. Create a family of Morlet wavelets ranging in frequency from 2 Hz to 30 Hz in five steps.

2. Select one electrode from the scalp EEG dataset and convolve each wavelet with EEG data from all trials from that electrode. Apply the Matlab function real to the convolution result, as in convol_result=real(convol_result). This will return the EEG data bandpass filtered at the peak frequency of the wavelet. You learn more about why this is in the next chapter.

3. Average the result of convolution over all trials and plot an ERP corresponding to each wavelet frequency. Each frequency should be in its own subplot.

4. Plot the broadband ERP (without any convolution). Thus, you will have six subplots in one figure. How do the wavelet-convolved ERPs compare with the broadband ERP? Are there dynamics revealed by the wavelet-convolved ERPs that are not apparent in the broadband ERP, and are there dynamics in the broadband ERP that are not apparent in the waveletconvolved ERPs? Base your answer on qualitative visual inspection of the results; statistics or other quantitative comparisons are not necessary.

%% Exercise 12
%% 1. Create a family of Morlet wavelets
load sampleEEGdata.mat
frequencies = linspace(2,30,5);
num_wavelets = length(frequencies);

srate = EEG.srate;
time = -1:1/srate:1;
wavelet_family = zeros(num_wavelets,length(time));

% Loop through frequencies and make a family of wavelets.
for fi=1:num_wavelets
 
    % create a sine wave at this frequency
    sinewave = exp(2*1i*pi*frequencies(fi).*time); % the "1i" makes it a complex wavelet
    
    % create a Gaussian window
    gaus_win = exp(-time.^2./(2*(6/(2*pi*frequencies(fi)))^2));
    
    % create wavelet via element-by-element multiplication of the sinewave and gaussian window
    wavelet_family(fi,:) = sinewave.*gaus_win;

end

%% 2. Get EEG data from electrode 47, convolve each wavelet with EEG data
% 
channel = 47;
n_conv = EEG.pnts + length(time) - 1;
convol_result = zeros(EEG.trials,num_wavelets,n_conv);
half_of_kernal_size = ceil((length(time)-1)/2);

for trial_i = 1:EEG.trials
    fft_EEG = fft(EEG.data(channel,:,trial_i),n_conv);
    for fi = 1:num_wavelets
        fft_wavelet = fft(wavelet_family(fi,:),n_conv);
        convol_result(trial_i,fi,:) = ifft(fft_EEG.*fft_wavelet,n_conv) * sqrt(6/(2*pi*frequencies(fi))) / 10;
    end
end
convol_result = real(convol_result(:,:,half_of_kernal_size+1:n_conv-half_of_kernal_size));
filtered_ERP = squeeze(mean(convol_result,1));		% 经过wavelet滤波后的ERP
broadband_ERP = squeeze(mean(EEG.data(channel,:,:),3));		 % 未滤波的ERP

%% plot
figure
for i = 1:num_wavelets
    subplot(2,3,i)
    plot(EEG.times,filtered_ERP(i,:));
    set(gca,'xlim',[-300, 1000])
    ylim("padded")
    xlabel('Time(ms)')
    title(['Wavelet filtered ERP at ' num2str(frequencies(i)) ' Hz'])
end

subplot(2,3,6)
plot(EEG.times,broadband_ERP);
ylim("padded")
xlabel('Time(ms)')
title('Unfiltered broadband ERP')

---

Exercises 13 | complex Morlet wavelets

1. Create a family of complex Morlet wavelets

Create a family of complex Morlet wavelets ranging in frequencies from 2 Hz to 30 Hz in five steps.

%% Exercises 13
% 1. Create a family of complex Morlet wavelets
load sampleEEGdata.mat
frex = linspace(2,30,5);
time = -2:1/EEG.srate:2;

% 在3和10之间生成5个对数分别的数，作为不同频率下wavelet包含的周期数目
numcycles = logspace(log10(3),log10(10),length(frex));
wavelet = zeros(length(frex), length(time));

% make complex Morlet wavelets
for fi = 1:length(frex)
    s = numcycles(fi)/(2*pi*frex(fi));
    wavelet(fi,:) = 1/sqrt(s*sqrt(pi)) * exp(-time.^2./(2*s^2)) .* exp(1i*2*pi*frex(fi).*time);
    
    % 检查两侧是否衰减到0
    % subplot(ceil(length(frex)/3), ceil(length(frex)/2), fi);
    % plot(time,real(wavelet));
end

2. Convolve each wavelet with EEG data

Convolve each wavelet with EEG data from all electrodes and from one trial.

3. Extract power and phase

Extract power and phase from the result of complex wavelet convolution and store in a time × frequency × electrodes × power/phase matrix (thus, a 640 × 5 × 64 × 2 matrix).

%% 2. Convolve each wavelet with EEG data from all electrode and from one trial
trial = 1;

% FFT parameters
n_wavelet            = length(time);
n_data               = EEG.pnts;
n_convolution        = n_wavelet+n_data-1;
half_of_wavelet_size = (length(time)-1)/2;

% 3. a time × frequency × electrodes × power/phase matrix (thus, a 640 × 5 × 64 × 2 matrix)
convolution_result_fft = zeros(n_convolution,length(frex),EEG.nbchan,2);
% FFT of wavelet and EEG data
for electrodei = 1:EEG.nbchan
    fft_data    = fft(squeeze(EEG.data(electrodei,:,trial)),n_convolution); % trial 1
    for fi = 1:length(frex)
        fft_wavelet = fft(wavelet(fi,:),n_convolution);
        convoli = ifft(fft_wavelet.*fft_data,n_convolution) * sqrt(s);
        convolution_result_fft(:,fi,electrodei,1) = abs(convoli.^2);
        convolution_result_fft(:,fi,electrodei,2) = angle(convoli);
    end
end

% cut off edges
convolution_result_fft = convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size,:,:,:);

4. Make topographical plots of power and phase

Make topographical plots of power and phase at 180 ms at all frequencies (hint: you may need to use the squeeze function to remove singleton dimensions). Arrange the plots in one figure with five columns for frequency and two rows for power/phase. Put labels in the plot so it is clear which topographical maps correspond to which frequencies.

%% make topographical features in power ans phase at 180 ms
figure
t = 180;
for fi = 1:length(frex)
    subplot(2,5,fi)
    topoplot(squeeze(convolution_result_fft(EEG.times>=t-500/EEG.srate & EEG.times <= t+500/EEG.srate,fi,:,1)),EEG.chanlocs,'electrodes','off','numcontour',0);
    title(['Pwr:' num2str(frex(fi)) 'Hz,' num2str(t) 'ms'])

    subplot(2,5,fi+5)
    topoplot(squeeze(convolution_result_fft(EEG.times>=t-500/EEG.srate & EEG.times <= t+500/EEG.srate,fi,:,2)),EEG.chanlocs,'electrodes','off','numcontour',0);
    title(['Phs:' num2str(frex(fi)) 'Hz,' num2str(t) 'ms'])
end

5. Repeat step 4 for activity at 360 ms.

%% make topographical features in power ans phase at 360 ms
figure
t = 360;
for fi = 1:length(frex)
    subplot(2,5,fi)
    topoplot(squeeze(convolution_result_fft(EEG.times>=t-500/EEG.srate & EEG.times <= t+500/EEG.srate,fi,:,1)),EEG.chanlocs,'electrodes','off','numcontour',0);
    title(['Pwr:' num2str(frex(fi)) 'Hz,' num2str(t) 'ms'])

    subplot(2,5,fi+5)
    topoplot(squeeze(convolution_result_fft(EEG.times>=t-500/EEG.srate & EEG.times <= t+500/EEG.srate,fi,:,2)),EEG.chanlocs,'electrodes','off','numcontour',0);
    title(['Phs:' num2str(frex(fi)) 'Hz,' num2str(t) 'ms'])
end

6. Think

Are there any prominent topographical features in power or in phase? - Yes.

Do these differ for different frequencies? - Yes.

Do power and phase have similar topographical distributions? - It doesn’t look very similar.

Is there any reason to suspect that they might have similar or different topographies?

(From ChatGPT for reference, GPT is not provided with corresponding pictures)

When analyzing EEG data, particularly through time-frequency decomposition, you might observe various topographical features in both power and phase distributions across the scalp. These features can provide insights into underlying neural processes. Let’s break down each question:

(1) Are there any prominent topographical features in power or in phase?

• Power Topography: Power topography represents the distribution of signal power across the scalp. Prominent features in power maps could indicate areas of strong neural activity, which might correspond to specific cognitive or sensory processes. For example, increased power in the occipital regions might be related to visual processing.
• Phase Topography: Phase topography shows the distribution of phase angles across the scalp. Prominent features in phase topography could indicate phase synchronization between regions, which might be associated with communication between different brain areas.

(2) Do these differ for different frequencies?

• Yes, topographical features often differ for different frequencies. Different frequencies in EEG are associated with different types of neural activity (e.g., alpha waves in the 8-12 Hz range are often related to relaxation and closed eyes, while beta waves in the 13-30 Hz range are associated with active thinking or focus). The power and phase topographies can change significantly across different frequency bands, reflecting the different functional roles of these frequencies.

(3) Do power and phase have similar topographical distributions?

• Power and phase distributions are often not identical. Power topography might highlight regions with strong amplitude, while phase topography might emphasize regions where the timing of oscillatory activity is synchronized across electrodes. However, there can be similarities, especially in regions where strong power is accompanied by consistent phase alignment (e.g., in tasks requiring coordinated neural activity across brain regions).

(4) Is there any reason to suspect that they might have similar or different topographies?

• Similar Topographies: Power and phase could have similar topographies in cases where high power is associated with phase-locked activity across trials or regions. For instance, during sensory processing, a stimulus might evoke both strong power and phase synchronization in relevant brain areas.
• Different Topographies: Conversely, power and phase might differ in their topographical distributions due to different underlying mechanisms. High power might occur in one region without strong phase alignment, especially if the neural activity is more localized or independent across trials. Additionally, phase synchronization can occur without strong power, especially in low-amplitude oscillations that are critical for certain cognitive processes.

Summary:

• Power and phase topographies can both reveal important aspects of neural activity, but they often provide different insights. Power maps tend to show where neural activity is strong, while phase maps can indicate where different regions are temporally synchronized. The topographies of power and phase may differ depending on the frequency of interest and the specific neural processes involved. Understanding these distributions and their relationships is key to interpreting EEG data in the context of neural dynamics.

7. Create a circular colormap

Because phase values are circular ( – π and + π are identical), most color maps are inappropriate because they suggest that – π and + π are very different values (represented, e.g., by blue and red colors). Create a circular colormap that can be used for phase values. You can do this by setting the red, green, and/or blue values to be a cosine function rather than a linear function. Recreate the phase topographical maps. Do they look any different with the new color maps?

%% 自定义用于相位的环形colormap，使得-pi和pi对应的颜色相同
figure

for fi = 1:length(frex)
    subplot(2,3,fi)
    topoplot(squeeze(convolution_result_fft(EEG.times>=t-500/EEG.srate & EEG.times <= t+500/EEG.srate,fi,:,2)),EEG.chanlocs,'electrodes','off','numcontour',0);
    title(['Phs:' num2str(frex(fi)) 'Hz,' num2str(t) 'ms'])
end

ax6 = subplot(2,3,6);

nColors = 128; % 调色板的颜色数量
cmap = hsv(nColors); % 使用 HSV colormap
cmap = cmap([1:end 1], :); % 将第一个颜色重复，使得 -pi 和 pi 一致
colormap(cmap);

% 添加环形 colorbar
cb = colorbar;
set(cb, 'Ticks', [-pi, -pi/2, 0, pi/2, pi], 'TickLabels', {'-\pi', '-\pi/2', '0', '\pi/2', '\pi'});
clim([-pi pi]);

axis off
cb.Position(1) = ax6.Position(1);
cb.Position(2) = ax6.Position(2);

---

Exercises 14 | filter-Hilbert

1. Pick two frequencies (e.g., 5 Hz and 25 Hz) and one electrode and perform complex Morlet wavelet convolution and filter-Hilbert using those two frequencies as the peak/center frequencies for all trials. Plot the resulting power and the bandpass-filtered signal (that is, the real component of the analytic signal) from each method. Plot one single trial (you can choose the trial randomly but plot the same trial for both methods) and then plot the average of all trials. Describe some similarities and differences between the results of the two time-frequency decomposition methods.

%% Exercises 14
%% 1. perform complex Morlet wavelet convolution and filter-Hilbert using frequencies 5 Hz and 25 Hz
load sampleEEGdata.mat
% two frequencies
frex = [5, 25]; % the center frequencies
% pick one electrode (Fcz) and one trial
electrode = find(strcmpi({EEG.chanlocs.labels},'Fcz'));
trial = 1;

data2filter = squeeze(EEG.data(electrode,:,trial));
data2filter_average = squeeze(mean(EEG.data(electrode,:,:),3));
%% perform complex Morlet wavelet convolution
% wavelet parameters
time = -1:1/EEG.srate:1;
numcycles = [4 5];
wavelet = zeros(1,length(time));

% FFT parameters
n_wavelet            = length(time);
n_data               = EEG.pnts;
n_convolution        = n_wavelet+n_data-1;
half_of_wavelet_size = (length(time)-1)/2;

conv_result_wavelet = zeros(length(frex),n_convolution);
conv_result_wavelet_2 = zeros(length(frex),n_convolution);
fft_data = fft(data2filter,n_convolution);
fft_data_2 = fft(data2filter_average,n_convolution);
%
for fi = 1:length(frex)
    s = numcycles(fi)/(2*pi*frex(fi));
    wavelet = exp(-time.^2./(2*s^2)) .* exp(1i*2*pi*frex(fi).*time);
    
    fft_wavelet = fft(wavelet,n_convolution);
    fft_wavelet = fft_wavelet ./ max(abs(fft_wavelet));
    conv_result_wavelet(fi,:) = ifft(fft_wavelet.*fft_data,n_convolution)*(1/sqrt(s*sqrt(pi)));
    conv_result_wavelet_2(fi,:) = ifft(fft_wavelet.*fft_data_2,n_convolution)*(1/sqrt(s*sqrt(pi)));
end

% cut off edges
conv_result_wavelet = conv_result_wavelet(:,half_of_wavelet_size+1:end-half_of_wavelet_size);
conv_result_wavelet_2 = conv_result_wavelet_2(:,half_of_wavelet_size+1:end-half_of_wavelet_size);
%% perform filter-Hilber
% parameters of filter
freqspread = [2 10]; % Hz +/- the center frequency
transwid   = [.01 .15];
nyquist = [1;1] .* EEG.srate/2;
% construct filter kernels
ffrequencies  = [ zeros(2,1) ([1 1]-transwid)'.*(frex-freqspread)' (frex-freqspread)' (frex+freqspread)' ([1 1]+transwid)'.*(frex+freqspread)' nyquist ]./nyquist(1);
idealresponse = [ 0 0 1 1 0 0 ];
filter_result = zeros(length(fi),EEG.pnts);
filter_result_2 = zeros(length(fi),EEG.pnts);
for fi = 1:length(frex)
    filterweights = firls(100,ffrequencies(fi,:),idealresponse);
    fft_filter = fft(filterweights);
    filter_result(fi,:) = hilbert(filtfilt(filterweights,1,data2filter));
    filter_result_2(fi,:) = hilbert(filtfilt(filterweights,1,data2filter_average));
end

figure
for i = 1:length(frex)
    subplot(2,1,i)
    plot(EEG.times,abs(conv_result_wavelet(i,:)).^2);
    hold on
    plot(EEG.times,abs(filter_result(i,:)).^2)
    set(gca,"XLim",[-500,1000])
    xlabel('Time (ms)')
    ylabel('Voltage (\muV)')
    title(['Center Frequency: ' num2str(frex(i)) ' Hz'])
    legend('wavelet','filter-Hilbert')
end

figure
for i = 1:length(frex)
    subplot(2,1,i)
    plot(EEG.times,abs(conv_result_wavelet_2(i,:)).^2);
    hold on
    plot(EEG.times,abs(filter_result_2(i,:)).^2)
    set(gca,"XLim",[-500,1000])
    xlabel('Time (ms)')
    ylabel('Voltage (\muV)')
    title(['Center Frequency: ' num2str(frex(i)) ' Hz'])
    legend('wavelet','filter-Hilbert')
end

1. Modify the wavelet and filter settings (but keep the peak/center frequencies the same) until these two methods produce very similar results. Next, modify the wavelet and filter settings (except the peak/center frequencies) to make the results different (stay within a reasonable range of parameter settings; they do not need to look dramatically different). Which parameters did you change to make the results look more similar versus more different? How different are the results, and would you consider this a meaningful difference? What does this difference tell you about when to use specific parameter settings for wavelet convolution and the filter-Hilbert method?

尝试改变的参数：

• wavelet：
  • numcycles：wavelet包含的周期数目参数。包含的周期数目越多，wavelet频带越窄，频率分辨率越高
• filter-Hilbert：
  • freqspread：带宽参数。带宽越大，频率分辨率越低
  • transwid：过渡区宽度参数。
  • firls()的第一个参数，the order of the filter：决定filter kernel频率响应的精度。

---

Exercises 15 | The Short-Time FFT

1. Compute the short-time FFT at each electrode and make topographical maps of theta band (around 6 Hz) power and alpha-band (around 10 Hz) power at 150 ms and 700 ms.

%% Compute the short-time FFT at each electrode
% and make topographical maps of theta band (around 6 Hz) power and alpha-band (around 10 Hz) power at 150 ms and 700 ms. 
load sampleEEGdata.mat

timewin = 400;      % in ms
times2save = -300:50:1000; % in ms
timepoint(1) = 150;   % ms
timepoint(2) = 700;   % ms
frequency(1) = 6;     % Hz
frequency(2) = 10;    % Hz

% convert ms to idx
times2saveidx = zeros(size(times2save));
for frexi=1:length(times2save)
    [junk,times2saveidx(frexi)]=min(abs(EEG.times-times2save(frexi)));
end
timewinidx = round(EEG.srate*timewin/1000);

% create hann taper
hann_win = 0.5*(1-cos(2*pi*(0:timewinidx-1)/(timewinidx-1)));

% define frequencies
frex = linspace(0,EEG.srate/2,floor(timewinidx/2)+1);

% initialize power output matrix
tf = zeros(EEG.nbchan,length(frex),length(times2save));

% loop over time points and perform FFT
for channi = 1:EEG.nbchan
    for timepointi=1:length(times2save)
        
        % extract time series data for this center time point
        tempdat = squeeze(EEG.data(channi,times2saveidx(timepointi)-floor(timewinidx/2):times2saveidx(timepointi)+floor(timewinidx/2)-mod(timewinidx+1,2),:)); % note: the 'mod' function here corrects for even or odd number of points
        
        % taper data (using bsxfun instead of repmat... note sizes of tempdat
        % and hann_win)
        taperdat = bsxfun(@times,tempdat,hann_win');
        
        fdat = fft(taperdat,[],1)/timewinidx; % 3rd input is to make sure fft is over time
        tf(channi,:,timepointi) = mean(abs(fdat(1:floor(timewinidx/2)+1,:)).^2,2); % average over trials
    end
end

[junk,frequencyidx(1)]=min(abs(frex-frequency(1)));
[junk,frequencyidx(2)]=min(abs(frex-frequency(2)));
[junk,timepointidx(1)]=min(abs(times2save-timepoint(1)));
[junk,timepointidx(2)]=min(abs(times2save-timepoint(2)));

figure
for timei = 1:2
    for frexi = 1:2
        subplot(2,2,(timei-1)*2+frexi)
        topoplot(squeeze(tf(:,frequencyidx(frexi),timepointidx(timei))),EEG.chanlocs);
        title([num2str(timepoint(timei)) 'ms, ' num2str(frequency(frexi)) 'Hz'])
    end
end

1. Select one electrode and one frequency and compute power over time at that electrode and that frequency using complex wavelet convolution, filter-Hilbert, and the short-time FFT. Plot the results of these three time-frequency decomposition methods in different subplots of one figure. Note that the scaling might be different because no baseline normalization has been applied. How visually similar are the results from these three methods? If the results from the three methods are different, how are they different, and what parameters do you think you could change in the three methods to make the results look more or less similar?

%% Select one electrode and one frequency and compute power over time at that electrode and that frequency
% using complex wavelet convolution, filter-Hilbert, and the short-time FFT.

load sampleEEGdata.mat
channel2plot = 'Fcz';
electrode = strcmpi(channel2plot,{EEG.chanlocs.labels});
frequency = 6;     % in Hz

times2save = -200:1000/EEG.srate:800;
% convert from ms to index
times2saveidx = zeros(size(times2save));
for i=1:length(times2save)
    [junk,times2saveidx(i)]=min(abs(EEG.times-times2save(i)));
end

%% complex wavelet convolution
time = -1:1/EEG.srate:1;
numcycle = 4;
s = numcycle/(2*pi*frequency);
wavelet = exp(-time.^2./(2*s^2)) .* exp(1i*2*pi*frequency.*time);

% FFT parameters
n_wavelet            = length(time);
n_data               = EEG.pnts;
n_convolution        = n_wavelet+n_data-1;
half_of_wavelet_size = (length(time)-1)/2;

% FFT of wavelet and EEG data
fft_wavelet = fft(wavelet,n_convolution);
convol = zeros(EEG.trials,n_convolution);
for triali = 1:EEG.trials
    fft_data    = fft(squeeze(EEG.data(electrode,:,triali)),n_convolution);
    convol(triali,:) = ifft(fft_wavelet.*fft_data,n_convolution);
end

convolution_result_fft = mean((abs(convol).^2),1);

% cut off edges
convolution_result_fft = convolution_result_fft(half_of_wavelet_size+1:end-half_of_wavelet_size);

%% filter-Hilbert
freqspread = 3; % Hz +/- the center frequency
transwid   = 0.01;
nyquist = EEG.srate/2;
% construct filter kernels
ffrequencies  = [ 0 (1-transwid)'.*(frequency-freqspread)' (frequency-freqspread)' (frequency+freqspread)' (1+transwid)'.*(frequency+freqspread)' nyquist ]./nyquist;
idealresponse = [ 0 0 1 1 0 0 ];

filterweights = firls(200,ffrequencies,idealresponse);
filter_result = zeros(EEG.trials,EEG.pnts);

for triali = 1:EEG.trials
    filter_result(triali,:) = hilbert(filtfilt(filterweights,1,EEG.data(electrode,:,triali)));
end
filter_result = mean(abs(filter_result.^2),1);
%% the short-time FFT
timewin = 400;      % in ms
timewinidx = round(EEG.srate*timewin/1000);

% create hann taper
hann_win = .5*(1-cos(2*pi*(0:timewinidx-1)/(timewinidx-1)));

% define frequencies
frex = linspace(0,EEG.srate/2,floor(timewinidx/2)+1);

% initialize power output matrix
tf = zeros(length(frex),length(times2save));

% loop over time points and perform FFT
for timepointi=1:length(times2save)
    
    % extract time series data for this center time point
    tempdat = squeeze(EEG.data(electrode,times2saveidx(timepointi)-floor(timewinidx/2):times2saveidx(timepointi)+floor(timewinidx/2)-mod(timewinidx+1,2),:)); % note: the 'mod' function here corrects for even or odd number of points
    
    % taper data (using bsxfun instead of repmat... note sizes of tempdat
    % and hann_win)
    taperdat = bsxfun(@times,tempdat,hann_win');
    
    fdat = fft(taperdat,[],1)/timewinidx; % 3rd input is to make sure fft is over time
    tf(:,timepointi) = mean(abs(fdat(1:floor(timewinidx/2)+1,:)).^2,2); % average over trials
end
[junk,frequencyidx]=min(abs(frex-frequency));
%% plot
figure
plot(times2save,tf(frequencyidx,:)./(max(tf(frequencyidx,:))),'b','LineWidth',1)
hold on
plot(times2save,convolution_result_fft(times2saveidx)./(max(convolution_result_fft(times2saveidx))),'Color',[0 0.8 0],'LineWidth',1)
plot(times2save,filter_result(times2saveidx)./(max(filter_result(times2saveidx))),'r','LineWidth',1)
legend('Short-Time FFT','complex wavelet','Hilbert-filter')
set(gca,'ylim',[0.2,1.3]) 
colormap parula

---

Exercises 16 | Multitaper

1. Pick one electrode and compute a time-frequency map of power using both the multitaper method and the short-time FFT. Store all of the power values for all of the trials.

2. Next, compute a time-frequency map of signal-to-noise ratio. The signal-to-noise ratio of power is discussed more in chapter 18, but it can be estimated as the average power at each time-frequency point across trials, divided by the standard deviation of power at each time-frequency point across trials.

3. Make time-frequency plots of power and signal-to-noise ratio from the two methods. Make another plot in which you subtract the signal-to-noise plots between the two methods. Are there any noticeable differences between the signal-to-noise results when the multitaper method versus the short-time FFT is used?

%% Pick one electrode and compute a time-frequency map of power using both the multitaper method and the short-time FFT.
load sampleEEGdata

channel2plot = 'p7';
channel2plotidx = strcmpi(channel2plot,{EEG.chanlocs.labels});

timewin = 400;      % in ms
times2save = -300:50:1000; % in ms
% convert from ms to index
timewinidx = round(EEG.srate*timewin/1000);
times2saveidx = zeros(size(times2save));
for i = 1:length(times2save)
    [~,times2saveidx(i)] = min(abs(EEG.times-times2save(i)));
end

baseline_range  = [-200 -00];
% find baselinetimepoints
baseidx = zeros(size(baseline_range));
[~,baseidx(1)] = min(abs(times2save-baseline_range(1)));
[~,baseidx(2)] = min(abs(times2save-baseline_range(2)));
% define frequencies
frex = linspace(0,EEG.srate/2,floor(timewinidx/2)+1);

%% multitaper

% 参数设置
N = timewinidx; % 样本点数
W = 4;    % 频带宽度
Fs = EEG.srate; % 采样率
NW = N * (W / (Fs/2));

% 生成tapers
[tapers, eigenvalues] = dpss(N, NW);

% 选择特征值大于0.9的tapers
tapersToUse = tapers(:, eigenvalues > 0.9);

% initialize output matrix
multitaper_tf = zeros(length(frex),EEG.trials,length(times2save));

% loop through time bins
for timepointi=1:length(times2save)  % 对应一个时间点
    
    % initialize power vector (over tapers)
    % 一个时间点处各频率的功率值，对应tf图上与时间轴垂直的一条线段
    taperpow = zeros(length(frex),EEG.trials);
    
    % loop through tapers
    for tapi = 1:size(tapersToUse,2)
        % extract time series data for this center time point
        tempdat = squeeze(EEG.data(channel2plotidx,times2saveidx(timepointi)-floor(timewinidx/2):times2saveidx(timepointi)+floor(timewinidx/2)-mod(timewinidx+1,2),:)); % note: the 'mod' function here corrects for even or odd number of points
        data      = bsxfun(@times,tempdat,tapersToUse(:,tapi));

        % timewindow截取、taper两端减弱后的数据进行Fourier变换，提取特定频率的功率值
        pow       = fft(data,timewinidx)/timewinidx;
        pow       = pow(1:length(frex),:);
        taperpow  = taperpow + pow.*conj(pow);
    end
    
    % finally, get power from closest frequency
    multitaper_tf(:,:,timepointi) = taperpow/tapi;
end

%% the short-time FFT
% create hann taper
hann_win = .5*(1-cos(2*pi*(0:timewinidx-1)/(timewinidx-1)));

% initialize power output matrix
fft_tf = zeros(length(frex),EEG.trials,length(times2save));

% loop over time points and perform FFT
for timepointi=1:length(times2save)
    
    % extract time series data for this center time point
    tempdat = squeeze(EEG.data(channel2plotidx,times2saveidx(timepointi)-floor(timewinidx/2):times2saveidx(timepointi)+floor(timewinidx/2)-mod(timewinidx+1,2),:)); % note: the 'mod' function here corrects for even or odd number of points
    
    % taper data (using bsxfun instead of repmat... note sizes of tempdat
    % and hann_win)
    taperdat = bsxfun(@times,tempdat,hann_win');
    
    fdat = fft(taperdat,[],1)/timewinidx; % 3rd input is to make sure fft is over time
    fft_tf(:,:,timepointi) =abs(fdat(1:floor(timewinidx/2)+1,:)).^2; % average over trials
end

%% compute a time-frequency map of signal-to-noise ratio.
multitaper_mean = squeeze(mean(multitaper_tf,2));
fft_mean = squeeze(mean(fft_tf,2));

multitaper_std = squeeze(std(multitaper_tf,[],2));
fft_std = squeeze(std(fft_tf,[],2));

snr_multitaper = multitaper_mean ./ multitaper_std;
snr_fft = fft_mean ./ fft_std;

%% plot
% db-correct
db_multitaper_tf = 10*log10( multitaper_mean ./ repmat(mean(multitaper_mean(:,baseidx(1):baseidx(2)),2),1,length(times2save)) );
% db-correct
db_fft_tf = 10*log10( fft_mean ./ repmat(mean(fft_mean(:,baseidx(1):baseidx(2)),2),1,length(times2save)) );

figure
colormap jet
subplot(2,2,1)
contourf(times2save,frex,db_multitaper_tf,40,'linecolor','none')
set(gca,'clim',[-3 3])
colorbar
title('Power of Multitaper')
xlabel('time (ms)')
ylabel('frequcy (Hz)')

subplot(2,2,2)
contourf(times2save,frex,db_fft_tf,40,'linecolor','none')
set(gca,'clim',[-3 3])
colorbar
title('Power of FFT')
xlabel('time (ms)')
ylabel('frequcy (Hz)')

subplot(2,2,3)
contourf(times2save,frex,snr_multitaper,40,'linecolor','none')
set(gca,'clim',[0.4 1.6])
colorbar
title('SNR of Multitaper')
xlabel('time (ms)')
ylabel('frequcy (Hz)')

subplot(2,2,4)
contourf(times2save,frex,snr_fft,40,'linecolor','none')
set(gca,'clim',[0.4 1.6])
colorbar
title('SNR of FFT')
xlabel('time (ms)')
ylabel('frequcy (Hz)')

figure
colormap jet
contourf(times2save,frex,snr_multitaper-snr_fft,40,'linecolor','none')
colorbar
set(gca,'clim',[-0.2 1.2])
title('SNR_{Multipaper} minus SNR_{FFT}')
xlabel('time (ms)')
ylabel('frequcy (Hz)')

1. Select two frequencies, one relatively low and one relatively high (e.g., 8 Hz and 60 Hz), and compare the power time series and signal-to-noise time series in these frequency bands from the two methods in a separate figure, using line plots. Comment on the differences if there are any.

%% Select two frequencies, one relatively low and one relatively high (e.g., 8 Hz and 60 Hz), 
% and compare the power time series and signal-to-noise time series in these frequency bands 
% from the two methods in a separate figure, using line plots.
frequncies = [8, 60];
frequnciesidx = zeros(size(frequncies));
[~,frequnciesidx(1)] = min(abs(frex-frequncies(1)));
[~,frequnciesidx(2)] = min(abs(frex-frequncies(2)));

figure
for i = 1:length(frequncies)
    subplot(2,2,(i-1)*2+1)
    plot(times2save,db_multitaper_tf(frequnciesidx(i),:))
    title(['Frequency: ' num2str(frequncies(i)) ' Hz, Power'])
    xlabel('time (ms)')
    ylabel('Power (dB)')

    hold on
    plot(times2save,db_fft_tf(frequnciesidx(i),:))
    legend('Multitaper','FFT')

    subplot(2,2,(i-1)*2+2)
    plot(times2save,snr_multitaper(frequnciesidx(i),:))
    title(['Frequency: ' num2str(frequncies(i)) ' Hz, SNR'])
    xlabel('time (ms)')
    ylabel('SNR')
    
    hold on
    plot(times2save,snr_fft(frequnciesidx(i),:))
    legend('Multitaper','FFT')
end

• Multitaper常用于低信噪比的情况，如高频活动或功率的单试次估计。通过使用多个taper，multitaper方法在频率轴上引入了一定的平滑效应。这种平滑可以减少频率分辨率的精细度，使得频谱变得更加连续和平滑。

---

Exercises 18 | Baseline Normalization

1. Select three frequency bands and compute time-varying power at each electrode in these three bands, using either complex wavelet convolution or filter-Hilbert. Compute and store both the baseline-corrected power and the raw non-baseline-corrected power. You can choose which time period and baseline normalization method to use.

2. Select five time points and create topographical maps of power with and without baseline normalization at each selected time-frequency point. You should have time in columns and with/without baseline normalization in rows. Use separate figures for each frequency. The color scaling should be the same for all plots over time within a frequency, but the color scaling should be different for with versus without baseline normalization and should also be different for each frequency.

3. Are there qualitative differences in the topographical distributions of power with compared to without baseline normalization? Are the differences more prominent in some frequency bands or at some time points? What might be causing these differences?

%% Select three frequency bands and compute time-varying power at each electrode in these three bands, using either complex wavelet convolution or filter-Hilbert.
load sampleEEGdata.mat
frequency = [6 15 60];

% perform filter-Hilber
% parameters of filter
freqspread = [4 4 4]; % Hz +/- the center frequency
transwid   = 0.1;
nyquist = [1;1;1] .* EEG.srate/2;

% construct filter kernels
ffrequencies  = [ zeros(3,1) (1-transwid).*(frequency-freqspread)' (frequency-freqspread)' (frequency+freqspread)' (1+transwid).*(frequency+freqspread)' nyquist ]./nyquist(1);
idealresponse = [ 0 0 1 1 0 0 ];
raw_result = zeros(EEG.nbchan,length(frequency),EEG.pnts,EEG.trials);

for electrodei = 1:EEG.nbchan
    data2filter = squeeze(EEG.data(electrodei,:,:));
    for fi = 1:length(frequency)
        filterweights = firls(100,ffrequencies(fi,:),idealresponse);
        raw_result(electrodei,fi,:,:) = hilbert(filtfilt(filterweights,1,data2filter));
    end
end
raw_result = squeeze(mean(abs(raw_result).^2,4));
% define baseline period
baselinetime = [-500 -200 ]; % in ms

% convert baseline window time to indices
[~,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[~,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));

% dB-correct
baseline_power = mean(raw_result(:,:,baselineidx(1):baselineidx(2)),3);
baseline_dbconverted = 10*log10( bsxfun(@rdivide,raw_result,baseline_power));

timepoints = 0:100:400;
for i = 1:length(timepoints)
    [~,timepointsidx(i)] = min(abs(EEG.times-timepoints(i)));
end

%% plot 

for fi = 1:length(frequency)
    figure
    sgtitle(['frequncy: ' num2str(frequency(fi)) ' Hz'])
    max_raw = ceil(max(max(abs(raw_result(:,fi,timepointsidx)))));
    max_baseline = ceil(max(max(abs(baseline_dbconverted(:,fi,timepointsidx)))));
    for ti = 1:length(timepoints)
        subplot(2,length(timepoints)+1,ti)
        topoplot(raw_result(:,fi,timepointsidx(ti)),EEG.chanlocs,'electrodes','off','numcontour',0);
        title(['Raw Power, ' num2str(timepoints(ti)) 'ms'])   
        set(gca,'clim',[-max_raw,max_raw])

        subplot(2,length(timepoints)+1,length(timepoints)+1+ti);
        topoplot(baseline_dbconverted(:,fi,timepointsidx(ti)),EEG.chanlocs,'electrodes','off','numcontour',0);
        title(['DB Power, ' num2str(timepoints(ti)) 'ms'])
        set(gca,'clim',[-max_baseline,max_baseline])
    end
    subplot(subplot(2,length(timepoints)+1,length(timepoints)+1))
    axis off
    colorbar('Ticks',[0, 0.25, 0.5, 0.75, 1],'TickLabels',[-max_raw, -max_raw/2, 0, max_raw/2, max_raw]);
    
    subplot(subplot(2,length(timepoints)+1,2*(length(timepoints)+1)))
    axis off
    colorbar('Ticks',[0, 0.25, 0.5, 0.75, 1],'TickLabels',[-max_baseline, -max_baseline/2, 0, max_baseline/2, max_baseline]);
end

---

Exercises 19 | ITPC (Intertrial Phase Clustering)

1. Pick three electrodes. Compute time-frequency plots of ITPC and decibel-corrected power for these electrodes, using either complex Morlet wavelet convolution or the filter-Hilbert method. Plot the results side by side for each electrode (power and ITPC in subplots; one figure for each electrode). Are the patterns of results from ITPC and power generally similar or generally different? Do the results look more similar at some electrodes and less similar at other electrodes?

%% Pick three electrodes. Compute time-frequency plots of ITPC and decibel-corrected power for these electrodes, using either complex Morlet wavelet convolution or the filter-Hilbert method. 
load sampleEEGdata.mat
electrode = ["fz", "p8", "oz"];
electrodeidx = zeros(size(electrode));
for i = 1:length(electrode) 
    electrodeidx(i) = find(strcmpi({EEG.chanlocs.labels},electrode(i)));
end
frequencies = logspace(log10(3),log10(30),30);
%% calculate power

% perform filter-Hilber
% parameters of filter

freqspread = linspace(1,6,length(frequencies)); % Hz +/- the center frequency
transwid   = 0.1;
nyquist = repmat(EEG.srate/2, size(frequencies));

% construct filter kernels
ffrequencies  = [ ones(length(frequencies),1) (1-transwid).*(frequencies-freqspread)' (frequencies-freqspread)' (frequencies+freqspread)' (1+transwid).*(frequencies+freqspread)' nyquist' ]./nyquist(1);
idealresponse = [ 0 0 1 1 0 0 ];
filter_result = zeros(length(electrode),length(frequencies),EEG.pnts,EEG.trials);

for electrodei = 1:length(electrode)
    data2filter = squeeze(EEG.data(electrodeidx(electrodei),:,:));
    for fi = 1:length(frequencies)
        filterweights = firls(100,ffrequencies(fi,:),idealresponse);
        filter_result(electrodei,fi,:,:) = hilbert(filtfilt(filterweights,1,data2filter));
    end
end
power_result = squeeze(mean(abs(filter_result).^2,4));
% define baseline period
baselinetime = [-500 -200 ]; % in ms

% convert baseline window time to indices
[~,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[~,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));

% dB-correct
baseline_power = mean(power_result(:,:,baselineidx(1):baselineidx(2)),3);
power_dbconverted = 10*log10( bsxfun(@rdivide,power_result,baseline_power));



%% calculate ITPC
ITPC = abs(mean(exp(1i*angle(filter_result)),4));

%% plot
time2plot = [-300, 1000];
time2plotidx = zeros(size(time2plot));
for i = 1:2
    [~,time2plotidx(i)] = min(abs(EEG.times-time2plot(i)));
end

figure
for electrodei = 1:length(electrode)
    subplot(2,3,electrodei)
    contourf(EEG.times(time2plotidx(1):time2plotidx(2)),frequencies,squeeze(power_dbconverted(electrodei,:,time2plotidx(1):time2plotidx(2))),30,'linecolor','none');
    colormap jet
    set(gca,'clim',[-3 3])
    title(['Power from ' num2str(electrode(electrodei))])

    subplot(2,3,3+electrodei)
    contourf(EEG.times(time2plotidx(1):time2plotidx(2)),frequencies,squeeze(ITPC(electrodei,:,time2plotidx(1):time2plotidx(2))),30,'linecolor','none');
    colormap jet
    title(['ITPC from ' num2str(electrode(electrodei))])
end

1. For each of these three electrodes, compute wITPCz using reaction time as the trial-varying modulator. Perform this analysis for all time-frequency points to generate time-frequency maps of the relationship between phase and reaction time. Do the time-frequency maps of wITPCz look different from the time-frequency maps of ITPC? Do you see any striking patterns in the ITPCz results, and do the results differ across the different electrodes (don ’ t worry about statistics, base your judgment on qualitative patterns)? How would you interpret the results if they were statistically significant?

---

Exercises 20 | Total, Phase-Locked, and Non-Phase-Locked Power and Intertrial Phase Consistency

Pick two frequencies and compute total and non-phase-locked power from each electrode over time at these two frequencies. Pick two time windows, one early and one late, of several hundreds of milliseconds each (e.g., 100 – 300 ms and 500 – 800 ms) and show topographical maps of total power, non-phase-locked power, and phase-locked power from the average of all time points within these windows. Are there striking topographical differences among these results? If so, are the differences bigger or smaller in the early or the late time window? Why might this be the case?

%% Pick two frequencies and compute total and non-phase-locked power from each electrode over time at these two frequencies. 
% compute total power
load sampleEEGdata.mat
frequency = [5 20];

% perform filter-Hilber
% parameters of filter
freqspread = [1.5 4]; % Hz +/- the center frequency
transwid   = 0.15;
nyquist = [1;1] .* EEG.srate/2;

% construct filter kernels
ffrequencies  = [ zeros(2,1) (1-transwid).*(frequency-freqspread)' (frequency-freqspread)' (frequency+freqspread)' (1+transwid).*(frequency+freqspread)' nyquist ]./nyquist(1);
idealresponse = [ 0 0 1 1 0 0 ];
raw_result = zeros(EEG.nbchan,length(frequency),EEG.pnts,EEG.trials);

for electrodei = 1:EEG.nbchan
    data2filter = squeeze(EEG.data(electrodei,:,:));
    for fi = 1:length(frequency)
        filterweights = firls(200,ffrequencies(fi,:),idealresponse);
        raw_result(electrodei,fi,:,:) = hilbert(filtfilt(filterweights,1,data2filter));
    end
end
raw_result = squeeze(mean(abs(raw_result).^2,4));
% define baseline period
baselinetime = [-500 -200 ]; % in ms

% convert baseline window time to indices
[~,baselineidx(1)]=min(abs(EEG.times-baselinetime(1)));
[~,baselineidx(2)]=min(abs(EEG.times-baselinetime(2)));

% dB-correct
baseline_power = mean(raw_result(:,:,baselineidx(1):baselineidx(2)),3);
baseline_dbconverted = 10*log10( bsxfun(@rdivide,raw_result,baseline_power));
timewindow = [100,300; 600, 900];
timewindowidx = zeros(size(timewindow));

for i = 1:length(timewindowidx)
    [~,timewindowidx(i,1)] = min(abs(EEG.times-timewindow(i,1)));
    [~,timewindowidx(i,2)] = min(abs(EEG.times-timewindow(i,2)));
end

total = zeros(length(timewindow),EEG.nbchan,length(frequency));

for i = 1:length(timewindow)
    total(i,:,:) = squeeze(mean(baseline_dbconverted(:,:,timewindowidx(i,1):timewindowidx(i,2)),3));
end
%% compute non-phase-locked power
% compute ERP
erp = squeeze(mean(EEG.data(:,:,:),3));

% compute induced power by subtracting ERP from each trial
induced_EEG = EEG.data(:,:,:) - repmat(erp,[1 1 EEG.trials]);

% filter the induced eeg data
raw_result_non_phase_locked = zeros(EEG.nbchan,length(frequency),EEG.pnts,EEG.trials);

for electrodei = 1:EEG.nbchan
    data2filter = squeeze(induced_EEG(electrodei,:,:));
    for fi = 1:length(frequency)
        filterweights = firls(200,ffrequencies(fi,:),idealresponse);
        raw_result_non_phase_locked(electrodei,fi,:,:) = hilbert(filtfilt(filterweights,1,data2filter));
    end
end
raw_result_non_phase_locked = squeeze(mean(abs(raw_result_non_phase_locked).^2,4));
% baseline normalization
baseline_dbconverted_non_phase_locked = 10*log10( bsxfun(@rdivide,raw_result_non_phase_locked,baseline_power));

non_phase_locked = zeros(length(timewindow),EEG.nbchan,length(frequency));

for i = 1:length(timewindow)
    non_phase_locked(i,:,:) = squeeze(mean(baseline_dbconverted_non_phase_locked(:,:,timewindowidx(i,1):timewindowidx(i,2)),3));
end

%% compute phase-locked power
phase_locked = total - non_phase_locked;
%% plot
figure
n = 0;  % 用于绘子图时计算行数
for fi = 1:length(frequency)
    for twi = 1:length(timewindow)
        subplot(4,3,3*n+1)
        topoplot(squeeze(total(twi,:,fi)),EEG.chanlocs,'electrodes','off','numcontour',0);
        set(gca,'clim',[-2,2])
        title(['Total (' num2str(frequency(fi)) 'Hz, ' num2str(timewindow(twi,1)) '-' num2str(timewindow(twi,2)) 'ms)'])


        subplot(4,3,3*n+2)
        topoplot(squeeze(non_phase_locked(twi,:,fi)),EEG.chanlocs,'electrodes','off','numcontour',0);
        set(gca,'clim',[-2,2])
        title(['Non-Phase-Locked (' num2str(frequency(fi)) 'Hz, ' num2str(timewindow(twi,1)) '-' num2str(timewindow(twi,2)) 'ms)'])

        subplot(4,3,3*n+3)
        topoplot(squeeze(phase_locked(twi,:,fi)),EEG.chanlocs,'electrodes','off','numcontour',0);
        set(gca,'clim',[-1,1])
        title(['Phase-Locked (' num2str(frequency(fi)) 'Hz, ' num2str(timewindow(twi,1)) '-' num2str(timewindow(twi,2)) 'ms)'])

        n = n+1;
    end
    colormap jet
end

---

Tips

1. 把以ms为单位的时间转化为时间的索引index可以用dsearchn函数：

   % convert requested times to indices
   times2saveidx = dsearchn(EEG.times',times2save');