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

Question

Olivia 2024-4-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2113696-i-am-struggling-trying-to-figure-out-this-code-i-have-it-completed-but-i-keep-getting-errors-can-yo

编辑： John D'Errico 2024-5-1

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

2 个评论
显示无隐藏无

Torsten 2024-4-30

We can neither run your code nor know the error message you are getting.

Image Analyst 2024-5-1

编辑：Image Analyst 2024-5-1

If you have any more questions, then attach your mat file with the paperclip icon after you read this:

TUTORIAL: How to ask a question (on Answers) and get a fast answer

Or, wait for Walter Roberson to see your post.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Voss 2024-5-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2113696-i-am-struggling-trying-to-figure-out-this-code-i-have-it-completed-but-i-keep-getting-errors-can-yo#answer_1450906

在 MATLAB Online 中打开

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

% 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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

John D'Errico 2024-5-1

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2113696-i-am-struggling-trying-to-figure-out-this-code-i-have-it-completed-but-i-keep-getting-errors-can-yo#answer_1450936

编辑：John D'Errico 2024-5-1

在 MATLAB Online 中打开

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.