[Eeglablist] Apologies, MATLAB is not my first language... (spectopo() code error)
Makoto Miyakoshi
mmiyakoshi at ucsd.edu
Sun Nov 22 10:26:35 PST 2020
Dear Ian,
I ran your original code with my epoched data but could not replicate your
problem.
I modified your code for you. Please try it out if it works.
By the way, my code runs 190 times faster than yours.
Notice how I avoid repetitive use of spectopo().
Your code contains 9500 iterations of spectopo() applications inside the
triple loop of 50 x 5 x 38 for epochs, freqBans, and electrodes,
respectively.
But the latter two does not have to be performed separately. You just need
to 'extract' your freq bins of interest post hoc as I showed in my code.
Thus my code contains only 50 iterations, saving 9450 iterations. My code
finished the job in a minute.
Note also the element-wise division applied between the 3-D tensors. This
process in the original code is almost wrong, since 2-D matrix-wise
normalization is repeatedly applied while the value is stored progressively
one by one. You should avoid this.
> Apologies, MATLAB is not my first language.
My English sounds fucking shit even after living in the English speaking
countries over 12 years and constantly embarrasses me. I hope you feel
better now.
Makoto
%%%%%%%%%%%%%
% 11/22/2020 Makoto. For Ian.
% EEG =
pop_loadset('filename','G00001_Ct_MMN.set','filepath','/data/mobi/Daisuke/p0100_icSelected_epoched/');
%
% freq_ranges = {[1 4],[5 8],[9 12],[13 30],[31 64]};
% psd_data_norm = zeros([50 61 5]);
%
% saved_shape_data = EEG.data(:,:,1:100);
% saved_baseline_data = EEG.data(:,:,201:300);
%
% meanPowerMicroV = zeros(61,5);
% meanPowerMicroV_baseline = zeros(61,5);
%
% for dataEpochIdx = 1:50 % each trial
% current_data = squeeze(saved_shape_data(:,:,dataEpochIdx));
% current_baseline_data =
squeeze(saved_baseline_data(:,:,dataEpochIdx));
% meanPowerMicroV = zeros([61 5]);
% meanPowerMicroV_baseline = zeros([61 5]);
% for i = 1:5 % for all frequency bands
% for channelIdx = 1:38
%
% % For shape
% disp(channelIdx)
% [psdOutDb(channelIdx,:), freq] =
spectopo(current_data(channelIdx, :), 0, 128, 'plot', 'off');
% lowerFreqIdx = find(freq==freq_ranges{i}(1));
% higherFreqIdx = find(freq==freq_ranges{i}(2));
% meanPowerMicroV(channelIdx,i) =
mean(10.^((psdOutDb(channelIdx, lowerFreqIdx:higherFreqIdx))/10), 2);
%
% % For corresponding baseline (this is the aforementioned
baseline section)
% [psdOutDb_baseline(channelIdx,:), freq_baseline] =
spectopo(current_baseline_data(channelIdx, :), 0, 128, 'plot', 'off');
% lowerFreqIdx = find(freq_baseline==freq_ranges{i}(1));
% higherFreqIdx = find(freq_baseline==freq_ranges{i}(2));
% meanPowerMicroV_baseline(channelIdx,i) =
mean(10.^((psdOutDb_baseline(channelIdx,
lowerFreqIdx:higherFreqIdx))/10),2);
% end
% meanPowerMicroV_norm = meanPowerMicroV ./
meanPowerMicroV_baseline; % normalize
% psd_data_norm(:,:,dataEpochIdx) = squeeze(meanPowerMicroV_norm);
% end
% end
%% Modified by Makoto.
EEG =
pop_loadset('filename','G00001_Ct_MMN.set','filepath','/data/mobi/Daisuke/p0100_icSelected_epoched/');
freq_ranges = {[1 4],[5 8],[9 12],[13 30],[31 64]};
psd_data_norm = zeros([50 61 5]);
dataEpochRange = 1:100;
baselineEpochRange = 201:300;
saved_shape_data = EEG.data(:,:,dataEpochRange);
saved_baseline_data = EEG.data(:,:,baselineEpochRange);
meanPowerMicroV = zeros(EEG.nbchan, 5, length(dataEpochRange));
meanPowerMicroV_baseline = zeros(EEG.nbchan, 5, length(baselineEpochRange));
% Obtain data-epoch PSD.
for dataEpochIdx = 1:length(dataEpochRange) % For each data epoch.
current_data = squeeze(saved_shape_data(:,:,dataEpochIdx));
[allElecPsd, freqBins] = spectopo(current_data, 0, 128, 'plot', 'off');
if dataEpochIdx == 1
freqRange1Idx = find(freqBins >= freq_ranges{1}(1) & freqBins <=
freq_ranges{1}(2));
freqRange2Idx = find(freqBins >= freq_ranges{2}(1) & freqBins <=
freq_ranges{2}(2));
freqRange3Idx = find(freqBins >= freq_ranges{3}(1) & freqBins <=
freq_ranges{3}(2));
freqRange4Idx = find(freqBins >= freq_ranges{4}(1) & freqBins <=
freq_ranges{4}(2));
freqRange5Idx = find(freqBins >= freq_ranges{5}(1) & freqBins <=
freq_ranges{5}(2));
end
meanPowerMicroV(:,:,dataEpochIdx) =
[mean(10.^((allElecPsd(:,freqRange1Idx))/10), 2) ...
mean(10.^((allElecPsd(:,freqRange2Idx))/10),
2) ...
mean(10.^((allElecPsd(:,freqRange3Idx))/10),
2) ...
mean(10.^((allElecPsd(:,freqRange4Idx))/10),
2) ...
mean(10.^((allElecPsd(:,freqRange5Idx))/10),
2)];
end
% Obtain baseline-epoch PSD.
for baselineEpochIdx = 1:length(dataEpochRange) % For each baseline epoch.
current_baseline_data =
squeeze(saved_baseline_data(:,:,baselineEpochIdx));
[allElecPsd, freqBins] = spectopo(current_baseline_data, 0, 128,
'plot', 'off');
if baselineEpochIdx == 1
freqRange1Idx = find(freqBins >= freq_ranges{1}(1) & freqBins <=
freq_ranges{1}(2));
freqRange2Idx = find(freqBins >= freq_ranges{2}(1) & freqBins <=
freq_ranges{2}(2));
freqRange3Idx = find(freqBins >= freq_ranges{3}(1) & freqBins <=
freq_ranges{3}(2));
freqRange4Idx = find(freqBins >= freq_ranges{4}(1) & freqBins <=
freq_ranges{4}(2));
freqRange5Idx = find(freqBins >= freq_ranges{5}(1) & freqBins <=
freq_ranges{5}(2));
end
meanPowerMicroV_baseline(:,:,baselineEpochIdx) =
[mean(10.^((allElecPsd(:,freqRange1Idx))/10), 2) ...
mean(10.^((allElecPsd(:,freqRange2Idx))/10), 2) ...
mean(10.^((allElecPsd(:,freqRange3Idx))/10), 2) ...
mean(10.^((allElecPsd(:,freqRange4Idx))/10), 2) ...
mean(10.^((allElecPsd(:,freqRange5Idx))/10), 2)];
end
% Normalize data-epoch PSD with baseline-epoch PSD.
psd_data_norm = meanPowerMicroV./meanPowerMicroV_baseline;
On Sun, Nov 22, 2020 at 8:19 AM Ian Jackson <ianjackso at reed.edu> wrote:
> Hi EEGLAB folks,
>
> I've repurposed code from Makoto's useful code from the documentation (
>
> https://sccn.ucsd.edu/wiki/Makoto's_useful_EEGLAB_code#How_to_extract_EEG_power_of_frequency_bands_.2806.2F06.2F2020_updated.29
> )
>
> *The purpose* is to generate mean power for 5 frequency bands over 61
> components of ICA data. Each mean power generated for a trial/epoch is
> normalized by dividing by the mean power of the corresponding baseline.
>
> *The issue* is that, when running the following code, I get an error
> *"Unable
> to perform assignment because the indices on the left side are not
> compatible with the size of the right side." *The error occurs at the point
> in the loop where j=3, i = 1 and channelIdx = 1, though the code works
> completely fine without the baseline section.
>
> ************************************
> *The code*
>
> %% do PSD for each trial and normalize by baseline
> freq_ranges = {[1 4],[5 8],[9 12],[13 30],[31 64]};
> psd_data_norm = zeros([50 61 5]);
>
> for j = 1:50 % each trial
> current_data = squeeze(saved_shape_data(j,:,:));
> current_baseline_data = squeeze(saved_baseline_data(j,:,:));
> meanPowerMicroV = zeros([61 5]);
> meanPowerMicroV_baseline = zeros([61 5]);
> for i = 1:5 % for all frequency bands
> for channelIdx = 1:61
>
> % For shape
> disp(channelIdx)
> [psdOutDb(channelIdx,:), freq] =
> spectopo(current_data(channelIdx, :), 0, 128, 'plot', 'off');
> lowerFreqIdx = find(freq==freq_ranges{i}(1));
> higherFreqIdx = find(freq==freq_ranges{i}(2));
> meanPowerMicroV(channelIdx,i) = mean(10.^((psdOutDb(channelIdx,
> lowerFreqIdx:higherFreqIdx))/10), 2);
>
> % For corresponding baseline (this is the aforementioned
> baseline section)
> [psdOutDb_baseline(channelIdx,:), freq_baseline] =
> spectopo(current_baseline_data(channelIdx, :), 0, 128, 'plot', 'off');
> lowerFreqIdx = find(freq_baseline==freq_ranges{i}(1));
> higherFreqIdx = find(freq_baseline==freq_ranges{i}(2));
> meanPowerMicroV_baseline(channelIdx,i) =
> mean(10.^((psdOutDb_baseline(channelIdx, lowerFreqIdx:higherFreqIdx))/10),
> 2);
> end
> meanPowerMicroV_norm = meanPowerMicroV ./ meanPowerMicroV_baseline; %
> normalize
> psd_data_norm(j,:,:) = squeeze(meanPowerMicroV_norm);
> end
> end
>
> ************************************
>
> Here is some more useful information about the variables at the point of
> error in the loop:
> size(saved_shape_data) = [50 61 576]
> size(saved_baseline_data) = [50 61 192]
> size(spectopo(current_data(channelIdx, :), 0, 128, 'plot', 'off')) = [1 65]
> size(freq) = [65 1] (it is an array [0,1,...,64])
> size(spectopo(current_baseline_data(channelIdx, :), 0, 128, 'plot', 'off'))
> = [1 65]
> size(freq_baseline) = [128 1] (it is an array [0,1,..,127])
>
> ************************************
>
> I suspect it has to do with freq_baseline due to its strangely large size,
> but I am at a loss as to how to fix the error. Any help would be greatly
> appreciated!
>
> The general upkeep with correspondence on this forum is great to see. Thank
> you and I look forward to hearing from you!
>
> Regards,
> Ian Jackson
> Reed College
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu
>
More information about the eeglablist
mailing list