load('eegdata.mat');
Error using load
Unable to find file or directory 'eegdata.mat'.
data = eegdata;
srate = 500;
N = size(data, 2); % Length of data on a channel; total number of data points
nchans = size(data, 1) ; % Number of channels, which is the size of the first dim of 'data'
time = (0:N-1) / srate; % Time vector
ntrials = 10; % Number of trials
% Lowpass Filter --------------
% You will create the gaussian kernel here for low-pass filtering. Use
% a steepness value of 0.5 and kernel length of 15.
steepness = 0.5; % kernel steepness
kernelLength = 15; % kernel length
incr = 2*steepness/(kernelLength-1);
kernel = exp(-(-steepness:incr:steepness).^2);
kernel = kernel./sum(kernel); % normalize kernel
% Filtering
% Use a loop to go through the data on each channel and
% use the 'conv' function to convolve the data with the gaussian
% kernel you created above. Ensure that the output has the same
% length as the input data, using the 'same' option. The fdata matrix
% will store the filtered data.
for i = 1:nchans
fdata(i,:) = conv(data(i,:), kernel, 'same');
end
% FIGURE 1: Plot Raw and Filtered Channel Data
% Plot the raw and filtered data ('data' and 'fdata') for each channel
% separately, in the first and second columns. Modify the x/y labels and
% the plot titles as shown in the figure in the question prompt.
figure('Name', 'Channel Data')
for i = 1:nchans
subplot(nchans, 2, (i-1)*2 + 1)
plot(time, data(i,:))
xlabel('Time (s)')
ylabel('Amplitude')
title(['Raw Data - Channel ', num2str(i)])
subplot(nchans, 2, (i-1)*2 + 2)
plot(time, fdata(i,:))
xlabel('Time (s)')
ylabel('Amplitude')
title(['Filtered Data - Channel ', num2str(i)])
end
% Epoch data --------------
% Epoch the data into epochs using the 'reshape' function.
% The number of epochs (trials) are indicated by the 'ntrials' variable.
% Data in each channel should be divided into 'ntrials' epochs and
% the resulting data should be stored in the 'efdata' matrix.
% The 'efdata' is a 3-dimensional matrix; with channels, data points, and trials
% as the dimensions, respectively.
efdata= reshape(fdata,nchans,[],ntrials);
% Average data across the epochs
% Calculate the mean data for each channel.
% You should take the average of all trials, separately for each channel.
% The 'mean' function allows you to indicate over which dimension
% should the mean be calculated. You should take the mean over trials.
mefdata = mean(efdata, 3);
eN = size(efdata, 2); % Calculate the length of a single epoch (number of data points in an epoch)
etime = 0:1/srate:(eN/srate)-1/srate; % Create the time vector for a single epoch
% Taper averaged epochs --------------
% As we discussed in class, edges in epoched data cause artifacts in
% frequency analysis. To avoid this problem, we will taper the data with
% a Hann window. Use the 'hann' function to create a hann window
% at the lenght of an epoch (you will need to transpose the output),
% then multiply the hann window with the averaged
% epoched data ('mefdata'). Remember, this is an element-wise
% multiplication, meaning that each point in the hann window will be
% multiplied with the matching point in the data, separately for each
% channel (use '.*'). Store the tapered data matrix in 'tmefdata'/
hanntaper = hann(eN); % Use the hann function and 'eN' as length. Transpose as well.
tmefdata = 'mean(efdata, 3)' .* 'hann(eN)' ; % multiply the mefdata with hanntaper. Use element-wise multiplication '.*'
% Frequency Analysis --------------
% Before you conduct the frequency analysis, go to the bottom of the script
% and follow the instructions to complete the 'myFourier' function. Once
% done, come back here to continue.
% Go through each channel and calculate the frequencies and powers for both
% the untapered (mefdata) and tapered (tmefdata) data.
for i=1:nchans
[mef_freqs(i,:),mef_pows(i,:)] = myFourier(mefdata(i,:)); % Left side is already done, complete the right side
[tmef_freqs(i,:),tmef_pows(i,:)] = myFourier(tmefdata(i,:)); % Left side is already done, complete the right side
end
% Time-Frequency Analysis--------------
% Before you conduct the time-frequency analysis, go to the bottom of the
% script and follow the instructions to complete the 'myMorlet' function.
% Once done, come back here to continue.
% Let's say we are interested in frequencies between 0 and 40 Hz
minFreq = 0;
maxFreq = 40;
% Calculate the frequency vectors and matching power values over time using
% the 'myMorlet' function. Remember, we don't use the epoched data for frequency analysis
% First use the continous filtered data, 'fdata' for the transformation to
% the time-frequency domain, then epoch the t-f data.
for i=1:nchans
[tf_freqs{i}, tf_pows{i}] = myMorlet(fdata(i,:), srate, minFreq, maxFreq); % Only complete the right side. Use 'fdata', 'srate', 'minFreq' and 'maxFreq' inside myMorlet
end
nfreqs = size(tf_freqs{1},1); % The number of frequencies
etf_pows = zeros(nchans, nfreqs, eN, ntrials); % Initiate matrix for epoched t-f data
metf_pows = zeros(nchans, nfreqs, eN); % Initiate matrix for averaged epoched t-f data
% Epoch the t-f data and calculate the mean t-f data
% This part is already done, though inspect it to see how it was done
for i = 1:nchans % calculate etf_pows and metf_pows
etf_pows(i,:,:,:) = reshape(tf_pows{i}, nfreqs,[], ntrials);
metf_pows(i,:,:) = squeeze(mean(etf_pows(i,:,:,:),4));
end
% Plots
nC = 4; % Number of types of plots
% FIGURE 2
figure('Name','Analysis Results', 'Position', [10 200 900 nC*200]) % Create figure 2, adjust its size
% Here you will create Figure 2, where averaged TD (time-domain), averaged untapered FD (frequency-domain),
% averaged Tapered FD, and the averaged TF (time-frequency domain) data are plotted, separately for each channel
% Use subplots in the loop to create the figure. The figure created should
% look like the one in the question (Figure 2).
for i = 1:nchans
subplot(nchans, nC, (i-1)*nC + 1) % Use 'nchans' and 'nC' variables to adjust the indices
plot(etime,mefdata(i,:)) % Plot epoch time vs mean epoched, filtered data.
title('Averaged Signal TD')
xlim([min(etime), max(etime)]) % Use min and max etime for the range of x values
ylabel(['Channel ', num2str(i)]) % The ylabel should show the channel number (see Fig. 2 in the question)
subplot(nchans, nC, (i-1)*nC + 2) % Use 'nchans' and 'nC' variables to adjust the indices
plot(tf_freqs{i}, tf_pows{i}) % Plots freqs vs pows for untapered data
title('Avg. Signal FD')
xlim([minFreq, maxFreq]) % Use minFreq and maxFreq as limits
subplot(nchans, nC, (i-1)*nC + 3) % Use 'nchans' and 'nC' variables to adjust the indices
plot(tf_freqs{i}, tmefdata(i, :)) % Plots freqs vs pows for tapered data
title('Avg. Tapered Signal FD')
xlim([minFreq, maxFreq]) % Use minFreq and maxFreq as limits
subplot(nchans,nC,(i-1)*nC + 4);
% Use etime, tf_freqs, and metf_pows for the contourf plot. Be careful
% about whether you use () or {} and the indices. Consider using the
% squeeze function as well. Number of contour levels is already set as
% 40, and 'linecolor' is set to 'none'.
contourf(etime, squeeze(tf_freqs{i}), squeeze(metf_pows(i, :, :)),40,'linecolor','none')
title('Signal TF')
end
% Coherence Analysis
% Here you will study the coherence in the signal between pairs of
% channels.
% Initialize the cell array to store coherence results
cohRes = cell(nchans, nchans);
% FIGURE 3
figure('Name','Coherence Analysis', 'Position', [10 200 800 800])
% Loop over each pair of channels
for i = 1:nchans
for j = (i+1):nchans
% Extract data for the current channel pair
chanA_data = data(i, :);
chanB_data = data(j, :);
% Perform spectral coherence analysis
% Check the 'mscohere' documentation. Use chanA_data and chanB_data
% for the two signals. Use the default window and nooverlap values
% (just put '[]' to use the default values). Check frequencies 0 to
% 40. Use srate for the sampling rate.
[cohVals, cohFreqs] = mscohere(chanA_data, chanB_data, [], [], 0:40, srate);
% Store coherence results in the cell array
cohRes{i, j} = cohVals;
% Plot coherence results in subplots
subplot(nchans, nchans, (i-1)*nchans + j); % Cool indexing magic here
plot(cohFreqs, cohVals); % Plot coherence frequencies vs coherence values
xlabel('Frequency (Hz)');
ylabel('Coherence');
title(['Chan. ', num2str(i), ' & ', num2str(j)]);
end
end
% My Fourier Transform Function
% The 'myFourier' function takes the signal and sampling rate as inputs
% and produces a vector frequencies and a vector of matching power values
% Use the 'fft' function to calculate the Fouerier coefficients.
function [freqs, pows] = myFourier(signal,srate)
n = length(signal); % length of signal
freqs = linspace(0, srate/2, n/2+1); % frequency vector
signalx = fft(signal); % conduct fft with the signal
pows = abs(signalx(1:n/2+1)).^2; % calculate power values based on the coefficients
end
% Morlet Convolution
% The 'myMorlet' function takes the signal, sampling rate, and min & max
% frequencies as inputs, and outputs the frequencies and matching powers
% over the time series. Use Matlab's 'cwt' function with the 'amor' option
% to conduct a Morlet wavelet transofrmation. Use the 'FrequencyLimits'
% option with 'minFreq' and 'maxFreq' as limits. Check the 'cwt'
% documentation for info on how to use these options.
function [freqs, pows] = myMorlet(signal,srate,minFreq,maxFreq)
[coefs, freqs] = cwt(signal, 'amor', srate, 'FrequencyLimits', [minFreq, maxFreq]); % Use the 'cwt' function with the 'amor' option to calculate the coefficients and frequencies
pows =abs(coefs).^2; % Use the 'abs' function to calculate ampltitudes from 'coefs' and take its square to calculate power
end
