Clear Filters
Clear Filters

I am struggling trying to figure out this code I have it completed but I keep getting errors, can you help?

4 views (last 30 days)
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

Answers (2)

Voss
Voss on 1 May 2024
Well, it doesn't take a crystal ball to see that this is not going to work (multiplying a 1x15 character vector by a 1x8 character vector)
tmefdata = 'mean(efdata, 3)' .* 'hann(eN)' ;
and it's likely not at all what was intended. My best guess at what was intended is
tmefdata = mean(efdata, 3) .* hann(eN).' ;
or, using some of the variables already calculated at that point in the code,
mefdata = mean(efdata, 3);
hanntaper = hann(eN).'; % Use the hann function and 'eN' as length. !!!!! Transpose as well. !!!!!
tmefdata = mefdata .* hanntaper ; % multiply the mefdata with hanntaper. Use element-wise multiplication '.*'
After that, you'll have some other errors when trying to plot some of the stuff. I made some more guesses.
% load('eegdata.mat');
eegdata = rand(3,1200);
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 = mefdata .* hanntaper ; % 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,:),srate); % Left side is already done, complete the right side
[tmef_freqs(i,:),tmef_pows(i,:)] = myFourier(tmefdata(i,:),srate); % 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(mef_freqs(i,:), mef_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(tmef_freqs(i,:), tmef_pows(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

John D'Errico
John D'Errico on 1 May 2024
Edited: John D'Errico on 1 May 2024
I need to laugh (sorry, but it is true) every time I see someone claim they have completely written a long piece of code, and claim it ti be correct, but they keep getting errors. That just means you were not complete in what you did. And of course in your code, I see some major errors in what you wrote.
More importantly, it tells me that you made a major mistake in HOW you are writing that code. What mistake? You need to write code SLOWLY, in pieces, segments, chunks. Make sure each piece works, before going on to the next step. Don't write a huge mess of code, and then just expect it to work.
So you might go this far:
load('eegdata.mat');
First, verify what was in that mat file! Use tools like whos. Make sure you completely understand what you are starting with. PLot it if necessary. In fact, plot EVERYTHING, even if you don't think it is necessary.
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
At this point, you might stop temporarily. Verify what is in those vectors. Insure they make sense. Again, PLOT THEM. As you get better at coding, you might bite off a little larger chunk at a time, but even I don't get too far ahead of things.
In this case, what do I see?
kernel
kernel =
Columns 1 through 7
0.056907 0.06081 0.06432 0.067343 0.069791 0.071594 0.072699
Columns 8 through 14
0.073071 0.072699 0.071594 0.069791 0.067343 0.06432 0.06081
Column 15
0.056907
Plot(kernel)
If you do as I did with your code, you would see this is NOT in fact a very complete Gaussian Kernel! It would appear you have given it far too little support. There is NOTHING out in the tails. If you look carefully at the peak versus min of that curve, you would see it is only a little more than 10%.
The point is, look carefully at what you have done. THINK about what you see.
...
When you are done, now you will know that EVERYTHING works, and does as expected. There will be no surprises.

Categories

Find more on Time-Frequency Analysis in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!