Make a 3D matrix*****SOS******

21

You should not crop to the smallest, because you don’t want to lose data, and resizing it could be confusing because you want the output shape to be time_points X n_electrodes. You should just leave it as a cell array as it can store inconsistent matrix sizes.

Try this code (replacing the variable direc with the path to your data):

direc = 'PATHTODATA';
files = dir(fullfile(direc,'*.mat'));
n_files = numel(files);
output_cell = cell(1,n_files);

for i = 1:n_files
loaded_file = load(fullfile(direc,files(i).name));
field_name = fieldnames(loaded_file);
data = loaded_file.(field_name{1});
output_cell{i} = data;
end
1 Like
22

Thanks in advance for your complete reponse!
I will try.

23

In this code, you did not make any 3D matrix
Am I right?

wanna pad all arrays to the size of the largest array (with what value?).
how can I do this strategy?

24

This could be dangerous down the line. It is probably best to keep the individual arrays as they are. Each subject’s data is indexed in the {nth} index of output_cell.

1 Like
25

yeep!
‘‘It is probably best to keep the individual arrays as they are.’’
in your code, without making 3D matrix, u did the same in a cell?
in this way no need to change the sizes, Am I right?
So how can I plot it?
A normal plot?

26

For more basic functions like plotting, you should refer to the MATLAB documentation, which is really well written. For example: 2-D line plot - MATLAB plot

You can also type help FUNCTIONNAME, e.g. help plot to see usage in the command window.

In general, you can extract a subject’s data with sub_data = output_cell{n} (where n is a number 1-12 corresponding to the subject), then extract a column with electrode_data = sub_data(:,k) where k is the number of the electrode you want to plot, then you can do plot(electrode_data). You can review the link I set above to do more customization of labels, if that’s something you would want to do.

The x-axis here just correspond with sample number, you would need the sampling rate to convert to seconds or milliseconds.

27

Dear Steven,
I cannot find a word to say really really thank you. You did so much to me.
Finally, it works :star_struck:

Bests,
Neda

1 Like
28

Hi Steven,
I’m filtering my data and extracting features. I’ve done for 7 similar files with same codes, but for one file I’ve got this error which is about ‘filtfilt’ function.
What should I do?

1
1889×363 8.91 KB

29

I would need to know more about what is being input into the function to help here.

30

%% Feature Extraction from EEG Rhythm in Time Domain, EO1, Active
clc
clear
close all
%% Load EO1, Active
load P01EO1
load P03EO1
load P05EO1
load P08EO1
load P09EO1
load P10EO1
load P11EO1
load P17EO1
load P18EO1
load P19EO1
load P21EO1
load P24EO1
load P25EO1
load P27EO1
%%
fs=256; % Sampling frequency(Hz)
Nf=36; %number of feautures
Nch= size(P01EO1,2); % number of columns(channels)
TfeaturesR_P01EO1= zeros(Nf,Nch);
TfeaturesR_P03EO1= zeros(Nf,Nch);
TfeaturesR_P05EO1= zeros(Nf,Nch);
TfeaturesR_P08EO1= zeros(Nf,Nch);
TfeaturesR_P09EO1= zeros(Nf,Nch);
TfeaturesR_P10EO1= zeros(Nf,Nch);
TfeaturesR_P11EO1= zeros(Nf,Nch);
TfeaturesR_P17EO1= zeros(Nf,Nch);
TfeaturesR_P18EO1= zeros(Nf,Nch);
TfeaturesR_P19EO1= zeros(Nf,Nch);
TfeaturesR_P21EO1= zeros(Nf,Nch);
TfeaturesR_P27EO1= zeros(Nf,Nch);
TfeaturesR_P25EO1= zeros(Nf,Nch);
TfeaturesR_P27EO1= zeros(Nf,Nch);
%% Step1: Denoising (50Hz)
% Design stop filter
fl= 49.9;
fh= 50.1;
order= 3;
wn= [fl fh]/ (fs/2);
type= ‘stop’;
[b1,a1]= butter(order,wn,type);
% Set up a for loop
for i=1:Nch
sig1= P01EO1(:,i);
sig2= P03EO1(:,i);
sig3= P05EO1(:,i);
sig4= P08EO1(:,i);
sig5= P09EO1(:,i);
sig6= P10EO1(:,i);
sig7= P11EO1(:,i);
sig8= P17EO1(:,i);
sig9= P18EO1(:,i);
sig10= P19EO1(:,i);
sig11= P21EO1(:,i);
sig12= P24EO1(:,i);
sig13= P25EO1(:,i);
sig14= P27EO1(:,i);
% Apply notch (stop) filter
sig1= filtfilt(b1,a1,sig1);
sig2= filtfilt(b1,a1,sig2);
sig3= filtfilt(b1,a1,sig3);
sig4= filtfilt(b1,a1,sig4);
sig5= filtfilt(b1,a1,sig5);
sig6= filtfilt(b1,a1,sig6);
sig7= filtfilt(b1,a1,sig7);
sig8= filtfilt(b1,a1,sig8);
sig9= filtfilt(b1,a1,sig9);
sig10= filtfilt(b1,a1,sig10);
sig11= filtfilt(b1,a1,sig11);
sig12= filtfilt(b1,a1,sig12);
sig13= filtfilt(b1,a1,sig13);
sig14= filtfilt(b1,a1,sig14);
%% Frequency Bands
% delta= 1-4 Hz,
% theta= 4-8 Hz,
% alpha= 8-12 Hz,
% beta= 12-20 Hz,
% low gamma= 20-25 Hz,
% gamma= 25-40 Hz
band= [1,4,8,12,20,25;
4,8,12,20,25,40];
%% Step2: Extract EEG Rhythms in Time Domain
for j=1:size(band,2); %first column:delta
% Design bandpass filter
wn= band(:,j)/ (fs/2);
[b2,a2]= butter(3,wn,‘bandpass’);
% Apply designed bandpass filter
sig_1= filtfilt(b2,a2,sig1);
sig_2= filtfilt(b2,a2,sig2);
sig_3= filtfilt(b2,a2,sig3);
sig_4= filtfilt(b2,a2,sig4);
sig_5= filtfilt(b2,a2,sig5);
sig_6= filtfilt(b2,a2,sig6);
sig_7= filtfilt(b2,a2,sig7);
sig_8= filtfilt(b2,a2,sig8);
sig_9= filtfilt(b2,a2,sig9);
sig_10= filtfilt(b2,a2,sig10);
sig_11= filtfilt(b2,a2,sig11);
sig_12= filtfilt(b2,a2,sig12);
sig_13= filtfilt(b2,a2,sig13);
sig_14= filtfilt(b2,a2,sig14);
tp_P01EO1(:,j)= myfeatureExtraction(sig_1);
tp_P03EO1(:,j)= myfeatureExtraction(sig_2);
tp_P05EO1(:,j)= myfeatureExtraction(sig_3);
tp_P08EO1(:,j)= myfeatureExtraction(sig_4);
tp_P09EO1(:,j)= myfeatureExtraction(sig_5);
tp_P10EO1(:,j)= myfeatureExtraction(sig_6);
tp_P11EO1(:,j)= myfeatureExtraction(sig_7);
tp_P17EO1(:,j)= myfeatureExtraction(sig_8);
tp_P18EO1(:,j)= myfeatureExtraction(sig_9);
tp_P19EO1(:,j)= myfeatureExtraction(sig_10);
tp_P21EO1(:,j)= myfeatureExtraction(sig_11);
tp_P24EO1(:,j)= myfeatureExtraction(sig_12);
tp_P25EO1(:,j)= myfeatureExtraction(sig_13);
tp_P27EO1(:,j)= myfeatureExtraction(sig_14);
end
%% Step3: Feature Extraction in time domain
TfeaturesR_P01EO1(:,i)= tp_P01EO1(:);
TfeaturesR_P03EO1(:,i)= tp_P03EO1(:);
TfeaturesR_P05EO1(:,i)= tp_P05EO1(:);
TfeaturesR_P08EO1(:,i)= tp_P08EO1(:);
TfeaturesR_P09EO1(:,i)= tp_P09EO1(:);
TfeaturesR_P10EO1(:,i)= tp_P10EO1(:);
TfeaturesR_P11EO1(:,i)= tp_P11EO1(:);
TfeaturesR_P17EO1(:,i)= tp_P17EO1(:);
TfeaturesR_P18EO1(:,i)= tp_P18EO1(:);
TfeaturesR_P19EO1(:,i)= tp_P19EO1(:);
TfeaturesR_P21EO1(:,i)= tp_P21EO1(:);
TfeaturesR_P24EO1(:,i)= tp_P24EO1(:);
TfeaturesR_P25EO1(:,i)= tp_P25EO1(:);
TfeaturesR_P27EO1(:,i)= tp_P27EO1(:);
% disp([‘iteration’,num2str(i)])
end
Tfeatures_Rhythms_EO1Active= [TfeaturesR_P01EO1;TfeaturesR_P03EO1;TfeaturesR_P05EO1;
TfeaturesR_P08EO1;TfeaturesR_P09EO1;TfeaturesR_P10EO1;
TfeaturesR_P11EO1;TfeaturesR_P17EO1;TfeaturesR_P18EO1;
TfeaturesR_P19EO1;TfeaturesR_P21EO1;TfeaturesR_P24EO1;
TfeaturesR_P25EO1;TfeaturesR_P27EO1];
save Tfeatures_allRhythms_EO1Active TfeaturesR_P01EO1 TfeaturesR_P03EO1 TfeaturesR_P05EO1 TfeaturesR_P08EO1 TfeaturesR_P09EO1 TfeaturesR_P10EO1 TfeaturesR_P11EO1 TfeaturesR_P17EO1 TfeaturesR_P18EO1 TfeaturesR_P19EO1 TfeaturesR_P21EO1 TfeaturesR_P24EO1 TfeaturesR_P25EO1 TfeaturesR_P27EO1
save Tfeatures_Rhythms_EO1Active Tfeatures_Rhythms_EO1Active

31

Can you just describe the variables that are being input to the filtfilt function throwing the error?

32

sig_10= filtfilt(b2,a2,sig10);

I’ve got the error in this line,
“sig10” is denoised signal which is an input

33

It looks like sig_10 may contain NaNs or Infs. You should investigate why it is happening, perhaps plotting sig10 before and after each filtering step can reveal where something went wrong. Here is how to replace those infinite values with a 0, which should at least stop the error from occuring. How to replace the array elements with "Nan" or "Inf" to either zero or some other number? -

34

I will look up it.
Thanks in advance for your answer.

Neda

35

Yesssssssssssssss,
Thatsss ittttttt,
It worksssssssssssssss :slight_smile: <3
Thhhhaaaaaaankssssssssssssssss many manyyyyy.
1st index of ‘sig10’ was NaN
I repaired by the one you suggested: sig10(isinf(sig10)|isnan(sig10)) = 0; and it works.

1 Like