diff --git a/.gitignore b/.gitignore index b853f42..431d2fa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ ch_names.mat +*.asv data +.DS_Store diff --git a/BaseRaw.m b/BaseRaw.m index 382dabd..0765d30 100644 --- a/BaseRaw.m +++ b/BaseRaw.m @@ -1,34 +1,49 @@ classdef BaseRaw < handle properties - data - subject - times - psd - freq + 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 - first_samp = 1 - last_samp + fs % Sampling Frequency + chanlocs % Channel locations + raw % Contains raw EEG, won't be changed end properties (Hidden = true) time_as_index + id = 1 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 + % Channels for analysis can be specified in the picks argument arguments EEG picks = (1: size(EEG.data, 1)) end - if nargin >= 1 - obj.data = EEG.data(picks, :); - obj.fs = EEG.srate; - obj.subject = EEG.subject; + if nargin >= 1 + obj.data = EEG.data(picks, :); % Rows = channels + obj.raw = EEG.data(picks, :); % Stores original data + obj.fs = EEG.srate; + obj.no_chan = size(obj.data,1); - 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,49 +51,286 @@ end [obj.psd, obj.freq] = pwelch(obj.data',[],[],256,obj.fs); - if ~isempty(EEG.chanlocs) - obj.chanlocs = EEG.chanlocs; + if isfield(EEG, 'chanlocs') + obj.chanlocs = EEG.chanlocs(1,picks); else obj.chanlocs = []; end - obj.last_samp = length(obj); obj.time_as_index = (1 : length(obj)); 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 (--> Channels as cols, rows as timeseries) + 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 + end + + + function crop(obj, start, stop) + % Crop the EEG signal at a specific index + arguments + obj + start {mustBeGreaterThanOrEqual(start, 1)} + stop {mustBeReal} = obj.times(end) + end - function r = sum_power(obj, l_freq, h_freq) - r = sum_freq_band(obj.psd, obj.freq, l_freq, h_freq); + obj.data = obj.data(:, start:stop); + times_vec = obj.times(:, start:stop); + obj.times = times_vec - min(times_vec); end - 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, :), ... - [],[],256,obj.fs); - - r = sum_freq_band(psd_window, freq_window, l_freq, h_freq); + + function reset_rawdata(obj) + obj.data = obj.raw; + obj.DSA = []; + obj.times = obj.create_times(); + 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) / ... + (obj.fs*noverlap)) + 1; + matrix_DSA = zeros(length(obj.freq), obj.no_chan,rowsNo); + 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(:,:,cnt) = powers; + cnt = cnt+1; + end + obj.DSA = matrix_DSA; end - function crop(obj, tmin, tmax) + + % 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 arguments obj - tmin {mustBeGreaterThanOrEqual(tmin, 0)} - tmax {mustBeReal} + time_window (1,1) double + noverlap (1,1) double + lfreq (1,1) = 8 + hfreq (1,1) = 13 + verbose logical = false end - start = tmin * obj.fs + 1; - stop = tmax * obj.fs + 1; + matrix_DSA = obj.create_DSA(time_window, noverlap); + r = obj.sum_freq_band(matrix_DSA, lfreq, hfreq); + meta = struct('time_window', time_window, 'noverlap', noverlap); - obj.first_samp = obj.time_as_index(1) + start - 1; - obj.last_samp = obj.time_as_index(1) + stop - 1; + if verbose + disp(meta); + end + end + + % Calculates PSD for a defined segment with artifacts + function [psd_window] = power_segment(obj, segment) + transposed = obj.data'; + art_mat = abs(transposed(segment,:)) > 100; + art_vec = ~any(art_mat); + 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); + 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!") - obj.time_as_index = obj.time_as_index(1, ... - start:stop); + 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 + %% Plotting + % Single Channgel DSA and Raw EEG + 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 + 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 + 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]') + set(gca,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(gca,'off') + end - obj.data = obj.data(:, start:stop); - times_vec = obj.times(:, start:stop); - obj.times = times_vec - min(times_vec); + 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: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 = obj.create_DSA(time_window, noverlap); + end + 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') + 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]') + 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:noverlap:size(obj.data,2)/obj.fs-time_window; + end + if isempty(matrix_DSA) + matrix_DSA = obj.create_DSA(obj,time_window, noverlap); + end + + plt = figure('WindowState','maximized','Name',obj.subject); + 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; + 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') + for j = 1:size(bands,3) + ax(j) = subplot(2,2,j); + scatter(xax,bands(:,i,j),sz,'filled') + ylabel('Power') + title(names{1,j}) + end + set(ax,'FontName','Arial','FontSize',12,'FontWeight','Bold', 'LineWidth', 1) + box(ax,'off') + xlim(ax,[xax(1),xax(end)]) + 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); + 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 + 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 + 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; + scrollplot + 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',5); + + hold off end + function r = apply_function(obj, func, channel_wise, inplace) arguments obj @@ -93,29 +345,62 @@ 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 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'); + 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) title(('3D channel location for #' + string(obj.subject))) xlabel('X') ylabel('Y') zlabel('Z') if return_figure; r = gcf(); end end + function r = length(obj) 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) - 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; end + function r = apply_for_signal(obj, func, channel_wise) if channel_wise 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', ... 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 diff --git a/plot_mult_topo.m b/plot_mult_topo.m new file mode 100644 index 0000000..660b845 --- /dev/null +++ b/plot_mult_topo.m @@ -0,0 +1,72 @@ +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 +% 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... ")); + +%% Create Figure and axes +figure('WindowState','maximized','Visible','on') +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; + 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)]) + F(i) = getframe(gcf); + if num_plots == 1 + cla(gca) + end + 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 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 diff --git a/sum_freq_band.m b/sum_freq_band.m deleted file mode 100644 index 76a6a4d..0000000 --- a/sum_freq_band.m +++ /dev/null @@ -1,5 +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; -datavec = sum(psd(idxs, :), 1);