From 5cb6252a0a879ce311d1ef246420aa8859d056b9 Mon Sep 17 00:00:00 2001 From: dnaligase <60392997+dnaligase@users.noreply.github.com> Date: Tue, 12 Apr 2022 13:20:09 +0200 Subject: [PATCH 01/21] BaseRaw upd --- .gitignore | 1 + BaseRaw.m | 33 +++++++++++++++++++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/.gitignore b/.gitignore index b853f42..c54bfb5 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ ch_names.mat +*.asv data diff --git a/BaseRaw.m b/BaseRaw.m index 382dabd..942b6cb 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -46,6 +46,39 @@ end end + + function [r, meta] = windowedPower(obj, time_window, noverlap, ... + lfreq, hfreq, verbose) + % calculate powers in a windowed fashion + arguments + obj + time_window (1,1) double + noverlap (1,1) double + lfreq (1,1) = 8 + hfreq (1,1) = 13 + verbose logical = true + end + + rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... + (obj.fs*time_window - obj.fs*noverlap)) + 1; + matrix = zeros(rowsNo, size(obj.data,1)); + + for j = 1:rowsNo + begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + stop = time_window*obj.fs + ... + (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + + powers = obj.sum_power_segment((begin:stop), lfreq, hfreq); + matrix(j, :) = powers; + end + + r = matrix; + meta = struct('time_window', time_window, 'noverlap', noverlap); + + if verbose + disp(meta); + end + end function r = sum_power(obj, l_freq, h_freq) r = sum_freq_band(obj.psd, obj.freq, l_freq, h_freq); From a92b14e8d77f4bc20d93fc4b0bc159150588b582 Mon Sep 17 00:00:00 2001 From: dnaligase <60392997+dnaligase@users.noreply.github.com> Date: Wed, 13 Apr 2022 12:26:38 +0200 Subject: [PATCH 02/21] upd scripts --- BaseRaw.m | 19 ++++++++++++++----- clean_and_viz_data.m | 5 +---- 2 files changed, 15 insertions(+), 9 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 942b6cb..f70b698 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -14,6 +14,7 @@ end properties (Hidden = true) time_as_index + id = 1 end methods @@ -26,9 +27,16 @@ if nargin >= 1 obj.data = EEG.data(picks, :); obj.fs = EEG.srate; - obj.subject = EEG.subject; - if isempty(EEG.times) + if ~isfield(EEG, 'subject') || isempty(EEG.subject) + obj.subject = string(obj.id); + increment_id = @(x) x + 1; + obj.id = increment_id(obj.id); + else + obj.subject = EEG.subject; + end + + if ~isfield(EEG, 'times') || isempty(EEG.times) obj.times = obj.create_times(); else % convert EEGLAB times to ms @@ -36,7 +44,7 @@ end [obj.psd, obj.freq] = pwelch(obj.data',[],[],256,obj.fs); - if ~isempty(EEG.chanlocs) + if isfield(EEG, 'chanlocs') obj.chanlocs = EEG.chanlocs; else obj.chanlocs = []; @@ -95,9 +103,9 @@ function crop(obj, tmin, tmax) arguments obj tmin {mustBeGreaterThanOrEqual(tmin, 0)} - tmax {mustBeReal} + tmax {mustBeReal} = obj.times(end) end - + start = tmin * obj.fs + 1; stop = tmax * obj.fs + 1; @@ -131,6 +139,7 @@ function crop(obj, tmin, tmax) obj return_figure logical = false end + assert(~isempty(obj.chanlocs), "Channel locations not provided.") figure; plot3([obj.chanlocs.X],[obj.chanlocs.Y],[obj.chanlocs.Z], 'o'); title(('3D channel location for #' + string(obj.subject))) diff --git a/clean_and_viz_data.m b/clean_and_viz_data.m index 1e65033..48bb686 100644 --- a/clean_and_viz_data.m +++ b/clean_and_viz_data.m @@ -1,7 +1,4 @@ % initialize EEGLAB variables -if ~exist("ALLEEG", "var") - eeglab -end load ch_names.mat fnames = dir("**/*_eyesclosed_*.hdf5"); @@ -27,7 +24,7 @@ EEG = pop_importdata('dataformat','array','nbchan',64,'subject',patientName(end),'data','data', ... 'srate',250,'pnts',0,'xmin',0, ... - 'chanlocs','/Volumes/EXTSTORAGE/Labrotation/64_channel.ced'); + 'chanlocs','64_channel.ced'); EEG = eeg_checkset( EEG ); EEG = pop_clean_rawdata(EEG, 'FlatlineCriterion','off', ... 'ChannelCriterion','off','LineNoiseCriterion','off', ... From aef238cb0a63d3e27155615ce4bea05828091b8f Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Wed, 13 Apr 2022 15:13:50 +0200 Subject: [PATCH 03/21] Update BaseRaw.m --- BaseRaw.m | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 942b6cb..a037ea1 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -67,7 +67,7 @@ begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; stop = time_window*obj.fs + ... (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - + powers = obj.sum_power_segment((begin:stop), lfreq, hfreq); matrix(j, :) = powers; end @@ -86,10 +86,17 @@ function r = sum_power_segment(obj, segment, l_freq, h_freq) transposed = obj.data'; % segment in datapoints as ``double`` - [psd_window, freq_window] = pwelch(transposed(segment, :), ... + art_mat = abs(transposed(segment,:)) > 100; + art_vec = ~any(art_mat); + psd_window = NaN(129,64); + freq_window = NaN(129,1); + r = NaN(1,64); + if sum(art_vec) > 0 + [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... [],[],256,obj.fs); r = sum_freq_band(psd_window, freq_window, l_freq, h_freq); + end end function crop(obj, tmin, tmax) arguments From 0eac3f39ea338b00c9720ce79aff71c8e9916042 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Thu, 14 Apr 2022 13:08:55 +0200 Subject: [PATCH 04/21] Added new medthod filterEEG --- BaseRaw.m | 55 ++++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 54 insertions(+), 1 deletion(-) diff --git a/BaseRaw.m b/BaseRaw.m index a815b11..e460d43 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -5,6 +5,7 @@ times psd freq + raw end properties (SetAccess = private) fs @@ -55,6 +56,51 @@ end + function [r, meta] = windowedPower(obj, time_window, noverlap, ... + lfreq, hfreq, verbose) + % calculate powers in a windowed fashion + arguments + obj + time_window (1,1) double + noverlap (1,1) double + lfreq (1,1) = 8 + hfreq (1,1) = 13 + verbose logical = true + end + + rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... + (obj.fs*time_window - obj.fs*noverlap)) + 1; + matrix = zeros(rowsNo, size(obj.data,1)); + + for j = 1:rowsNo + begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + stop = time_window*obj.fs + ... + (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + + powers = obj.sum_power_segment((begin:stop), lfreq, hfreq); + matrix(j, :) = powers; + end + + r = matrix; + meta = struct('time_window', time_window, 'noverlap', noverlap); + + if verbose + disp(meta); + end + end + + function filterEEG(obj,hi,lo) + obj.raw = obj.data; + if ~isempty(hi) + obj.data = hifi(obj.data', 1e6/obj.fs, hi)'; + end + if ~isempty(lo) + obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; + end + end + + + function [r, meta] = windowedPower(obj, time_window, noverlap, ... lfreq, hfreq, verbose) % calculate powers in a windowed fashion @@ -149,6 +195,10 @@ function crop(obj, tmin, tmax) assert(~isempty(obj.chanlocs), "Channel locations not provided.") figure; plot3([obj.chanlocs.X],[obj.chanlocs.Y],[obj.chanlocs.Z], 'o'); + for i = 1:64 + locs(i) = {obj.chanlocs(i).labels}; + end + text([obj.chanlocs.X],[obj.chanlocs.Y],[obj.chanlocs.Z],locs,'VerticalAlignment','bottom','HorizontalAlignment','right','FontSize',6) title(('3D channel location for #' + string(obj.subject))) xlabel('X') ylabel('Y') @@ -162,7 +212,10 @@ function crop(obj, tmin, tmax) methods (Access = protected) function r = create_times(obj) - r = reshape(double( 0 : length(obj.data) ) * 1/obj.fs, 1, []); + time_step = 1/obj.fs; + endpoint = length(obj.data)/obj.fs; + r = [time_step:time_step:endpoint]; +% r = reshape(double( 1/obj.fs : length(obj.data) ) * 1/obj.fs, 1, []); end function r = apply_for_signal(obj, func, channel_wise) From 6fd9d410dfd636629532527d60b248a00358e75c Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Thu, 14 Apr 2022 13:19:15 +0200 Subject: [PATCH 05/21] Remove windowedPower as it was accidentaly copied twice --- BaseRaw.m | 34 +--------------------------------- 1 file changed, 1 insertion(+), 33 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index e460d43..74f2588 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -56,39 +56,7 @@ end - function [r, meta] = windowedPower(obj, time_window, noverlap, ... - lfreq, hfreq, verbose) - % calculate powers in a windowed fashion - arguments - obj - time_window (1,1) double - noverlap (1,1) double - lfreq (1,1) = 8 - hfreq (1,1) = 13 - verbose logical = true - end - - rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... - (obj.fs*time_window - obj.fs*noverlap)) + 1; - matrix = zeros(rowsNo, size(obj.data,1)); - - for j = 1:rowsNo - begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - stop = time_window*obj.fs + ... - (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - - powers = obj.sum_power_segment((begin:stop), lfreq, hfreq); - matrix(j, :) = powers; - end - - r = matrix; - meta = struct('time_window', time_window, 'noverlap', noverlap); - - if verbose - disp(meta); - end - end - + function filterEEG(obj,hi,lo) obj.raw = obj.data; if ~isempty(hi) From 07ea25003c43e19925142b8672eb1030bbbe021c Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Thu, 14 Apr 2022 16:54:52 +0200 Subject: [PATCH 06/21] Added function --- BaseRaw.m | 2 -- plot_mult_topo.m | 60 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 2 deletions(-) create mode 100644 plot_mult_topo.m diff --git a/BaseRaw.m b/BaseRaw.m index 74f2588..01a9238 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -66,8 +66,6 @@ function filterEEG(obj,hi,lo) obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; end end - - function [r, meta] = windowedPower(obj, time_window, noverlap, ... lfreq, hfreq, verbose) diff --git a/plot_mult_topo.m b/plot_mult_topo.m new file mode 100644 index 0000000..7b62294 --- /dev/null +++ b/plot_mult_topo.m @@ -0,0 +1,60 @@ +function plot_mult_topo(mat,chans,desc,lims,time_window,noverlap,fname,fps) +% mat: Rows: Timeseries +% Cols: Channels +% If more than two dimensions are provided, topoplots +% will be plottet into subplots +% chans: Channel Locations +% desc: Description/ Title of the plot +% lims: Limits of the colorbar, either als two column +% vector or matrix if different limits should be +% applied for the different plots +%time_window and noverlap: settings of the calculations +% fname filename of the produced video 'Test.avi' +% fps: Frames per second, e.g. 10 + + +%% Define number of subplots +num_plots = size(mat,3); +rows = ceil(num_plots/2); +cols = ceil(num_plots/rows); +%% Check input arguments +if size(lims,1) == 1 + lims = repmat(lims,[num_plots,1]); +end +%% Define number of frames +framesNo = size(mat,1); +% allocate memory for future figures +F = struct('cdata', cell(1,framesNo), 'colormap', cell(1,framesNo)); +pb = CmdLineProgressBar(("Preparing "+framesNo+" frames using "+noverlap+"s overlap and "+time_window+"s window size... ")); + +%% Create Figure and axes +figure('WindowState','maximized','Visible','off') +for i = 1:framesNo + pb.print(i,framesNo) + for j = 1:num_plots + ax(j) = subplot(rows,cols,j); + sgtitle([num2str(time_window + (time_window - noverlap) * (i-1)), ' s']) + title(desc{1,j}); + topoplot(mat(i,:,j),chans); + colorbar + caxis([lims(j,1),lims(j,2)]) + F(i) = getframe(gcf); + end +end + +writerObj = VideoWriter(fname); +writerObj.FrameRate = fps; +disp("Creating " + length(F) / writerObj.FrameRate + " s video..."); +% open the video writer +open(writerObj); +% write the frames to the video +for i=1:length(F) + % convert the image to a frame + frame = F(i); + writeVideo(writerObj, frame); +end +% close the writer object +close(writerObj); +disp('Done.') +disp(strcat('Find it as ', pwd, filesep, fname)) +end \ No newline at end of file From a5a643a5c11813c9c91ccb3d2179d443b8641990 Mon Sep 17 00:00:00 2001 From: dnaligase <60392997+dnaligase@users.noreply.github.com> Date: Thu, 14 Apr 2022 17:11:43 +0200 Subject: [PATCH 07/21] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c54bfb5..431d2fa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ ch_names.mat *.asv data +.DS_Store From d83ead115b8b3dfeb73b1ff9cfede343cd7bff02 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Fri, 15 Apr 2022 12:34:11 +0200 Subject: [PATCH 08/21] Add field to Base Raw Add number of channel property, adapt chanlocs based on user picks. Change time calculation in the topoplot --- BaseRaw.m | 10 ++++++---- plot_mult_topo.m | 6 +++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 01a9238..e0b8fe2 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -6,6 +6,7 @@ psd freq raw + no_chan end properties (SetAccess = private) fs @@ -28,6 +29,7 @@ if nargin >= 1 obj.data = EEG.data(picks, :); obj.fs = EEG.srate; + obj.no_chan = size(obj.data,1); if ~isfield(EEG, 'subject') || isempty(EEG.subject) obj.subject = string(obj.id); @@ -46,7 +48,7 @@ [obj.psd, obj.freq] = pwelch(obj.data',[],[],256,obj.fs); if isfield(EEG, 'chanlocs') - obj.chanlocs = EEG.chanlocs; + obj.chanlocs = EEG.chanlocs(1,picks); else obj.chanlocs = []; end @@ -108,9 +110,9 @@ function filterEEG(obj,hi,lo) % segment in datapoints as ``double`` art_mat = abs(transposed(segment,:)) > 100; art_vec = ~any(art_mat); - psd_window = NaN(129,64); + psd_window = NaN(129,obj.no_chan); freq_window = NaN(129,1); - r = NaN(1,64); + r = NaN(1,obj.no_chan); if sum(art_vec) > 0 [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... [],[],256,obj.fs); @@ -161,7 +163,7 @@ function crop(obj, tmin, tmax) assert(~isempty(obj.chanlocs), "Channel locations not provided.") figure; plot3([obj.chanlocs.X],[obj.chanlocs.Y],[obj.chanlocs.Z], 'o'); - for i = 1:64 + for i = 1:obj.no_chan locs(i) = {obj.chanlocs(i).labels}; end text([obj.chanlocs.X],[obj.chanlocs.Y],[obj.chanlocs.Z],locs,'VerticalAlignment','bottom','HorizontalAlignment','right','FontSize',6) diff --git a/plot_mult_topo.m b/plot_mult_topo.m index 7b62294..6360b86 100644 --- a/plot_mult_topo.m +++ b/plot_mult_topo.m @@ -1,4 +1,4 @@ -function plot_mult_topo(mat,chans,desc,lims,time_window,noverlap,fname,fps) +function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) % mat: Rows: Timeseries % Cols: Channels % If more than two dimensions are provided, topoplots @@ -25,7 +25,7 @@ function plot_mult_topo(mat,chans,desc,lims,time_window,noverlap,fname,fps) framesNo = size(mat,1); % allocate memory for future figures F = struct('cdata', cell(1,framesNo), 'colormap', cell(1,framesNo)); -pb = CmdLineProgressBar(("Preparing "+framesNo+" frames using "+noverlap+"s overlap and "+time_window+"s window size... ")); +pb = CmdLineProgressBar(("Preparing "+framesNo+ " Frames... ")); %% Create Figure and axes figure('WindowState','maximized','Visible','off') @@ -33,7 +33,7 @@ function plot_mult_topo(mat,chans,desc,lims,time_window,noverlap,fname,fps) pb.print(i,framesNo) for j = 1:num_plots ax(j) = subplot(rows,cols,j); - sgtitle([num2str(time_window + (time_window - noverlap) * (i-1)), ' s']) + sgtitle([num2str(time_vec(i)), ' s']) title(desc{1,j}); topoplot(mat(i,:,j),chans); colorbar From 3f21d08a9fba1e968fc9b52cf573eee4fe4bd0ac Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Fri, 15 Apr 2022 14:14:41 +0200 Subject: [PATCH 09/21] add pvpmod --- BaseRaw.m | 89 ++++++++++++++++++++++++++++++------------------------- pvpmod.m | 40 +++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 41 deletions(-) create mode 100644 pvpmod.m diff --git a/BaseRaw.m b/BaseRaw.m index e0b8fe2..0c209c9 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -28,6 +28,7 @@ if nargin >= 1 obj.data = EEG.data(picks, :); + obj.raw = EEG.data(picks, :); obj.fs = EEG.srate; obj.no_chan = size(obj.data,1); @@ -57,19 +58,41 @@ end end - - + function filterEEG(obj,hi,lo) - obj.raw = obj.data; - if ~isempty(hi) - obj.data = hifi(obj.data', 1e6/obj.fs, hi)'; + if ~isempty(hi) && ~isempty(lo) + obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; + obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; + elseif ~isempty(lo) + obj.data = lofi(obj.raw', 1e6/obj.fs, lo)'; + elseif ~isempty(hi) + obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; end - if ~isempty(lo) - obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; + end + + function crop(obj, tmin, tmax) + arguments + obj + tmin {mustBeGreaterThanOrEqual(tmin, 0)} + tmax {mustBeReal} = obj.times(end) end + + start = tmin * obj.fs + 1; + stop = tmax * obj.fs + 1; + + obj.first_samp = obj.time_as_index(1) + start - 1; + obj.last_samp = obj.time_as_index(1) + stop - 1; + + obj.time_as_index = obj.time_as_index(1, ... + start:stop); + + obj.data = obj.data(:, start:stop); + times_vec = obj.times(:, start:stop); + obj.times = times_vec - min(times_vec); + end - - function [r, meta] = windowedPower(obj, time_window, noverlap, ... + + function [r, r1, meta] = windowedPower(obj, time_window, noverlap, ... lfreq, hfreq, verbose) % calculate powers in a windowed fashion arguments @@ -78,23 +101,25 @@ function filterEEG(obj,hi,lo) noverlap (1,1) double lfreq (1,1) = 8 hfreq (1,1) = 13 - verbose logical = true + verbose logical = false end - + rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... (obj.fs*time_window - obj.fs*noverlap)) + 1; - matrix = zeros(rowsNo, size(obj.data,1)); - + matrix = zeros(rowsNo, obj.no_chan); + matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); for j = 1:rowsNo begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; stop = time_window*obj.fs + ... (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - - powers = obj.sum_power_segment((begin:stop), lfreq, hfreq); + + [powers,cur_DSA] = obj.sum_power_segment((begin:stop), lfreq, hfreq); matrix(j, :) = powers; + matrix_DSA(:,:,j) = cur_DSA; end r = matrix; + r1 = matrix_DSA; meta = struct('time_window', time_window, 'noverlap', noverlap); if verbose @@ -105,7 +130,8 @@ function filterEEG(obj,hi,lo) function r = sum_power(obj, l_freq, h_freq) r = sum_freq_band(obj.psd, obj.freq, l_freq, h_freq); end - function r = sum_power_segment(obj, segment, l_freq, h_freq) + + function [r,r1] = sum_power_segment(obj, segment, l_freq, h_freq) transposed = obj.data'; % segment in datapoints as ``double`` art_mat = abs(transposed(segment,:)) > 100; @@ -114,33 +140,13 @@ function filterEEG(obj,hi,lo) freq_window = NaN(129,1); r = NaN(1,obj.no_chan); if sum(art_vec) > 0 - [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... - [],[],256,obj.fs); - - r = sum_freq_band(psd_window, freq_window, l_freq, h_freq); + [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... + [],[],256,obj.fs); + r = sum_freq_band(psd_window, freq_window, l_freq, h_freq); end + r1 = psd_window; end - function crop(obj, tmin, tmax) - arguments - obj - tmin {mustBeGreaterThanOrEqual(tmin, 0)} - tmax {mustBeReal} = obj.times(end) - end - - start = tmin * obj.fs + 1; - stop = tmax * obj.fs + 1; - obj.first_samp = obj.time_as_index(1) + start - 1; - obj.last_samp = obj.time_as_index(1) + stop - 1; - - obj.time_as_index = obj.time_as_index(1, ... - start:stop); - - obj.data = obj.data(:, start:stop); - times_vec = obj.times(:, start:stop); - obj.times = times_vec - min(times_vec); - - end function r = apply_function(obj, func, channel_wise, inplace) arguments obj @@ -155,6 +161,7 @@ function crop(obj, tmin, tmax) r = obj.apply_for_signal(func, channel_wise); end end + function r = plot_channels_3d(obj, return_figure) arguments obj @@ -183,7 +190,7 @@ function crop(obj, tmin, tmax) time_step = 1/obj.fs; endpoint = length(obj.data)/obj.fs; r = [time_step:time_step:endpoint]; -% r = reshape(double( 1/obj.fs : length(obj.data) ) * 1/obj.fs, 1, []); + % r = reshape(double( 1/obj.fs : length(obj.data) ) * 1/obj.fs, 1, []); end function r = apply_for_signal(obj, func, channel_wise) diff --git a/pvpmod.m b/pvpmod.m new file mode 100644 index 0000000..1c4c20d --- /dev/null +++ b/pvpmod.m @@ -0,0 +1,40 @@ +function pvpmod(x,parNames) +% PVPMOD - evaluate parameter/value pairs +% pvpmod(x,parNames) assigns the value x(i+1) to the parameter defined by +% the string x(i) in the calling workspace. This is useful to evaluate +% contents in an mfile, e.g. to change default settings of any +% variable initialized before pvpmod(x) is called. +% Before these assignments the code checks whether the inputs are arranged +% in pairs with the first element of each pair being a char or string. If +% optional input argument parNames (a cell array of char arrays) is +% specified, it also ensures that each of the parameter names in x +% are part of the allowable set given in parNames. + + +% (c) U. Egert 1998 +% Modified by H. Hentschke 2017 + +% new: +% 0. check inputs +if nargin>1 + validateattributes(parNames,{'cell'},{'nonempty'}); +end +if ~isempty(x) + % 0. check inputs + validateattributes(x,{'cell'},{'nonempty'}); + % 1. check that we're really dealing with pairs in x (i.e. varargin has + % an even number of elements) + if rem(numel(x),2)~=0 + error('expecting parameter/value pairs') + end + for k = 1:2:numel(x) + % 2. make sure the first input within a pair is a char or string + validateattributes(x{k},{'char','string'},{'nonempty'}); + if nargin>1 + % 3. make sure strings match + validatestring(x{k},parNames); + end + % now make assignments + assignin('caller', x{k}, x{k+1}); + end +end \ No newline at end of file From b1d66a117ae48f3a9ba1c6914636a459e76c36e5 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Fri, 15 Apr 2022 15:20:14 +0200 Subject: [PATCH 10/21] Changes to BaseRaw Introduce test functions which basically do the same however give us more flexibility. Start a simple plotting function --- BaseRaw.m | 80 +++++++++++++++++++++++++++++++++++++++++++++++-- sum_freq_band.m | 5 ++++ 2 files changed, 83 insertions(+), 2 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 0c209c9..f749e29 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -18,7 +18,7 @@ time_as_index id = 1 end - +%% Methods for Basic EEG Processing and Visualization methods function obj = BaseRaw(EEG, picks) arguments @@ -60,6 +60,7 @@ end function filterEEG(obj,hi,lo) + % Apply Basic High and Low Pass filter if ~isempty(hi) && ~isempty(lo) obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; @@ -71,6 +72,8 @@ function filterEEG(obj,hi,lo) end function crop(obj, tmin, tmax) + % Crop at specific timepoint + % To Do: Crop at index, not at second arguments obj tmin {mustBeGreaterThanOrEqual(tmin, 0)} @@ -147,6 +150,80 @@ function crop(obj, tmin, tmax) r1 = psd_window; end + %Test functions + function matrix_DSA = create_DSA(obj,time_window, noverlap) + %Creates DSA of Power Spectrum for window and noverlap + rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... + (obj.fs*time_window - obj.fs*noverlap)) + 1; + matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); + for j = 1:rowsNo + begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + stop = time_window*obj.fs + ... + (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; + + [powers] = obj.power_segment((begin:stop)); + matrix_DSA(:,:,j) = powers; + end + + end + function [r, meta] = windowedPower1(obj, time_window, noverlap, ... + lfreq, hfreq, verbose) + % calculate powers in a windowed fashion + arguments + obj + time_window (1,1) double + noverlap (1,1) double + lfreq (1,1) = 8 + hfreq (1,1) = 13 + verbose logical = false + end + + matrix_DSA = obj.create_DSA(time_window, noverlap); + r = sum_freq_band(matrix_DSA, obj.freq, lfreq, hfreq); + meta = struct('time_window', time_window, 'noverlap', noverlap); + + if verbose + disp(meta); + end + end + + function [psd_window] = power_segment(obj, segment) + transposed = obj.data'; + % segment in datapoints as ``double`` + art_mat = abs(transposed(segment,:)) > 100; + art_vec = ~any(art_mat); + psd_window = NaN(129,obj.no_chan); + if sum(art_vec) > 0 + [psd_window(:,art_vec), ~] = pwelch(transposed(segment, art_vec), ... + [],[],256,obj.fs); + end + end + + function plt = plot_single_channel(obj,ch,time_window, noverlap) + % Plot single channel DSA and Raw EEG + % To Do: + % Display Channel name, check time_vector and position of + % colorbar + plt = figure('WindowState','maximized'); + % Calculate DSA + matrix_DSA = create_DSA(obj,time_window, noverlap); + subplot(2,1,1) + xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,ch,:)))) + shading flat + shading interp + colormap turbo + colorbar + caxis([-20 20]) + ylim([0 47]) + xlim([xax(1),xax(end)]) + ylabel('Frequency [Hz]') + xlabel('Time [s]') + subplot(2,1,2) + plot(obj.times,obj.data(ch,:)) + end + +% Test over function r = apply_function(obj, func, channel_wise, inplace) arguments obj @@ -190,7 +267,6 @@ function crop(obj, tmin, tmax) time_step = 1/obj.fs; endpoint = length(obj.data)/obj.fs; r = [time_step:time_step:endpoint]; - % r = reshape(double( 1/obj.fs : length(obj.data) ) * 1/obj.fs, 1, []); end function r = apply_for_signal(obj, func, channel_wise) diff --git a/sum_freq_band.m b/sum_freq_band.m index 76a6a4d..7081c4f 100644 --- a/sum_freq_band.m +++ b/sum_freq_band.m @@ -2,4 +2,9 @@ assert(size(psd, 1) == size(freq, 1), "psd and freqs dims don't match!") idxs = freq > l_freq & freq < h_freq; +if size(psd,3) == 1 datavec = sum(psd(idxs, :), 1); +elseif size(psd,3) > 1 % Works on 3D matrix + datavec = sum(psd(idxs, :,:), 1); + datavec = squeeze(datavec)'; +end \ No newline at end of file From b4a935723ba0722e848baf364273f88f0b043aab Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Sat, 16 Apr 2022 17:43:16 +0200 Subject: [PATCH 11/21] Add new functions --- BaseRaw.m | 117 ++++++++++++++++++++++++++++++++++++++++++++--- plot_mult_topo.m | 14 ++++-- sum_freq_band.m | 10 ---- 3 files changed, 122 insertions(+), 19 deletions(-) delete mode 100644 sum_freq_band.m diff --git a/BaseRaw.m b/BaseRaw.m index f749e29..8b998b0 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -199,16 +199,31 @@ function crop(obj, tmin, tmax) end end - function plt = plot_single_channel(obj,ch,time_window, noverlap) + function datavec = sum_freq_band(obj, psd, l_freq, h_freq) + assert(size(psd, 1) == size(obj.freq, 1), "psd and freqs dims don't match!") + + idxs = obj.freq >= l_freq & obj.freq < h_freq; + if size(psd,3) == 1 + datavec = sum(psd(idxs, :), 1); + elseif size(psd,3) > 1 % Works on 3D matrix + datavec = sum(psd(idxs, :,:), 1); + datavec = squeeze(datavec)'; + end + end + + function plt = plot_single_channel(obj,ch,time_window, noverlap,xax) % Plot single channel DSA and Raw EEG % To Do: - % Display Channel name, check time_vector and position of - % colorbar - plt = figure('WindowState','maximized'); + % check time_vector and position of the colorbar + plt = figure('WindowState','maximized','Name',obj.subject); + ch_name = obj.chanlocs(ch).labels; + if isempty(xax) + xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + end + sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') % Calculate DSA matrix_DSA = create_DSA(obj,time_window, noverlap); subplot(2,1,1) - xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,ch,:)))) shading flat shading interp @@ -218,9 +233,99 @@ function crop(obj, tmin, tmax) ylim([0 47]) xlim([xax(1),xax(end)]) ylabel('Frequency [Hz]') - xlabel('Time [s]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') subplot(2,1,2) plot(obj.times,obj.data(ch,:)) + ylabel('Amplitude [mV]') + xlabel('Time [s]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') + end + + function tg = plot_all_channels(obj,time_window, noverlap,xax) + % Plot channels DSA and Raw EEG + % To Do: + % Display Channel name, check time_vector and position of + % colorbar + if isempty(xax) + xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + end + plt = figure('WindowState','maximized','Name',obj.subject); + % Calculate DSA + matrix_DSA = create_DSA(obj,time_window, noverlap); + tg = uitabgroup(plt); % tabgroup + for i = 1:obj.no_chan + ch_name = obj.chanlocs(i).labels; + thistab = uitab(tg,"Title",ch_name); % build iith tab + axes('Parent',thistab); % somewhere to plot + sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') + subplot(2,1,1) + pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,i,:)))) + shading flat + shading interp + colormap turbo + colorbar + caxis([-20 20]) + ylim([0 47]) + xlim([xax(1),xax(end)]) + ylabel('Frequency [Hz]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') + subplot(2,1,2) + plot(obj.times,obj.data(i,:)) + ylabel('Amplitude [mV]') + xlabel('Time [s]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') + end + end + + function tg = plot_power_bands(obj,matrix_DSA,time_window,noverlap,xax) + % User can supply DSA_matrix (e.g. Averaged DSA over all + % patients), if nothing is supplied function operates on single + % patient basis + % X-Axis can be supplied,otherwise xax will be calculated based + % on window and noverlap + if isempty(xax) + xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + end + if isempty(matrix_DSA) + matrix_DSA = create_DSA(obj,time_window, noverlap); + end + plt = figure('WindowState','maximized','Name',obj.subject); + l = sum_freq_band(matrix_DSA, obj.freq, 1, 6.99); + a = sum_freq_band(matrix_DSA, obj.freq, 7, 12.99); + b = sum_freq_band(matrix_DSA, obj.freq, 13, 29.99); + g = sum_freq_band(matrix_DSA, obj.freq, 30, 46); + % Create 4 plots with bandpower over time + tg = uitabgroup(plt); % tabgroup + sz = 5; + for i = 1:obj.no_chan + ch_name = obj.chanlocs(i).labels; + thistab = uitab(tg,"Title",ch_name); % build iith tab + axes('Parent',thistab); % somewhere to plot + sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') + ax(1) = subplot(2,2,1); + scatter(xax,l(:,i),sz,'filled') + ylabel('Power') + title('Low Frequency Band') + ax(2) = subplot(2,2,2); + scatter(xax,a(:,i),sz,'filled') + ylabel('Power') + title('Alpha Band') + ax(3) = subplot(2,2,3); + scatter(xax,b(:,i),sz,'filled') + ylabel('Power') + title('Beta Band') + ax(4) = subplot(2,2,4); + scatter(xax,g(:,i),sz,'filled') + ylabel('Power') + title('Gamma Band') + set(ax,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(ax,'off') + xlim(ax,[xax(1),xax(end)]) + end end % Test over diff --git a/plot_mult_topo.m b/plot_mult_topo.m index 6360b86..e061476 100644 --- a/plot_mult_topo.m +++ b/plot_mult_topo.m @@ -32,9 +32,17 @@ function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) for i = 1:framesNo pb.print(i,framesNo) for j = 1:num_plots - ax(j) = subplot(rows,cols,j); - sgtitle([num2str(time_vec(i)), ' s']) - title(desc{1,j}); + if num_plots == 1 + t1 = [num2str(time_vec(i)), ' s']; + t2 = desc; + t = [t1;t2]; + title(t); + else + ax(j) = subplot(rows,cols,j); + sgtitle([num2str(time_vec(i)), ' s']) + title(desc{1,j}); + + end topoplot(mat(i,:,j),chans); colorbar caxis([lims(j,1),lims(j,2)]) diff --git a/sum_freq_band.m b/sum_freq_band.m deleted file mode 100644 index 7081c4f..0000000 --- a/sum_freq_band.m +++ /dev/null @@ -1,10 +0,0 @@ -function datavec = sum_freq_band(psd, freq, l_freq, h_freq) -assert(size(psd, 1) == size(freq, 1), "psd and freqs dims don't match!") - -idxs = freq > l_freq & freq < h_freq; -if size(psd,3) == 1 -datavec = sum(psd(idxs, :), 1); -elseif size(psd,3) > 1 % Works on 3D matrix - datavec = sum(psd(idxs, :,:), 1); - datavec = squeeze(datavec)'; -end \ No newline at end of file From e4f1c4d8d5517851f65a3404ec2e58a1016886fd Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Sun, 17 Apr 2022 10:38:55 +0200 Subject: [PATCH 12/21] Add artefact matrix plot --- BaseRaw.m | 45 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/BaseRaw.m b/BaseRaw.m index 8b998b0..a047704 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -328,7 +328,50 @@ function crop(obj, tmin, tmax) end end -% Test over + function fig = plot_artefact_matrix(obj,time_window, threshold,time_vec) + + size_data = size(obj.data,2); + no_epochs = floor(size_data(1)/(time_window*obj.fs)); + M_woA = zeros(obj.no_chan, no_epochs); % Matrix to contain time series without artefacts + if isempty(time_vec) + time_vec = 1:time_window:no_epochs*time_window; + end + for ch = 1:obj.no_chan % loop over (the 64) channels + + channel_data = obj.data(ch, :); % data set for channel under consideration + + for j = 1:no_epochs % loop over (the 4s) epochs + + first_sample = (j-1)*time_window*obj.fs + 1; + last_sample = j*time_window*obj.fs; + epoch_data = abs(channel_data(first_sample:last_sample)); + artifacts = epoch_data >= threshold; + if (any(artifacts)) % set values for Matrix M = 1 if artifacts found + M_woA(ch, j) = 1; + end + end + end + + % Plot matrix M + figure('WindowState','maximized','Name',obj.subject) + [r, c] = size(M_woA); % Get the matrix size + imagesc((1:c)+0.5, (1:r)+0.5, M_woA); % Plot the image + hold on; + colormap(gray); % Use a gray colormap + xlabel('Number of Epoch') + axis equal + for i = 1:obj.no_chan + locs(i) = {obj.chanlocs(i).labels}; + end + %y_labels = {channels{1}:channels{r}}; + set(gca,'XTick', 1:2:c, 'YTick', 1:r, ... % Change axes properties + 'YTickLabel', locs,... + 'XLim', [1 c+1], 'YLim', [1 r+1], ... + 'GridLineStyle', '-', 'XGrid', 'on', 'YGrid', 'on','FontSize',6); + hold off + end + + % Test over function r = apply_function(obj, func, channel_wise, inplace) arguments obj From 11ca648688e04e47e96bddb92d3de83062bad343 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Sun, 17 Apr 2022 17:18:49 +0200 Subject: [PATCH 13/21] Integrate sum_freq_band to BaseRaw --- BaseRaw.m | 14 +++++++------- plot_mult_topo.m | 3 ++- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index a047704..d2e8e86 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -131,7 +131,7 @@ function crop(obj, tmin, tmax) end function r = sum_power(obj, l_freq, h_freq) - r = sum_freq_band(obj.psd, obj.freq, l_freq, h_freq); + r = obj.sum_freq_band(obj.psd, l_freq, h_freq); end function [r,r1] = sum_power_segment(obj, segment, l_freq, h_freq) @@ -145,7 +145,7 @@ function crop(obj, tmin, tmax) if sum(art_vec) > 0 [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... [],[],256,obj.fs); - r = sum_freq_band(psd_window, freq_window, l_freq, h_freq); + r = obj.sum_freq_band(psd_window, l_freq, h_freq); end r1 = psd_window; end @@ -179,7 +179,7 @@ function crop(obj, tmin, tmax) end matrix_DSA = obj.create_DSA(time_window, noverlap); - r = sum_freq_band(matrix_DSA, obj.freq, lfreq, hfreq); + r = obj.sum_freq_band(matrix_DSA, lfreq, hfreq); meta = struct('time_window', time_window, 'noverlap', noverlap); if verbose @@ -294,10 +294,10 @@ function crop(obj, tmin, tmax) matrix_DSA = create_DSA(obj,time_window, noverlap); end plt = figure('WindowState','maximized','Name',obj.subject); - l = sum_freq_band(matrix_DSA, obj.freq, 1, 6.99); - a = sum_freq_band(matrix_DSA, obj.freq, 7, 12.99); - b = sum_freq_band(matrix_DSA, obj.freq, 13, 29.99); - g = sum_freq_band(matrix_DSA, obj.freq, 30, 46); + l = obj.sum_freq_band(matrix_DSA, 1, 6.99); + a = obj.sum_freq_band(matrix_DSA, 7, 12.99); + b = obj.sum_freq_band(matrix_DSA, 13, 29.99); + g = obj.sum_freq_band(matrix_DSA, 30, 46); % Create 4 plots with bandpower over time tg = uitabgroup(plt); % tabgroup sz = 5; diff --git a/plot_mult_topo.m b/plot_mult_topo.m index e061476..11d0c80 100644 --- a/plot_mult_topo.m +++ b/plot_mult_topo.m @@ -28,7 +28,7 @@ function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) pb = CmdLineProgressBar(("Preparing "+framesNo+ " Frames... ")); %% Create Figure and axes -figure('WindowState','maximized','Visible','off') +figure('WindowState','maximized','Visible','on') for i = 1:framesNo pb.print(i,framesNo) for j = 1:num_plots @@ -47,6 +47,7 @@ function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) colorbar caxis([lims(j,1),lims(j,2)]) F(i) = getframe(gcf); + cla(gca) end end From e8263339c20e0416983f47de39fb5000941cdb15 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Tue, 19 Apr 2022 13:37:26 +0200 Subject: [PATCH 14/21] Bug Fixes of overlap calculation There was a small mistake in the calculation of noverlap and the timewindows but it should work now --- BaseRaw.m | 29 ++++++++++++++++------------- plot_mult_topo.m | 3 +++ 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index d2e8e86..af4ead0 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -7,6 +7,7 @@ freq raw no_chan + DSA end properties (SetAccess = private) fs @@ -154,18 +155,18 @@ function crop(obj, tmin, tmax) function matrix_DSA = create_DSA(obj,time_window, noverlap) %Creates DSA of Power Spectrum for window and noverlap rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... - (obj.fs*time_window - obj.fs*noverlap)) + 1; + (obj.fs*noverlap)) + 1; matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); - for j = 1:rowsNo - begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - stop = time_window*obj.fs + ... - (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - + cnt = 1; + for j=1:obj.fs*noverlap:length(obj.data)-time_window*obj.fs + begin = j; + stop = j+time_window*obj.fs-1; [powers] = obj.power_segment((begin:stop)); - matrix_DSA(:,:,j) = powers; - end - + matrix_DSA(:,:,cnt) = powers; + cnt = cnt+1; + end end + function [r, meta] = windowedPower1(obj, time_window, noverlap, ... lfreq, hfreq, verbose) % calculate powers in a windowed fashion @@ -218,7 +219,7 @@ function crop(obj, tmin, tmax) plt = figure('WindowState','maximized','Name',obj.subject); ch_name = obj.chanlocs(ch).labels; if isempty(xax) - xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + xax = 0:noverlap:size(obj.data,2)/obj.fs-time_window; end sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') % Calculate DSA @@ -243,17 +244,19 @@ function crop(obj, tmin, tmax) box(gca,'off') end - function tg = plot_all_channels(obj,time_window, noverlap,xax) + function tg = plot_all_channels(obj,matrix_DSA,time_window, noverlap,xax) % Plot channels DSA and Raw EEG % To Do: % Display Channel name, check time_vector and position of % colorbar if isempty(xax) - xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + xax = 0:noverlap:size(obj.data,2)/obj.fs-time_window; end plt = figure('WindowState','maximized','Name',obj.subject); % Calculate DSA + if isempty(matrix_DSA) matrix_DSA = create_DSA(obj,time_window, noverlap); + end tg = uitabgroup(plt); % tabgroup for i = 1:obj.no_chan ch_name = obj.chanlocs(i).labels; @@ -288,7 +291,7 @@ function crop(obj, tmin, tmax) % X-Axis can be supplied,otherwise xax will be calculated based % on window and noverlap if isempty(xax) - xax = 0:time_window-noverlap:size(obj.data,2)/obj.fs-time_window; + xax = 0:noverlap:size(obj.data,2)/obj.fs-time_window; end if isempty(matrix_DSA) matrix_DSA = create_DSA(obj,time_window, noverlap); diff --git a/plot_mult_topo.m b/plot_mult_topo.m index 11d0c80..660b845 100644 --- a/plot_mult_topo.m +++ b/plot_mult_topo.m @@ -32,6 +32,7 @@ function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) for i = 1:framesNo pb.print(i,framesNo) for j = 1:num_plots + if num_plots == 1 t1 = [num2str(time_vec(i)), ' s']; t2 = desc; @@ -47,7 +48,9 @@ function plot_mult_topo(mat,chans,desc,lims,time_vec,fname,fps) colorbar caxis([lims(j,1),lims(j,2)]) F(i) = getframe(gcf); + if num_plots == 1 cla(gca) + end end end From 0014692884fd44a855ba0ca7788a67a9fe1d11ca Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Wed, 20 Apr 2022 08:21:18 +0200 Subject: [PATCH 15/21] Update BaseRaw.m --- BaseRaw.m | 91 ++++++++----------------------------------------------- 1 file changed, 12 insertions(+), 79 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index af4ead0..d6090d5 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -12,8 +12,6 @@ properties (SetAccess = private) fs chanlocs - first_samp = 1 - last_samp end properties (Hidden = true) time_as_index @@ -28,7 +26,7 @@ end if nargin >= 1 - obj.data = EEG.data(picks, :); + obj.data = EEG.data(picks, :); % Rows = channels obj.raw = EEG.data(picks, :); obj.fs = EEG.srate; obj.no_chan = size(obj.data,1); @@ -54,7 +52,6 @@ else obj.chanlocs = []; end - obj.last_samp = length(obj); obj.time_as_index = (1 : length(obj)); end @@ -62,6 +59,7 @@ function filterEEG(obj,hi,lo) % Apply Basic High and Low Pass filter + % Filter works on columns if ~isempty(hi) && ~isempty(lo) obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; @@ -72,86 +70,23 @@ function filterEEG(obj,hi,lo) end end - function crop(obj, tmin, tmax) - % Crop at specific timepoint - % To Do: Crop at index, not at second + function crop(obj, start, stop) + % Crop at specific index arguments obj - tmin {mustBeGreaterThanOrEqual(tmin, 0)} - tmax {mustBeReal} = obj.times(end) + start {mustBeGreaterThanOrEqual(start, 1)} + stop {mustBeReal} = obj.times(end) end - start = tmin * obj.fs + 1; - stop = tmax * obj.fs + 1; - - obj.first_samp = obj.time_as_index(1) + start - 1; - obj.last_samp = obj.time_as_index(1) + stop - 1; - - obj.time_as_index = obj.time_as_index(1, ... - start:stop); - obj.data = obj.data(:, start:stop); times_vec = obj.times(:, start:stop); obj.times = times_vec - min(times_vec); - - end - - function [r, r1, meta] = windowedPower(obj, time_window, noverlap, ... - lfreq, hfreq, verbose) - % calculate powers in a windowed fashion - arguments - obj - time_window (1,1) double - noverlap (1,1) double - lfreq (1,1) = 8 - hfreq (1,1) = 13 - verbose logical = false - end - - rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... - (obj.fs*time_window - obj.fs*noverlap)) + 1; - matrix = zeros(rowsNo, obj.no_chan); - matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); - for j = 1:rowsNo - begin = (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - stop = time_window*obj.fs + ... - (time_window*obj.fs - noverlap*obj.fs) * (j-1) + 1; - - [powers,cur_DSA] = obj.sum_power_segment((begin:stop), lfreq, hfreq); - matrix(j, :) = powers; - matrix_DSA(:,:,j) = cur_DSA; - end - - r = matrix; - r1 = matrix_DSA; - meta = struct('time_window', time_window, 'noverlap', noverlap); - - if verbose - disp(meta); - end end function r = sum_power(obj, l_freq, h_freq) r = obj.sum_freq_band(obj.psd, l_freq, h_freq); end - function [r,r1] = sum_power_segment(obj, segment, l_freq, h_freq) - transposed = obj.data'; - % segment in datapoints as ``double`` - art_mat = abs(transposed(segment,:)) > 100; - art_vec = ~any(art_mat); - psd_window = NaN(129,obj.no_chan); - freq_window = NaN(129,1); - r = NaN(1,obj.no_chan); - if sum(art_vec) > 0 - [psd_window(:,art_vec), freq_window] = pwelch(transposed(segment, art_vec), ... - [],[],256,obj.fs); - r = obj.sum_freq_band(psd_window, l_freq, h_freq); - end - r1 = psd_window; - end - - %Test functions function matrix_DSA = create_DSA(obj,time_window, noverlap) %Creates DSA of Power Spectrum for window and noverlap rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... @@ -167,7 +102,7 @@ function crop(obj, tmin, tmax) end end - function [r, meta] = windowedPower1(obj, time_window, noverlap, ... + function [r, meta] = windowedPower(obj, time_window, noverlap, ... lfreq, hfreq, verbose) % calculate powers in a windowed fashion arguments @@ -193,7 +128,8 @@ function crop(obj, tmin, tmax) % segment in datapoints as ``double`` art_mat = abs(transposed(segment,:)) > 100; art_vec = ~any(art_mat); - psd_window = NaN(129,obj.no_chan); + no_rows = size(obj.psd,1); + psd_window = NaN(no_rows,obj.no_chan); if sum(art_vec) > 0 [psd_window(:,art_vec), ~] = pwelch(transposed(segment, art_vec), ... [],[],256,obj.fs); @@ -331,14 +267,11 @@ function crop(obj, tmin, tmax) end end - function fig = plot_artefact_matrix(obj,time_window, threshold,time_vec) + function fig = plot_artefact_matrix(obj,time_window, threshold) size_data = size(obj.data,2); no_epochs = floor(size_data(1)/(time_window*obj.fs)); M_woA = zeros(obj.no_chan, no_epochs); % Matrix to contain time series without artefacts - if isempty(time_vec) - time_vec = 1:time_window:no_epochs*time_window; - end for ch = 1:obj.no_chan % loop over (the 64) channels channel_data = obj.data(ch, :); % data set for channel under consideration @@ -356,7 +289,7 @@ function crop(obj, tmin, tmax) end % Plot matrix M - figure('WindowState','maximized','Name',obj.subject) + fig = figure('WindowState','maximized','Name',obj.subject); [r, c] = size(M_woA); % Get the matrix size imagesc((1:c)+0.5, (1:r)+0.5, M_woA); % Plot the image hold on; @@ -417,7 +350,7 @@ function crop(obj, tmin, tmax) function r = create_times(obj) time_step = 1/obj.fs; endpoint = length(obj.data)/obj.fs; - r = [time_step:time_step:endpoint]; + r = time_step:time_step:endpoint; end function r = apply_for_signal(obj, func, channel_wise) From 378d02c43629cabfc3aeed689eed8b08a41d1cc9 Mon Sep 17 00:00:00 2001 From: dnaligase <60392997+dnaligase@users.noreply.github.com> Date: Fri, 22 Apr 2022 16:33:40 +0200 Subject: [PATCH 16/21] Create connectivity.m --- connectivity.m | 104 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 104 insertions(+) create mode 100644 connectivity.m diff --git a/connectivity.m b/connectivity.m new file mode 100644 index 0000000..8eddbe4 --- /dev/null +++ b/connectivity.m @@ -0,0 +1,104 @@ +fnames = dir('*.mat'); +% retrieve size of loading data w/o loading it to memory +data_size = whos('-file', fnames(end).name).size(); +matrix = zeros(length(fnames), data_size(1), data_size(2)); +%% +ntrials = 0; % number of trials +nobs = size(matrix,3); % number of observations per trial +regmode = 'LWR'; % VAR model estimation regression mode ('OLS', 'LWR' or empty for default) +icregmode = 'LWR'; % information criteria regression mode ('OLS', 'LWR' or empty for default) +morder = 'AIC'; % model order to use ('actual', 'AIC', 'BIC' or supplied numerical value) +momax = 3; % maximum model order for model order estimation +acmaxlags = []; % maximum autocovariance lags (empty for automatic calculation) +tstat = 'F'; % statistical test for MVGC: 'F' for Granger's F-test (default) or 'chi2' for Geweke's chi2 test +alpha = 0.05; % significance level for significance test +mhtc = 'FDR'; % multiple hypothesis test correction (see routine 'significance') +fs = 250; % sample rate (Hz) +fres = 250; % frequency resolution (empty for automatic calculation) +seed = 0; + +bands = struct('delta', {1, 4}, 'alpha', {7, 13}, 'beta', {13, 30}); +%% +% SET PARAMS +event_time = 120; % in seconds +patient_select = [13]; +from_electrodes = [(2:4), 18]; +time_window = 2; % in seconds +noverlap = 0; % in seconds +k_neighbors = 5; % k-nearest neighbors for smoothing +band = 'delta'; % choose from ['delta' | 'alpha' | 'beta'] +%% +for i = patient_select + path = fullfile(fnames(i).folder, fnames(i).name); + matrix(i, :, :) = importdata(path); +end + +rowsNo = fix((size(matrix, 3) - fs*time_window) / (fs*time_window - fs*noverlap)) + 1; +matrix = matrix(:, :, 1:size(matrix,3)-rem(size(matrix,3), time_window*fs)); +matrix = reshape(matrix, length(fnames), data_size(1), [], rowsNo); +% Z-scale our data across epochs (3rd dim) +matrix = normalize(matrix, 3); +%% +% bivariate granger +event = event_time / time_window; +traces = nan(rowsNo,1); +pairs = nchoosek(from_electrodes, 2); +lfreq = getfield(bands(1), band); +hfreq = getfield(bands(2), band); +f = figure(); + +for k = 1:size(pairs, 1) + % pick channel pair + matrix_pick = matrix(patient_select, pairs(k,:), :, :); + + slices = zeros(rowsNo, size(matrix_pick, 2)^2 - size(matrix_pick, 2)); + grangers_time = zeros(2,2,rowsNo); + Fints_time = zeros(2,2,rowsNo); + for slice_id = 1:size(matrix_pick, 4) + grangers = nan(size(matrix_pick, 2)); + Fints = nan(size(matrix_pick, 2)); + for patient_id = 1:size(matrix_pick, 1) + data = squeeze(matrix_pick(patient_id, :, :, slice_id)); + [A,SIG] = tsdata_to_var(data,momax); + + f = var_to_spwcgc(A,SIG,fres); + Fint = smvgc_to_mvgc(f, ... + [ ... + lfreq / (fs/2), ... + hfreq / (fs/2) ... + ]); + + [F, ~] = var_to_pwcgc(A,SIG,data,regmode,tstat); + grangers = cat(3, grangers, F); + Fints = cat(3, Fints, Fint); + + end + grangers_mean = mean(grangers, 3, 'omitnan'); + grangers_time(:, :, slice_id) = grangers_mean; + Fints_mean = median(Fints, 3, 'omitnan'); + Fints_time(:, :, slice_id) = Fints_mean; + end + +% PLOTS + + for i = 1:2 + for j = 1:2 + datplot = squeeze(Fints_time(j, i, :)); + % normalize + datplot_scaled = datplot / prctile(datplot(1:event), 95); + % apply moving average to smooth signal + datplot_scaled_movm = movmean(datplot_scaled, k_neighbors); + % plot log-y + semilogy(datplot_scaled_movm, Color='black') + hold on + if ~all(isnan(datplot_scaled_movm)) + traces = cat(2, traces, datplot_scaled_movm); + end + end + end +end + +xline(event) +title(band, ('across '+string(length(patient_select))+' patient(s)')) +plot(median(traces, 2, 'omitnan'), LineWidth=3, Color=[1.0000 0.6471 0]) +hold off \ No newline at end of file From 75e44a633d8e1c34537c7783bab01fffe672dd9c Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Mon, 25 Apr 2022 12:22:09 +0200 Subject: [PATCH 17/21] Just some comments --- BaseRaw.m | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index d6090d5..11903bf 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -1,17 +1,17 @@ classdef BaseRaw < handle properties - data - subject - times - psd - freq - raw - no_chan - DSA + data % Conatins current EEG signal + subject % Subject Name + times % + psd % Power Spectrum for the full EEG signal/"data" + freq % Frequency vector + no_chan % Number of channels currently selected + DSA % Density Spectral Array end properties (SetAccess = private) - fs - chanlocs + fs % Sampling Frequency + chanlocs % Channel locations + raw % Contains raw EEG, won't be changed end properties (Hidden = true) time_as_index @@ -25,10 +25,10 @@ picks = (1: size(EEG.data, 1)) end - if nargin >= 1 + if nargin >= 1 obj.data = EEG.data(picks, :); % Rows = channels - obj.raw = EEG.data(picks, :); - obj.fs = EEG.srate; + obj.raw = EEG.data(picks, :); % Stores original data + obj.fs = EEG.srate; % obj.no_chan = size(obj.data,1); if ~isfield(EEG, 'subject') || isempty(EEG.subject) From fb804f4256e494035a0c3bfbb178cf6b034a62e1 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Tue, 26 Apr 2022 09:07:15 +0200 Subject: [PATCH 18/21] some comments and deletion of unused function --- BaseRaw.m | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 11903bf..7b31b77 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -20,6 +20,9 @@ %% Methods for Basic EEG Processing and Visualization methods function obj = BaseRaw(EEG, picks) + % Create the base object and add some information about the input EEG object + % EEG.data and EEG.srate need to be supplied by the user + % Channels for analysis can be specified in the picks argument arguments EEG picks = (1: size(EEG.data, 1)) @@ -28,7 +31,7 @@ if nargin >= 1 obj.data = EEG.data(picks, :); % Rows = channels obj.raw = EEG.data(picks, :); % Stores original data - obj.fs = EEG.srate; % + obj.fs = EEG.srate; obj.no_chan = size(obj.data,1); if ~isfield(EEG, 'subject') || isempty(EEG.subject) @@ -71,7 +74,7 @@ function filterEEG(obj,hi,lo) end function crop(obj, start, stop) - % Crop at specific index + % Crop the EEG signal at a specific index arguments obj start {mustBeGreaterThanOrEqual(start, 1)} @@ -83,12 +86,8 @@ function crop(obj, start, stop) obj.times = times_vec - min(times_vec); end - function r = sum_power(obj, l_freq, h_freq) - r = obj.sum_freq_band(obj.psd, l_freq, h_freq); - end - function matrix_DSA = create_DSA(obj,time_window, noverlap) - %Creates DSA of Power Spectrum for window and noverlap + % Creates DSA of Power Spectrum for window and noverlap rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... (obj.fs*noverlap)) + 1; matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); @@ -341,6 +340,7 @@ function crop(obj, start, stop) zlabel('Z') if return_figure; r = gcf(); end end + function r = length(obj) r = length(obj.data); end From da3908987c89388dcd3941334f94b8c61af62417 Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Thu, 28 Apr 2022 12:46:26 +0200 Subject: [PATCH 19/21] Some comments --- BaseRaw.m | 68 +++++++++++++++++++++++++++++-------------------------- 1 file changed, 36 insertions(+), 32 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 7b31b77..8979938 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -19,6 +19,7 @@ end %% Methods for Basic EEG Processing and Visualization methods + % Constructor function, initializes the object function obj = BaseRaw(EEG, picks) % Create the base object and add some information about the input EEG object % EEG.data and EEG.srate need to be supplied by the user @@ -59,7 +60,8 @@ end end - + + %% Basic Low and Highpass Filter, EEG Cropping function filterEEG(obj,hi,lo) % Apply Basic High and Low Pass filter % Filter works on columns @@ -72,7 +74,8 @@ function filterEEG(obj,hi,lo) obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; end end - + + function crop(obj, start, stop) % Crop the EEG signal at a specific index arguments @@ -85,7 +88,9 @@ function crop(obj, start, stop) times_vec = obj.times(:, start:stop); obj.times = times_vec - min(times_vec); end - + + %% EEG Processing + % Create DSA matrix function matrix_DSA = create_DSA(obj,time_window, noverlap) % Creates DSA of Power Spectrum for window and noverlap rowsNo = fix((size(obj.data, 2) - obj.fs*time_window) / ... @@ -98,9 +103,11 @@ function crop(obj, start, stop) [powers] = obj.power_segment((begin:stop)); matrix_DSA(:,:,cnt) = powers; cnt = cnt+1; - end + end + obj.DSA = matrix_DSA; end - + + % Extracts Power in a specific frequency range for a specific window and overlap function [r, meta] = windowedPower(obj, time_window, noverlap, ... lfreq, hfreq, verbose) % calculate powers in a windowed fashion @@ -121,10 +128,10 @@ function crop(obj, start, stop) disp(meta); end end - + + % Calculates PSD for a defined segment with artifacts function [psd_window] = power_segment(obj, segment) transposed = obj.data'; - % segment in datapoints as ``double`` art_mat = abs(transposed(segment,:)) > 100; art_vec = ~any(art_mat); no_rows = size(obj.psd,1); @@ -134,7 +141,7 @@ function crop(obj, start, stop) [],[],256,obj.fs); end end - + % Sum power within a defines frequency range function datavec = sum_freq_band(obj, psd, l_freq, h_freq) assert(size(psd, 1) == size(obj.freq, 1), "psd and freqs dims don't match!") @@ -146,7 +153,8 @@ function crop(obj, start, stop) datavec = squeeze(datavec)'; end end - + %% Plotting + % Single Channgel DSA and Raw EEG function plt = plot_single_channel(obj,ch,time_window, noverlap,xax) % Plot single channel DSA and Raw EEG % To Do: @@ -158,7 +166,9 @@ function crop(obj, start, stop) end sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') % Calculate DSA + if isempty(obj.DSA) matrix_DSA = create_DSA(obj,time_window, noverlap); + end subplot(2,1,1) pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,ch,:)))) shading flat @@ -189,7 +199,7 @@ function crop(obj, start, stop) end plt = figure('WindowState','maximized','Name',obj.subject); % Calculate DSA - if isempty(matrix_DSA) + if isempty(matrix_DSA) && isempty(obj.DSA) matrix_DSA = create_DSA(obj,time_window, noverlap); end tg = uitabgroup(plt); % tabgroup @@ -228,14 +238,17 @@ function crop(obj, start, stop) if isempty(xax) xax = 0:noverlap:size(obj.data,2)/obj.fs-time_window; end - if isempty(matrix_DSA) + + if isempty(matrix_DSA) matrix_DSA = create_DSA(obj,time_window, noverlap); end + plt = figure('WindowState','maximized','Name',obj.subject); - l = obj.sum_freq_band(matrix_DSA, 1, 6.99); - a = obj.sum_freq_band(matrix_DSA, 7, 12.99); - b = obj.sum_freq_band(matrix_DSA, 13, 29.99); - g = obj.sum_freq_band(matrix_DSA, 30, 46); + bands(:,:,1) = obj.sum_freq_band(matrix_DSA, 1, 6.99); + bands(:,:,2) = obj.sum_freq_band(matrix_DSA, 7, 12.99); + bands(:,:,3) = obj.sum_freq_band(matrix_DSA, 13, 29.99); + bands(:,:,4) = obj.sum_freq_band(matrix_DSA, 30, 46); + names = {'Low Frequency', 'Alpha Band', 'Beta Band', 'High Frequency'}; % Create 4 plots with bandpower over time tg = uitabgroup(plt); % tabgroup sz = 5; @@ -244,22 +257,12 @@ function crop(obj, start, stop) thistab = uitab(tg,"Title",ch_name); % build iith tab axes('Parent',thistab); % somewhere to plot sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') - ax(1) = subplot(2,2,1); - scatter(xax,l(:,i),sz,'filled') + for j = 1:size(bands,3) + ax(j) = subplot(2,2,j); + scatter(xax,bands(:,i,j),sz,'filled') ylabel('Power') - title('Low Frequency Band') - ax(2) = subplot(2,2,2); - scatter(xax,a(:,i),sz,'filled') - ylabel('Power') - title('Alpha Band') - ax(3) = subplot(2,2,3); - scatter(xax,b(:,i),sz,'filled') - ylabel('Power') - title('Beta Band') - ax(4) = subplot(2,2,4); - scatter(xax,g(:,i),sz,'filled') - ylabel('Power') - title('Gamma Band') + title(names{1,j}) + end set(ax,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) box(ax,'off') xlim(ax,[xax(1),xax(end)]) @@ -292,6 +295,7 @@ function crop(obj, start, stop) [r, c] = size(M_woA); % Get the matrix size imagesc((1:c)+0.5, (1:r)+0.5, M_woA); % Plot the image hold on; + scrollplot colormap(gray); % Use a gray colormap xlabel('Number of Epoch') axis equal @@ -302,11 +306,11 @@ function crop(obj, start, stop) set(gca,'XTick', 1:2:c, 'YTick', 1:r, ... % Change axes properties 'YTickLabel', locs,... 'XLim', [1 c+1], 'YLim', [1 r+1], ... - 'GridLineStyle', '-', 'XGrid', 'on', 'YGrid', 'on','FontSize',6); + 'GridLineStyle', '-', 'XGrid', 'on', 'YGrid', 'on','FontSize',5); + hold off end - % Test over function r = apply_function(obj, func, channel_wise, inplace) arguments obj From bad914ba255afae7b62298e0586a73baa6e2806e Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Thu, 28 Apr 2022 23:12:32 +0200 Subject: [PATCH 20/21] Update BaseRaw.m --- BaseRaw.m | 76 ++++++++++++++++++++++++++++++------------------------- 1 file changed, 41 insertions(+), 35 deletions(-) diff --git a/BaseRaw.m b/BaseRaw.m index 8979938..9900302 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -64,9 +64,9 @@ %% Basic Low and Highpass Filter, EEG Cropping function filterEEG(obj,hi,lo) % Apply Basic High and Low Pass filter - % Filter works on columns + % Filter works on columns (--> Channels as cols, rows as timeseries) if ~isempty(hi) && ~isempty(lo) - obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; + obj.data = hifi(obj.raw', 1e6/obj.fs, hi)'; obj.data = lofi(obj.data', 1e6/obj.fs, lo)'; elseif ~isempty(lo) obj.data = lofi(obj.raw', 1e6/obj.fs, lo)'; @@ -88,6 +88,12 @@ function crop(obj, start, stop) times_vec = obj.times(:, start:stop); obj.times = times_vec - min(times_vec); end + + function reset_rawdata(obj) + obj.data = obj.raw; + obj.DSA = []; + obj.times = obj.create_times(); + end %% EEG Processing % Create DSA matrix @@ -155,33 +161,22 @@ function crop(obj, start, stop) end %% Plotting % Single Channgel DSA and Raw EEG - function plt = plot_single_channel(obj,ch,time_window, noverlap,xax) + function ax = plot_single_channel(obj,ch,time_window, noverlap,xax) % Plot single channel DSA and Raw EEG % To Do: % check time_vector and position of the colorbar - plt = figure('WindowState','maximized','Name',obj.subject); - ch_name = obj.chanlocs(ch).labels; if isempty(xax) xax = 0:noverlap:size(obj.data,2)/obj.fs-time_window; end + figure('WindowState','maximized','Name',obj.subject); + ch_name = obj.chanlocs(ch).labels; sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') % Calculate DSA - if isempty(obj.DSA) - matrix_DSA = create_DSA(obj,time_window, noverlap); - end - subplot(2,1,1) - pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,ch,:)))) - shading flat - shading interp - colormap turbo - colorbar - caxis([-20 20]) - ylim([0 47]) - xlim([xax(1),xax(end)]) - ylabel('Frequency [Hz]') - set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) - box(gca,'off') - subplot(2,1,2) + create_DSA(obj,time_window, noverlap); + ax(1) = subplot(2,1,1); + plot_data = squeeze(obj.DSA(:,ch,:)); + ax(1) = obj.plot_DSA(ax(1),plot_data,obj.freq,xax); + ax(2) = subplot(2,1,2); plot(obj.times,obj.data(ch,:)) ylabel('Amplitude [mV]') xlabel('Time [s]') @@ -199,8 +194,8 @@ function crop(obj, start, stop) end plt = figure('WindowState','maximized','Name',obj.subject); % Calculate DSA - if isempty(matrix_DSA) && isempty(obj.DSA) - matrix_DSA = create_DSA(obj,time_window, noverlap); + if isempty(matrix_DSA) + matrix_DSA = obj.create_DSA(time_window, noverlap); end tg = uitabgroup(plt); % tabgroup for i = 1:obj.no_chan @@ -208,18 +203,9 @@ function crop(obj, start, stop) thistab = uitab(tg,"Title",ch_name); % build iith tab axes('Parent',thistab); % somewhere to plot sgtitle(ch_name,'FontName','Arial','FontSize',12,'FontWeight','Bold') - subplot(2,1,1) - pcolor(xax,obj.freq,10*log10(squeeze(matrix_DSA(:,i,:)))) - shading flat - shading interp - colormap turbo - colorbar - caxis([-20 20]) - ylim([0 47]) - xlim([xax(1),xax(end)]) - ylabel('Frequency [Hz]') - set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) - box(gca,'off') + ax(1) = subplot(2,1,1); + plot_data = squeeze(matrix_DSA(:,i,:)); + ax(1) = obj.plot_DSA(ax(1),plot_data,obj.freq,xax); subplot(2,1,2) plot(obj.times,obj.data(i,:)) ylabel('Amplitude [mV]') @@ -349,7 +335,26 @@ function crop(obj, start, stop) r = length(obj.data); end end + methods(Static) + function axis = plot_DSA(axis,DSA_mat,freq,xax) + if isempty(xax) + xax = 1:size(DSA_mat,2); + end + pcolor(gca,xax,freq,10*log10(DSA_mat)); + shading flat + shading interp + colormap turbo + colorbar + caxis(gca,[-15 20]) + ylim(gca,[0 47]) + xlim(gca,[xax(1),xax(end)]) + ylabel('Frequency [Hz]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') + end + + end methods (Access = protected) function r = create_times(obj) time_step = 1/obj.fs; @@ -357,6 +362,7 @@ function crop(obj, start, stop) r = time_step:time_step:endpoint; end + function r = apply_for_signal(obj, func, channel_wise) if channel_wise From 0c6e93db578ce7c5c043ebc76bd3e05495550cad Mon Sep 17 00:00:00 2001 From: JulianOstertag Date: Mon, 2 May 2022 10:47:09 +0200 Subject: [PATCH 21/21] Update BaseRaw.m --- BaseRaw.m | 40 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) diff --git a/BaseRaw.m b/BaseRaw.m index 9900302..0765d30 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -226,7 +226,7 @@ function reset_rawdata(obj) end if isempty(matrix_DSA) - matrix_DSA = create_DSA(obj,time_window, noverlap); + matrix_DSA = obj.create_DSA(obj,time_window, noverlap); end plt = figure('WindowState','maximized','Name',obj.subject); @@ -255,6 +255,40 @@ function reset_rawdata(obj) end end + function plot_mult_DSA(obj,DSA_mat,freq,window,shift,xax) + arguments + obj + DSA_mat + freq + window = 2 + shift = 1 + xax = 1 + end + if isempty(DSA_mat) + DSA_mat = obj.create_DSA(window,shift); + end + if isempty(xax) + xax = 1:size(DSA_mat,3); + end + plt = figure('WindowState','maximized'); + tg = uitabgroup(plt); + num_rows = 4;num_cols = 2; + ch=1; + for num_tab = 1:size(DSA_mat,2)/(num_rows*num_cols) + tab = uitab(tg,'Title',strcat('Channel #',num2str(ch),'-',num2str(ch+num_rows*num_cols-1))); + axes('Parent',tab); + for pos = 1:num_rows*num_cols + ch_name = obj.chanlocs(ch).labels; + ax(pos) = subplot(num_rows,num_cols,pos); + plot_data = squeeze(DSA_mat(:,ch,:)); + ax(pos) = obj.plot_DSA(gca,plot_data,freq,xax); + title(ch_name) + set(gca,'FontName','Arial','FontSize',6,'FontWeight','Bold', 'LineWidth', 1) + ch = ch + 1; + end + end + end + function fig = plot_artefact_matrix(obj,time_window, threshold) size_data = size(obj.data,2); @@ -334,6 +368,8 @@ function reset_rawdata(obj) function r = length(obj) r = length(obj.data); end + + end methods(Static) @@ -354,6 +390,8 @@ function reset_rawdata(obj) box(gca,'off') end + + end methods (Access = protected) function r = create_times(obj)