diff --git a/Measures/FOOOFer.m b/Measures/FOOOFer.m index 9ab5bea..794b2c2 100644 --- a/Measures/FOOOFer.m +++ b/Measures/FOOOFer.m @@ -1,199 +1,201 @@ -%% FOOOFer -% This function parameterizes the neural power spectra using the Haller et -% al. method. -% -% FOOOFer(fs, cf, nEpochs, dt, inDir, tStart, outTypes) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch -% outTypes is the list of variables to save (offset, exponent, -% peak parameters, gaussian parameters, error, r squared, fooofed -% spectrum, bg fit, power spectrum), writen in an array as strings -% in the same order as their output directories in the variable -% outDirs (es. outDirs=["C:\offset\", "C:\exponent", "C:\bg_fit"], -% outTypes=["offset", "exponent", "bg fit"]) -% -% NB: peak and gaussian params will contain a matrix where every group of 3 -% numbers are relative to one peak (CF, Amp, BW) and the others zeros -% are utilized to export only one matrix for each subject - - -function FOOOFer(fs, cf, nEpochs, dt, inDir, tStart, outTypes, maxPeaks) - switch nargin - case 5 - tStart = 0; - outTypes = ""; - maxPeaks = (cf(end)-cf(1))*2; - case 6 - outTypes = ""; - maxPeaks = (cf(end)-cf(1))*2; - case 7 - maxPeaks = (cf(end)-cf(1))*2; - end - - f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - - dt = fs*dt; - tStart = tStart*fs+1; - - askList = 'Insert the background measures separated by a comma'; - - inDir = path_check(inDir); - cases = define_cases(inDir); - outTypes = string(outTypes); - - if strcmp(outTypes, "") - outTypes = []; - sup = string(inputdlg(askList)); - if stcmp(sup, "exit") - return - end - sup = split(sup,","); - for i = 1:length(sup) - if contains(sup(i, 1), ["Exp", "exp", "EXP"]) - outTypes = [outTypes, "exp"]; - elseif contains(sup(i, 1), ["off", "OFF", "Off"]) - outTypes = [outTypes, "off"]; - elseif contains(sup(i, 1), ["peak", "PEAK", "Peak"]) - outTypes = [outTypes, "peak_params"]; - elseif contains(sup(i, 1), ["gau", "GAU", "Gau"]) - outTypes = [outTypes, "gaussian_params"]; - elseif contains(sup(i, 1), ["err", "ERR", "Err"]) - outTypes = [outTypes, "error"]; - elseif contains(sup(i, 1), ["squar", "Squar", "SQUAR"]) - outTypes = [outTypes, "r_squared"]; - elseif contains(sup(i, 1), ["foo", "FOO", "Foo"]) - outTypes = [outTypes, "fooofed_spectrum"]; - elseif contains(sup(i, 1), ["bg", "BG", "Bg"]) - outTypes = [outTypes, "bg_fit"]; - elseif contains(sup(i, 1), ["pow", "POW", "Pow"]) - outTypes = [outTypes, "power_spectrum"]; - end - end - end - - for i = 1:length(outTypes) - if contains(outTypes(i), ["Exp","exp","EXP"]) - outTypes(i) = "exponent"; - elseif contains(outTypes(i), ["off","OFF", "Off"]) - outTypes(i) = "offset"; - elseif contains(outTypes(i), ["peak","PEAK", "Peak"]) - outTypes(i) = "peak_params"; - elseif contains(outTypes(i), ["gau","GAU", "Gau"]) - outTypes(i) = "gaussian_params"; - elseif contains(outTypes(i), ["err","ERR", "Err"]) - outTypes(i) = "error"; - elseif contains(outTypes(i), ["squar","Squar", "SQUAR"]) - outTypes(i) = "r_squared"; - elseif contains(outTypes(i), ["foo","FOO", "Foo"]) - outTypes(i) = "fooofed_spectrum"; - elseif contains(outTypes(i), ["bg","BG", "Bg", "fit", "FIT", ... - "Fit"]) - outTypes(i) = "bg_fit"; - elseif contains(outTypes(i), ["pow","POW", "Pow"]) - outTypes(i) = "power_spectrum"; - end - end - - % initial setting - settings = struct(); - settings.max_n_peaks = maxPeaks; - [time_series, fsOld] = load_data(strcat(inDir, cases(1).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - time_series = time_series(:, tStart:end); - data = squeeze(time_series(1, 1:dt)); - [~, w] = pwelch(data, [], 0, [], fs); - sup = [find(w > cf(1), 1), find(w > cf(1), 1)-1]; - [~, y] = min([w(find(w > cf(1), 1))-cf(1), ... - cf(1)-w(find(w > cf(1), 1)-1)]); - inferior = sup(y); - sup = [find(w > cf(end), 1), find(w > cf(end), 1)-1]; - [~,y] = min([w(find(w > cf(end), 1))-cf(end), ... - cf(end)-w(find(w > cf(end), 1)-1)]); - superior = sup(y); - f_range_ind = [inferior superior]; - f_range = [w(inferior) w(superior)]; - clear time_series - - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - time_series = time_series(:, tStart:end); - nLoc = size(time_series, 1); - - setup_data = zeros(nEpochs, nLoc); - offset.data = setup_data; - offset.locations = locations; - offset.chanlocs = chanlocs; - exponent = offset; - error = offset; - r_squared = offset; - setup_data = zeros(nEpochs, nLoc, maxPeaks*3); - peak_params.data = setup_data; - peak_params.locations = locations; - gaussian_params = peak_params; - setup_data = zeros(nEpochs, nLoc, f_range_ind(2)-f_range_ind(1)+1); - fooofed_spectrum.data = setup_data; - fooofed_spectrum.locations = locations; - bg_fit = fooofed_spectrum; - power_spectrum = bg_fit; - for k = 1:nEpochs - for j = 1:nLoc - data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); - [pxx, w] = pwelch(data, [], 0, [], fs); - fooof_results = fooof(w', pxx, f_range, settings, 1); - - offset.data(k, j) = fooof_results.background_params(1); - exponent.data(k, j) = ... - fooof_results.background_params(2); - - r_squared.data(k, j) = fooof_results.r_squared; - error.data(k, j) = fooof_results.error; - fooofed_spectrum.data(k, j, :) = ... - fooof_results.fooofed_spectrum; - bg_fit.data(k, j, :) = fooof_results.bg_fit; - power_spectrum.data(k, j, :) = ... - fooof_results.power_spectrum; - - sizePeaks = size(fooof_results.peak_params, 1)*3; - peak_params.data(k, j, 1:sizePeaks) = ... - squeeze(reshape(... - fooof_results.peak_params', 1, sizePeaks)); - gaussian_params.data(k, j, 1:sizePeaks) = ... - squeeze(reshape(... - fooof_results.gaussian_params', 1, sizePeaks)); - end - end - - for s = 1:length(outTypes) - outDir = path_check(subdir(inDir, outTypes(s))); - name = split(cases(i).name, '\'); - if length(name) == 1 - name = split(cases(i).name, '\'); - end - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - save(fullfile_check(filename), outTypes(s)); - end - clear time_series - catch - end - waitbar(i/length(cases), f) - end - close(f) +%% FOOOFer +% This function parameterizes the neural power spectra using the Haller et +% al. method. +% +% FOOOFer(fs, cf, nEpochs, dt, inDir, tStart, outTypes) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch +% outTypes is the list of variables to save (offset, exponent, +% peak parameters, gaussian parameters, error, r squared, fooofed +% spectrum, bg fit, power spectrum), writen in an array as strings +% in the same order as their output directories in the variable +% outDirs (es. outDirs=["C:\offset\", "C:\exponent", "C:\bg_fit"], +% outTypes=["offset", "exponent", "bg fit"]) +% +% NB: peak and gaussian params will contain a matrix where every group of 3 +% numbers are relative to one peak (CF, Amp, BW) and the others zeros +% are utilized to export only one matrix for each subject + + +function vargout = FOOOFer(fs, cf, nEpochs, dt, inDir, tStart, outTypes, maxPeaks) + switch nargin + case 5 + tStart = 0; + outTypes = ""; + maxPeaks = (cf(end)-cf(1))*2; + case 6 + outTypes = ""; + maxPeaks = (cf(end)-cf(1))*2; + case 7 + maxPeaks = (cf(end)-cf(1))*2; + end + + f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + + vargout = -1; + dt = fs*dt; + tStart = tStart*fs+1; + + askList = 'Insert the background measures separated by a comma'; + + inDir = path_check(inDir); + cases = define_cases(inDir); + outTypes = string(outTypes); + + if strcmp(outTypes, "") + outTypes = []; + sup = string(inputdlg(askList)); + if stcmp(sup, "exit") + return + end + sup = split(sup,","); + for i = 1:length(sup) + if contains(sup(i, 1), ["Exp", "exp", "EXP"]) + outTypes = [outTypes, "exp"]; + elseif contains(sup(i, 1), ["off", "OFF", "Off"]) + outTypes = [outTypes, "off"]; + elseif contains(sup(i, 1), ["peak", "PEAK", "Peak"]) + outTypes = [outTypes, "peak_params"]; + elseif contains(sup(i, 1), ["gau", "GAU", "Gau"]) + outTypes = [outTypes, "gaussian_params"]; + elseif contains(sup(i, 1), ["err", "ERR", "Err"]) + outTypes = [outTypes, "error"]; + elseif contains(sup(i, 1), ["squar", "Squar", "SQUAR"]) + outTypes = [outTypes, "r_squared"]; + elseif contains(sup(i, 1), ["foo", "FOO", "Foo"]) + outTypes = [outTypes, "fooofed_spectrum"]; + elseif contains(sup(i, 1), ["bg", "BG", "Bg"]) + outTypes = [outTypes, "bg_fit"]; + elseif contains(sup(i, 1), ["pow", "POW", "Pow"]) + outTypes = [outTypes, "power_spectrum"]; + end + end + end + + for i = 1:length(outTypes) + if contains(outTypes(i), ["Exp","exp","EXP"]) + outTypes(i) = "exponent"; + elseif contains(outTypes(i), ["off","OFF", "Off"]) + outTypes(i) = "offset"; + elseif contains(outTypes(i), ["peak","PEAK", "Peak"]) + outTypes(i) = "peak_params"; + elseif contains(outTypes(i), ["gau","GAU", "Gau"]) + outTypes(i) = "gaussian_params"; + elseif contains(outTypes(i), ["err","ERR", "Err"]) + outTypes(i) = "error"; + elseif contains(outTypes(i), ["squar","Squar", "SQUAR"]) + outTypes(i) = "r_squared"; + elseif contains(outTypes(i), ["foo","FOO", "Foo"]) + outTypes(i) = "fooofed_spectrum"; + elseif contains(outTypes(i), ["bg","BG", "Bg", "fit", "FIT", ... + "Fit"]) + outTypes(i) = "bg_fit"; + elseif contains(outTypes(i), ["pow","POW", "Pow"]) + outTypes(i) = "power_spectrum"; + end + end + + % initial setting + settings = struct(); + settings.max_n_peaks = maxPeaks; + [time_series, fsOld] = load_data(strcat(inDir, cases(1).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + time_series = time_series(:, tStart:end); + data = squeeze(time_series(1, 1:dt)); + [~, w] = pwelch(data, [], 0, [], fs); + sup = [find(w > cf(1), 1), find(w > cf(1), 1)-1]; + [~, y] = min([w(find(w > cf(1), 1))-cf(1), ... + cf(1)-w(find(w > cf(1), 1)-1)]); + inferior = sup(y); + sup = [find(w > cf(end), 1), find(w > cf(end), 1)-1]; + [~,y] = min([w(find(w > cf(end), 1))-cf(end), ... + cf(end)-w(find(w > cf(end), 1)-1)]); + superior = sup(y); + f_range_ind = [inferior superior]; + f_range = [w(inferior) w(superior)]; + clear time_series + + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + time_series = time_series(:, tStart:end); + nLoc = size(time_series, 1); + + setup_data = zeros(nEpochs, nLoc); + offset.data = setup_data; + offset.locations = locations; + offset.chanlocs = chanlocs; + exponent = offset; + error = offset; + r_squared = offset; + setup_data = zeros(nEpochs, nLoc, maxPeaks*3); + peak_params.data = setup_data; + peak_params.locations = locations; + gaussian_params = peak_params; + setup_data = zeros(nEpochs, nLoc, f_range_ind(2)-f_range_ind(1)+1); + fooofed_spectrum.data = setup_data; + fooofed_spectrum.locations = locations; + bg_fit = fooofed_spectrum; + power_spectrum = bg_fit; + for k = 1:nEpochs + for j = 1:nLoc + data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); + [pxx, w] = pwelch(data, [], 0, [], fs); + fooof_results = fooof(w', pxx, f_range, settings, 1); + + offset.data(k, j) = fooof_results.background_params(1); + exponent.data(k, j) = ... + fooof_results.background_params(2); + + r_squared.data(k, j) = fooof_results.r_squared; + error.data(k, j) = fooof_results.error; + fooofed_spectrum.data(k, j, :) = ... + fooof_results.fooofed_spectrum; + bg_fit.data(k, j, :) = fooof_results.bg_fit; + power_spectrum.data(k, j, :) = ... + fooof_results.power_spectrum; + + sizePeaks = size(fooof_results.peak_params, 1)*3; + peak_params.data(k, j, 1:sizePeaks) = ... + squeeze(reshape(... + fooof_results.peak_params', 1, sizePeaks)); + gaussian_params.data(k, j, 1:sizePeaks) = ... + squeeze(reshape(... + fooof_results.gaussian_params', 1, sizePeaks)); + vargout = 0; + end + end + + for s = 1:length(outTypes) + outDir = path_check(subdir(inDir, outTypes(s))); + name = split(cases(i).name, '\'); + if length(name) == 1 + name = split(cases(i).name, '\'); + end + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + save(fullfile_check(filename), outTypes(s)); + end + clear time_series + catch + end + waitbar(i/length(cases), f) + end + close(f) end \ No newline at end of file diff --git a/Measures/PSDr.m b/Measures/PSDr.m index 225070a..e1c697a 100644 --- a/Measures/PSDr.m +++ b/Measures/PSDr.m @@ -1,97 +1,102 @@ -%% PSDr -% This function computes the Power Spectral Density relative over a total -% band defined by the user, using the Welch's power formulation. -% -% PSDr(fs, cf, nEpochs, dt, inDir, tStart, relBand) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 4 8 13 30 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch (0 as default) -% relBand is the total band used to obtain the relative band -% ([min(min(cf),1) max(max(cf),40)] as default, the first value is -% the minimum between 1 and the first cut frequency and the second -% value is the maximum between 40 and the last cut frequency) - - -function PSDr(fs, cf, nEpochs, dt, inDir, tStart, relBand) - switch nargin - case 5 - tStart = 0; - relBand = [min(min(cf), 1) max(max(cf), 40)]; - case 6 - relBand = [min(min(cf), 1) max(max(cf), 40)]; - end - - f = waitbar(0,'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - - cfstart = 0; - cfstop = length(cf)-1; - nBands = cfstop-cfstart; - dt = fs*dt; - tStart = tStart*fs+1; - inDir = path_check(inDir); - if iscell(relBand) - aux_relBand = [str2double(string(relBand{1})), ... - str2double(string(relBand{2}))]; - relBand = aux_relBand; - end - - cases = define_cases(inDir); - - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - time_series = time_series(:, tStart:end); - psdr.data = zeros(nBands, nEpochs, size(time_series, 1)); - psdr.locations = locations; - psdr.chanlocs = chanlocs; - for k = 1:nEpochs - - for j = 1:size(time_series, 1) - data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); - [pxx, w] = pwelch(data, [], [], [], fs); - bandPower = zeros(nBands, 1); - - for b = 1:nBands - [infft, supft] = band_boundaries(w, ... - cf(b+cfstart), cf(b+cfstart+1)); - - bandPower(b, 1) = sum(pxx(infft:supft)); - end - - [infft, supft] = band_boundaries(w, relBand(1), ... - relBand(end)); - totalPower = sum(pxx(infft:supft)); - - for b = 1:nBands - psdr.data(b, k, j) = bandPower(b, 1)/totalPower; - end - end - end - - outDir = path_check(subdir(inDir, 'PSDr')); - name = split(cases(i).name, filesep); - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - - save(fullfile_check(filename), 'psdr'); - catch - end %end try - waitbar(i/length(cases), f) - end - close(f) +%% PSDr +% This function computes the Power Spectral Density relative over a total +% band defined by the user, using the Welch's power formulation. +% +% PSDr(fs, cf, nEpochs, dt, inDir, tStart, relBand) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 4 8 13 30 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch (0 as default) +% relBand is the total band used to obtain the relative band +% ([min(min(cf),1) max(max(cf),40)] as default, the first value is +% the minimum between 1 and the first cut frequency and the second +% value is the maximum between 40 and the last cut frequency) + + +function vargout = PSDr(fs, cf, nEpochs, dt, inDir, tStart, relBand) + switch nargin + case 5 + tStart = 0; + relBand = [min(min(cf), 1) max(max(cf), 40)]; + case 6 + relBand = [min(min(cf), 1) max(max(cf), 40)]; + end + + f = waitbar(0,'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + + vargout = -1; + cfstart = 0; + cfstop = length(cf)-1; + nBands = cfstop-cfstart; + dt = fs*dt; + tStart = tStart*fs+1; + inDir = path_check(inDir); + if iscell(relBand) + aux_relBand = [str2double(string(relBand{1})), ... + str2double(string(relBand{2}))]; + relBand = aux_relBand; + end + + cases = define_cases(inDir); + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + time_series = time_series(:, tStart:end); + psdr.data = zeros(nBands, nEpochs, size(time_series, 1)); + psdr.locations = locations; + psdr.chanlocs = chanlocs; + + for k = 1:nEpochs + for j = 1:size(time_series, 1) + + data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); + [pxx, w] = pwelch(data, [], [], [], fs); + bandPower = zeros(nBands, 1); + + for b = 1:nEpochs + [infft, supft] = band_boundaries(w, ... + cf(b+cfstart), cf(b+cfstart+1)); + + bandPower(b, 1) = sum(pxx(infft:supft)); + vargout = 0; + end + + [infft, supft] = band_boundaries(w, relBand(1), ... + relBand(end)); + + totalPower = sum(pxx(infft:supft)); + + for b = 1:nBands + psdr.data(b, k, j) = bandPower(b, 1)/totalPower; + end + + end + end + + outDir = path_check(subdir(inDir, 'PSDr')); + name = split(cases(i).name, filesep); + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + + save(fullfile_check(filename), 'psdr'); + catch + end + %end try + waitbar(i/length(cases), f) + end + close(f) end \ No newline at end of file diff --git a/Measures/autocorrelation_measures.m b/Measures/autocorrelation_measures.m index 7167403..13eaeff 100644 --- a/Measures/autocorrelation_measures.m +++ b/Measures/autocorrelation_measures.m @@ -1,142 +1,144 @@ -%% autocorrelation_measures -% This function computes the autocorrelation evaluation related to EEG -% channels or ROIs. -% -% autocorrelation_measures(fs, cf, nEpochs, dt, inDir, outDirs, tStart, ... -% outTypes, filter_name) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch -% outTypes is the list of autocorrelation measures to compute (hurst -% exponent) -% filter_name is the name of the filtering function (athena_filter as -% default) - - -function autocorrelation_measures(fs, cf, nEpochs, dt, inDir, tStart, ... - outTypes, filter_name) - - if nargin < 6 - tStart = 0; - end - if nargin < 7 - outTypes = ""; - end - if nargin < 8 - filter_name = 'athena_filter'; - end - if nargin < 9 - m = 2; - end - if nargin < 10 - r = 0.2; - end - - f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - askList = 'Insert the autocorrelation measures separated by a comma'; - - filter_handle = eval(strcat('@', filter_name)); - dt = fs*dt; - tStart = tStart*fs+1; - nBands = length(cf)-1; - - hurst_names = ["hurst_exponent", "hurst"]; - - inDir = path_check(inDir); - cases = define_cases(inDir); - outTypes = string(outTypes); - - ctrl = 0; - if length(outTypes) == 1 - if strcmp(outTypes,'') - ctrl = 1; - end - end - - while ctrl == 1 - outTypes = []; - sup = string(inputdlg(askList)); - if strcmp(sup, 'exit') - return - end - sup = split(sup, ','); - for i = 1:length(sup) - if contains(sup(i, 1), hurst_names) - outTypes = [outTypes, "Hurst"]; - end - end - if length(outDirs) == length(outTypes) - ctrl = 0; - end - end - - for i = 1:length(outTypes) - if contains(outTypes(i), hurst_names) - outTypes(i) = "Hurst"; - end - end - - [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - check = check_filtering(time_series, dt, tStart, fs, cf, ... - filter_handle); - if check - close(f) - return; - end - - for c = 1:length(outTypes) - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - nLoc = size(time_series, 1); - he.data = zeros(nBands, nEpochs, nLoc); - he.locations = locations; - he.chanlocs = chanlocs; - for j = 1:nBands - for k = 1:nEpochs - for loc = 1:nLoc - ti = ((k-1)*dt)+tStart; - tf = k*dt+tStart-1; - data = filter_handle(... - time_series(loc, ti:tf), fs, cf(j), ... - cf(j+1)); - data = data'; - if strcmpi(outTypes(c), "Hurst") - he.data(j, k, loc) = ... - hurst_exponent(data); - end - end - end - end - outDir = path_check(subdir(inDir, outTypes(c))); - name = split(cases(i).name, '\'); - if length(name) == 1 - name = split(cases(i).name, '\'); - end - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - save(filename, 'he'); - catch - end %end try - waitbar((i+(c-1)*length(cases))/... - (length(cases)*length(outTypes)), f) - end - end - close(f) +%% autocorrelation_measures +% This function computes the autocorrelation evaluation related to EEG +% channels or ROIs. +% +% autocorrelation_measures(fs, cf, nEpochs, dt, inDir, outDirs, tStart, ... +% outTypes, filter_name) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch +% outTypes is the list of autocorrelation measures to compute (hurst +% exponent) +% filter_name is the name of the filtering function (athena_filter as +% default) + + +function vargout = autocorrelation_measures(fs, cf, nEpochs, dt, inDir, tStart, ... + outTypes, filter_name) + + if nargin < 6 + tStart = 0; + end + if nargin < 7 + outTypes = ""; + end + if nargin < 8 + filter_name = 'athena_filter'; + end + if nargin < 9 + m = 2; + end + if nargin < 10 + r = 0.2; + end + + vargout = -1; + f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + askList = 'Insert the autocorrelation measures separated by a comma'; + + filter_handle = eval(strcat('@', filter_name)); + dt = fs*dt; + tStart = tStart*fs+1; + nBands = length(cf)-1; + + hurst_names = ["hurst_exponent", "hurst"]; + + inDir = path_check(inDir); + cases = define_cases(inDir); + outTypes = string(outTypes); + + ctrl = 0; + if length(outTypes) == 1 + if strcmp(outTypes,'') + ctrl = 1; + end + end + + while ctrl == 1 + outTypes = []; + sup = string(inputdlg(askList)); + if strcmp(sup, 'exit') + return + end + sup = split(sup, ','); + for i = 1:length(sup) + if contains(sup(i, 1), hurst_names) + outTypes = [outTypes, "Hurst"]; + end + end + if length(outDirs) == length(outTypes) + ctrl = 0; + end + end + + for i = 1:length(outTypes) + if contains(outTypes(i), hurst_names) + outTypes(i) = "Hurst"; + end + end + + [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + check = check_filtering(time_series, dt, tStart, fs, cf, ... + filter_handle); + if check + close(f) + return; + end + + for c = 1:length(outTypes) + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + nLoc = size(time_series, 1); + he.data = zeros(nBands, nEpochs, nLoc); + he.locations = locations; + he.chanlocs = chanlocs; + for j = 1:nBands + for k = 1:nEpochs + for loc = 1:nLoc + ti = ((k-1)*dt)+tStart; + tf = k*dt+tStart-1; + data = filter_handle(... + time_series(loc, ti:tf), fs, cf(j), ... + cf(j+1)); + data = data'; + if strcmpi(outTypes(c), "Hurst") + he.data(j, k, loc) = ... + hurst_exponent(data); + vargout = 0; + end + end + end + end + outDir = path_check(subdir(inDir, outTypes(c))); + name = split(cases(i).name, '\'); + if length(name) == 1 + name = split(cases(i).name, '\'); + end + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + save(filename, 'he'); + catch + end %end try + waitbar((i+(c-1)*length(cases))/... + (length(cases)*length(outTypes)), f) + end + end + close(f) end \ No newline at end of file diff --git a/Measures/connectivity.m b/Measures/connectivity.m index c80d80c..4d90594 100644 --- a/Measures/connectivity.m +++ b/Measures/connectivity.m @@ -1,201 +1,211 @@ -%% connectivity -% This function computes the connectivity between EEG channels or ROIs, -% using the phase locking value (PLV), the phase lag index (PLI), the -% magnitude-squared coherence (MSC), the amplitude envelope correlation -% corrected (also called orthogonalized, AECo), the amplitude envelope -% correlation not corrected (AEC), the mutual information (MI), or the -% weighted phase lag index (wPLI). -% -% connectivity(fs, cf, nEpochs, dt, inDir, tStart, outTypes, filter_name) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch -% outTypes is the list of connectivity types to compute (PLI, PLV, etc.) -% filter_name is the name of the filtering function (athena_filter as -% default) - -function connectivity(fs, cf, nEpochs, dt, inDir, tStart, outTypes, ... - filter_name) - - switch nargin - case 5 - tStart=0; - outTypes=""; - filter_name = 'athena_filter'; - case 6 - outTypes=""; - filter_name = 'athena_filter'; - case 7 - filter_name = 'athena_filter'; - end - - f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - askList = 'Insert the connectivity measures separated by a comma'; - - filter_handle = eval(strcat('@', filter_name)); - dt = fs*dt; - tStart = tStart*fs+1; - nBands = length(cf)-1; - - PLInames = ["pli", "PLI", "Pli"]; - PLVnames = ["plv", "PLV", "Plv"]; - AECnames = ["aec", "AEC", "Aec"]; - AECOnames = ["aeco", "aec_o", "AECo", "AEC_o", "Aeco", "Aec_o", ... - "AECc", "AEC_c", "aecc", "aec_c", "Aecc", "Aec_c"]; - MSCnames = ["coherence", "MSC", "coh", "msc", "COH", "Coherence"]; - ICOHnames = ["ICOH", "Coherency", "coherency"]; - MInames = ["mutual_information", "MI", "Mutual_Information"]; - CCnames = ["correlation", "correlation_coefficient", "CC"]; - wPLInames = ["wPLI", "w_PLI"]; - - inDir = path_check(inDir); - cases = define_cases(inDir); - outTypes = string(outTypes); - - ctrl = 0; - if length(outTypes) == 1 - if strcmp(outTypes,'') - ctrl = 1; - end - end - - while ctrl == 1 - outTypes = []; - sup = string(inputdlg(askList)); - if strcmp(sup, 'exit') - return - end - sup = split(sup, ','); - for i = 1:length(sup) - if contains(sup(i, 1), wPLInames) - outTypes = [outTypes, "wPLI"]; - elseif contains(sup(i, 1), PLVnames) - outTypes = [outTypes, "PLV"]; - elseif contains(sup(i, 1), AECOnames) - outTypes = [outTypes, "AECo"]; - elseif contains(sup(i, 1), AECnames) - outTypes = [outTypes, "AEC"]; - elseif contains(sup(i, 1), MSCnames) - outTypes = [outTypes, "coherence"]; - elseif contains(sup(i, 1), ICOHnames) - outTypes = [outTypes, "ICOH"]; - elseif contains(sup(i, 1), MInames) - outTypes = [outTypes, "mutual_information"]; - elseif contains(sup(i, 1), CCnames) - outTypes = [outTypes, "correlation_coefficient"]; - elseif contains(sup(i, 1), PLInames) - outTypes = [outTypes, "PLI"]; - end - end - if length(outDirs) == length(outTypes) - ctrl = 0; - end - end - - for i = 1:length(outTypes) - if contains(outTypes(i), wPLInames) - outTypes(i) = "wPLI"; - elseif contains(outTypes(i), PLVnames) - outTypes(i) = "PLV"; - elseif contains(outTypes(i), AECOnames) - outTypes(i) = "AECo"; - elseif contains(outTypes(i), AECnames) - outTypes(i) = "AEC"; - elseif contains(outTypes(i), ICOHnames) - outTypes(i) = "ICOH"; - elseif contains(outTypes(i), MSCnames) - outTypes(i) = "coherence"; - elseif contains(outTypes(i), MInames) - outTypes(i) = "mutual_information"; - elseif contains(outTypes(i), CCnames) - outTypes(i) = "correlation_coefficient"; - elseif contains(outTypes(i), PLInames) - outTypes(i) = "PLI"; - end - end - - [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - check = check_filtering(time_series, dt, tStart, fs, cf, ... - filter_handle); - if check - close(f) - return; - end - - for c = 1:length(outTypes) - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - nLoc = size(time_series, 1); - conn.data = zeros(nBands, nEpochs, nLoc, nLoc); - conn.locations = locations; - conn.chanlocs = chanlocs; - for j = 1:nBands - for k = 1:nEpochs - ti = ((k-1)*dt)+tStart; - tf = k*dt+tStart-1; - data = filter_handle(time_series(:, ti:tf), fs, ... - cf(j), cf(j+1)); - data = data'; - if strcmpi(outTypes(c), "PLI") - conn.data(j, k, :, :) = phase_lag_index(data); - elseif strcmpi(outTypes(c), "PLV") - conn.data(j, k, :, :) = ... - phase_locking_value(data); - elseif strcmpi(outTypes(c), "AECo") - conn.data(j, k, :, :) = ... - amplitude_envelope_correlation_orth(data); - elseif strcmpi(outTypes(c), "AEC") - conn.data(j, k, :, :) = ... - amplitude_envelope_correlation(data); - elseif strcmpi(outTypes(c), "coherence") - conn.data(j, k, :, :) = ... - magnitude_squared_coherence(data); - elseif strcmpi(outTypes(c), "ICOH") - conn.data(j, k, :, :) = ... - imaginary_coherency(data); - elseif strcmpi(outTypes(c), "mutual_information") - conn.data(j, k, :, :) = ... - mutual_information(data, 'max'); - elseif strcmpi(outTypes(c), ... - "correlation_coefficient") - conn.data(j, k, :, :) = ... - correlation_coefficient(data); - elseif strcmpi(outTypes(c), "wPLI") - conn.data(j, k, :, :) = ... - weighted_phase_lag_index(data); - end - end - end - outDir = path_check(subdir(inDir, outTypes(c))); - name = split(cases(i).name, '\'); - if length(name) == 1 - name = split(cases(i).name, '\'); - end - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - save(fullfile_check(filename), 'conn'); - catch - end %end try - waitbar((i+(c-1)*length(cases))/... - (length(cases)*length(outTypes)), f) - end - end - close(f) +%% connectivity +% This function computes the connectivity between EEG channels or ROIs, +% using the phase locking value (PLV), the phase lag index (PLI), the +% magnitude-squared coherence (MSC), the amplitude envelope correlation +% corrected (also called orthogonalized, AECo), the amplitude envelope +% correlation not corrected (AEC), the mutual information (MI), or the +% weighted phase lag index (wPLI). +% +% connectivity(fs, cf, nEpochs, dt, inDir, tStart, outTypes, filter_name) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch +% outTypes is the list of connectivity types to compute (PLI, PLV, etc.) +% filter_name is the name of the filtering function (athena_filter as +% default) + +function vargout = connectivity(fs, cf, nEpochs, dt, inDir, tStart, outTypes, ... + filter_name) + + switch nargin + case 5 + tStart=0; + outTypes=""; + filter_name = 'athena_filter'; + case 6 + outTypes=""; + filter_name = 'athena_filter'; + case 7 + filter_name = 'athena_filter'; + end + + vargout = -1; + f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + askList = 'Insert the connectivity measures separated by a comma'; + + filter_handle = eval(strcat('@', filter_name)); + dt = fs*dt; + tStart = tStart*fs+1; + nBands = length(cf)-1; + + PLInames = ["pli", "PLI", "Pli"]; + PLVnames = ["plv", "PLV", "Plv"]; + AECnames = ["aec", "AEC", "Aec"]; + AECOnames = ["aeco", "aec_o", "AECo", "AEC_o", "Aeco", "Aec_o", ... + "AECc", "AEC_c", "aecc", "aec_c", "Aecc", "Aec_c"]; + MSCnames = ["coherence", "MSC", "coh", "msc", "COH", "Coherence"]; + ICOHnames = ["ICOH", "Coherency", "coherency"]; + MInames = ["mutual_information", "MI", "Mutual_Information"]; + CCnames = ["correlation", "correlation_coefficient", "CC"]; + wPLInames = ["wPLI", "w_PLI"]; + + inDir = path_check(inDir); + cases = define_cases(inDir); + outTypes = string(outTypes); + + ctrl = 0; + if length(outTypes) == 1 + if strcmp(outTypes,'') + ctrl = 1; + end + end + + while ctrl == 1 + outTypes = []; + sup = string(inputdlg(askList)); + if strcmp(sup, 'exit') + return + end + sup = split(sup, ','); + for i = 1:length(sup) + if contains(sup(i, 1), wPLInames) + outTypes = [outTypes, "wPLI"]; + elseif contains(sup(i, 1), PLVnames) + outTypes = [outTypes, "PLV"]; + elseif contains(sup(i, 1), AECOnames) + outTypes = [outTypes, "AECo"]; + elseif contains(sup(i, 1), AECnames) + outTypes = [outTypes, "AEC"]; + elseif contains(sup(i, 1), MSCnames) + outTypes = [outTypes, "coherence"]; + elseif contains(sup(i, 1), ICOHnames) + outTypes = [outTypes, "ICOH"]; + elseif contains(sup(i, 1), MInames) + outTypes = [outTypes, "mutual_information"]; + elseif contains(sup(i, 1), CCnames) + outTypes = [outTypes, "correlation_coefficient"]; + elseif contains(sup(i, 1), PLInames) + outTypes = [outTypes, "PLI"]; + end + end + if length(outDirs) == length(outTypes) + ctrl = 0; + end + end + + for i = 1:length(outTypes) + if contains(outTypes(i), wPLInames) + outTypes(i) = "wPLI"; + elseif contains(outTypes(i), PLVnames) + outTypes(i) = "PLV"; + elseif contains(outTypes(i), AECOnames) + outTypes(i) = "AECo"; + elseif contains(outTypes(i), AECnames) + outTypes(i) = "AEC"; + elseif contains(outTypes(i), ICOHnames) + outTypes(i) = "ICOH"; + elseif contains(outTypes(i), MSCnames) + outTypes(i) = "coherence"; + elseif contains(outTypes(i), MInames) + outTypes(i) = "mutual_information"; + elseif contains(outTypes(i), CCnames) + outTypes(i) = "correlation_coefficient"; + elseif contains(outTypes(i), PLInames) + outTypes(i) = "PLI"; + end + end + + [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + check = check_filtering(time_series, dt, tStart, fs, cf, ... + filter_handle); + if check + close(f) + return; + end + + for c = 1:length(outTypes) + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + nLoc = size(time_series, 1); + conn.data = zeros(nBands, nEpochs, nLoc, nLoc); + conn.locations = locations; + conn.chanlocs = chanlocs; + for j = 1:nBands + for k = 1:nEpochs + ti = ((k-1)*dt)+tStart; + tf = k*dt+tStart-1; + data = filter_handle(time_series(:, ti:tf), fs, ... + cf(j), cf(j+1)); + data = data'; + if strcmpi(outTypes(c), "PLI") + conn.data(j, k, :, :) = phase_lag_index(data); + vargout = 0; + elseif strcmpi(outTypes(c), "PLV") + conn.data(j, k, :, :) = ... + phase_locking_value(data); + vargout = 0; + elseif strcmpi(outTypes(c), "AECo") + conn.data(j, k, :, :) = ... + amplitude_envelope_correlation_orth(data); + vargout = 0; + elseif strcmpi(outTypes(c), "AEC") + conn.data(j, k, :, :) = ... + amplitude_envelope_correlation(data); + vargout = 0; + elseif strcmpi(outTypes(c), "coherence") + conn.data(j, k, :, :) = ... + magnitude_squared_coherence(data); + vargout = 0; + elseif strcmpi(outTypes(c), "ICOH") + conn.data(j, k, :, :) = ... + imaginary_coherency(data); + vargout = 0; + elseif strcmpi(outTypes(c), "mutual_information") + conn.data(j, k, :, :) = ... + mutual_information(data, 'max'); + vargout = 0; + elseif strcmpi(outTypes(c), ... + "correlation_coefficient") + conn.data(j, k, :, :) = ... + correlation_coefficient(data); + vargout = 0; + elseif strcmpi(outTypes(c), "wPLI") + conn.data(j, k, :, :) = ... + weighted_phase_lag_index(data); + vargout = 0; + end + end + end + outDir = path_check(subdir(inDir, outTypes(c))); + name = split(cases(i).name, '\'); + if length(name) == 1 + name = split(cases(i).name, '\'); + end + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + save(fullfile_check(filename), 'conn'); + catch + end %end try + waitbar((i+(c-1)*length(cases))/... + (length(cases)*length(outTypes)), f) + end + end + close(f) end \ No newline at end of file diff --git a/Measures/spectral_entropy.m b/Measures/spectral_entropy.m index f79b498..63fafbe 100644 --- a/Measures/spectral_entropy.m +++ b/Measures/spectral_entropy.m @@ -1,74 +1,75 @@ -%% spectral_entropy -% This function computes the Spectral Entropy, on the single frequency -% bands and on the single epochs. -% -% spectral_entropy(fs, cf, nEpochs, dt, inDir, outDir, tStart) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 4 8 13 30 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch (0 as default) - - -function spectral_entropy(fs, cf, nEpochs, dt, inDir, tStart) - switch nargin - case 5 - tStart = 0; - end - - f = waitbar(0,'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - - cfstart = 0; - cfstop = length(cf)-1; - nBands = cfstop-cfstart; - dt = fs*dt; - tStart = tStart*fs+1; - inDir = path_check(inDir); - - cases = define_cases(inDir); - - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - time_series = time_series(:, tStart:end); - se.data = zeros(nBands, nEpochs, size(time_series, 1)); - se.locations = locations; - se.chanlocs = chanlocs; - for k = 1:nEpochs - - for j = 1:size(time_series, 1) - data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); - for b = 1:nBands - [p, fp, tp] = pspectrum(data, fs, 'spectrogram'); - se.data(b, k, j) = pentropy(p, fp, tp, ... - 'Instantaneous', false, ... - 'FrequencyLimits', [cf(b), cf(b+1)]); - end - end - end - - outDir = path_check(subdir(inDir, 'PEntropy')); - name = split(cases(i).name, filesep); - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - - save(fullfile_check(filename), 'se'); - catch - end - waitbar(i/length(cases), f) - end - close(f) +%% +% This function computes the Spectral Entropy, on the single frequency +% bands and on the single epochs. +% +% spectral_entropy(fs, cf, nEpochs, dt, inDir, outDir, tStart) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 4 8 13 30 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch (0 as default) + + +function vargout = spectral_entropy(fs, cf, nEpochs, dt, inDir, tStart) + switch nargin + case 5 + tStart = 0; + end + + f = waitbar(0,'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + + vargout = -1; + cfstart = 0; + cfstop = length(cf)-1; + nBands = cfstop-cfstart; + dt = fs*dt; + tStart = tStart*fs+1; + inDir = path_check(inDir); + cases = define_cases(inDir); + + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + time_series = time_series(:, tStart:end); + se.data = zeros(nBands, nEpochs, size(time_series, 1)); + se.locations = locations; + se.chanlocs = chanlocs; + for k = 1:nEpochs + + for j = 1:size(time_series, 1) + data = squeeze(time_series(j, dt*(k-1)+1:k*dt)); + for b = 1:nBands + [p, fp, tp] = pspectrum(data, fs, ' '); + se.data(b, k, j) = pentropy(p, fp, tp, ... + 'Instantaneous', false, ... + 'FrequencyLimits', [cf(b), cf(b+1)]); + vargout = 0; + end + end + end + + outDir = path_check(subdir(inDir, 'PEntropy')); + name = split(cases(i).name, filesep); + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + + save(fullfile_check(filename), 'se'); + catch + end + waitbar(i/length(cases), f) + end + close(f) end \ No newline at end of file diff --git a/Measures/statistical_information.m b/Measures/statistical_information.m index 372b122..b8a51d8 100644 --- a/Measures/statistical_information.m +++ b/Measures/statistical_information.m @@ -1,176 +1,183 @@ -%% statistical_information -% This function computes statistical measures on EEG channels or ROIs. -% -% statistical_information(fs, cf, nEpochs, dt, inDir, outDirs, tStart, ... -% outTypes, filter_name) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch -% outTypes is the list of statistical measures to compute (median, mean, -% std, variance, skewness, kurtosis) -% filter_name is the name of the filtering function (athena_filter as -% default) - - -function statistical_information(fs, cf, nEpochs, dt, inDir, tStart, ... - outTypes, filter_name) - - if nargin < 6 - tStart = 0; - end - if nargin < 7 - outTypes = ""; - end - if nargin < 8 - filter_name = 'athena_filter'; - end - if nargin < 9 - m = 2; - end - if nargin < 10 - r = 0.2; - end - - f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - askList = 'Insert the statistical measures separated by a comma'; - - filter_handle = eval(strcat('@', filter_name)); - dt = fs*dt; - tStart = tStart*fs+1; - nBands = length(cf)-1; - - variance_names = ["variance", "var"]; - mean_names = ["mean"]; - median_names = ["median"]; - std_names = ["std", "standard_deviation"]; - kurtosis_names = ["kurtosis"]; - skewness_names = ["skewness"]; - - inDir = path_check(inDir); - cases = define_cases(inDir); - outTypes = string(outTypes); - - ctrl = 0; - if length(outTypes) == 1 - if strcmp(outTypes,'') - ctrl = 1; - end - end - - while ctrl == 1 - outTypes = []; - sup = string(inputdlg(askList)); - if strcmp(sup, 'exit') - return - end - sup = split(sup, ','); - for i = 1:length(sup) - if contains(sup(i, 1), variance_names) - outTypes = [outTypes, "variance"]; - elseif contains(sup(i, 1), median_names) - outTypes = [outTypes, "median"]; - elseif contains(sup(i, 1), mean_names) - outTypes = [outTypes, "mean"]; - elseif contains(sup(i, 1), kurtosis_names) - outTypes = [outTypes, "kurtosis"]; - elseif contains(sup(i, 1), skewness_names) - outTypes = [outTypes, "skewness"]; - elseif contains(sup(i, 1), std_names) - outTypes = [outTypes, "Standard_deviation"]; - end - end - if length(outDirs) == length(outTypes) - ctrl = 0; - end - end - - for i = 1:length(outTypes) - if contains(outTypes(i), variance_names) - outTypes(i) = "Variance"; - elseif contains(outTypes(i), median_names) - outTypes(i) = "Median"; - elseif contains(outTypes(i), mean_names) - outTypes(i) = "Mean"; - elseif contains(outTypes(i), kurtosis_names) - outTypes(i) = "Kurtosis"; - elseif contains(outTypes(i), skewness_names) - outTypes(i) = "Skewness"; - elseif contains(outTypes(i), std_names) - outTypes(i) = "Standard_deviation"; - end - end - - [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - check = check_filtering(time_series, dt, tStart, fs, cf, ... - filter_handle); - if check - close(f) - return; - end - - for c = 1:length(outTypes) - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - nLoc = size(time_series, 1); - sa.data = zeros(nBands, nEpochs, nLoc); - sa.locations = locations; - sa.chanlocs = chanlocs; - for j = 1:nBands - for k = 1:nEpochs - for loc = 1:nLoc - ti = ((k-1)*dt)+tStart; - tf = k*dt+tStart-1; - data = filter_handle(... - time_series(loc, ti:tf), fs, cf(j), ... - cf(j+1)); - data = data'; - if strcmpi(outTypes(c), "Variance") - sa.data(j, k, loc) = var(data); - elseif strcmpi(outTypes(c), ... - "Standard_deviation") - sa.data(j, k, loc) = std(data); - elseif strcmpi(outTypes(c), "Mean") - sa.data(j, k, loc) = mean(data); - elseif strcmpi(outTypes(c), "Median") - sa.data(j, k, loc) = median(data); - elseif strcmpi(outTypes(c), "Kurtosis") - sa.data(j, k, loc) = kurtosis(data); - elseif strcmpi(outTypes(c), "Skewness") - sa.data(j, k, loc) = skewness(data); - end - end - end - end - outDir = path_check(subdir(inDir, outTypes(c))); - name = split(cases(i).name, '\'); - if length(name) == 1 - name = split(cases(i).name, '\'); - end - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - save(fullfile_check(filename), 'sa'); - catch - end %end try - waitbar((i+(c-1)*length(cases))/... - (length(cases)*length(outTypes)), f) - end - end - close(f) +%% statistical_information +% This function computes statistical measures on EEG channels or ROIs. +% +% statistical_information(fs, cf, nEpochs, dt, inDir, outDirs, tStart, ... +% outTypes, filter_name) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch +% outTypes is the list of statistical measures to compute (median, mean, +% std, variance, skewness, kurtosis) +% filter_name is the name of the filtering function (athena_filter as +% default) + + +function vargout = statistical_information(fs, cf, nEpochs, dt, inDir, tStart, ... + outTypes, filter_name) + + if nargin < 6 + tStart = 0; + end + if nargin < 7 + outTypes = ""; + end + if nargin < 8 + filter_name = 'athena_filter'; + end + if nargin < 9 + m = 2; + end + if nargin < 10 + r = 0.2; + end + + f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + askList = 'Insert the statistical measures separated by a comma'; + + vargout = -1; + filter_handle = eval(strcat('@', filter_name)); + dt = fs*dt; + tStart = tStart*fs+1; + nBands = length(cf)-1; + + variance_names = ["variance", "var"]; + mean_names = ["mean"]; + median_names = ["median"]; + std_names = ["std", "standard_deviation"]; + kurtosis_names = ["kurtosis"]; + skewness_names = ["skewness"]; + + inDir = path_check(inDir); + cases = define_cases(inDir); + outTypes = string(outTypes); + + ctrl = 0; + if length(outTypes) == 1 + if strcmp(outTypes,'') + ctrl = 1; + end + end + + while ctrl == 1 + outTypes = []; + sup = string(inputdlg(askList)); + if strcmp(sup, 'exit') + return + end + sup = split(sup, ','); + for i = 1:length(sup) + if contains(sup(i, 1), variance_names) + outTypes = [outTypes, "variance"]; + elseif contains(sup(i, 1), median_names) + outTypes = [outTypes, "median"]; + elseif contains(sup(i, 1), mean_names) + outTypes = [outTypes, "mean"]; + elseif contains(sup(i, 1), kurtosis_names) + outTypes = [outTypes, "kurtosis"]; + elseif contains(sup(i, 1), skewness_names) + outTypes = [outTypes, "skewness"]; + elseif contains(sup(i, 1), std_names) + outTypes = [outTypes, "Standard_deviation"]; + end + end + if length(outDirs) == length(outTypes) + ctrl = 0; + end + end + + for i = 1:length(outTypes) + if contains(outTypes(i), variance_names) + outTypes(i) = "Variance"; + elseif contains(outTypes(i), median_names) + outTypes(i) = "Median"; + elseif contains(outTypes(i), mean_names) + outTypes(i) = "Mean"; + elseif contains(outTypes(i), kurtosis_names) + outTypes(i) = "Kurtosis"; + elseif contains(outTypes(i), skewness_names) + outTypes(i) = "Skewness"; + elseif contains(outTypes(i), std_names) + outTypes(i) = "Standard_deviation"; + end + end + + [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + check = check_filtering(time_series, dt, tStart, fs, cf, ... + filter_handle); + if check + close(f) + return; + end + + for c = 1:length(outTypes) + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + nLoc = size(time_series, 1); + sa.data = zeros(nBands, nEpochs, nLoc); + sa.locations = locations; + sa.chanlocs = chanlocs; + for j = 1:nBands + for k = 1:nEpochs + for loc = 1:nLoc + ti = ((k-1)*dt)+tStart; + tf = k*dt+tStart-1; + data = filter_handle(... + time_series(loc, ti:tf), fs, cf(j), ... + cf(j+1)); + data = data'; + if strcmpi(outTypes(c), "Variance") + sa.data(j, k, loc) = var(data); + vargout = 0; + elseif strcmpi(outTypes(c), ... + "Standard_deviation") + sa.data(j, k, loc) = std(data); + vargout = 0; + elseif strcmpi(outTypes(c), "Mean") + sa.data(j, k, loc) = mean(data); + vargout = 0; + elseif strcmpi(outTypes(c), "Median") + sa.data(j, k, loc) = median(data); + vargout = 0; + elseif strcmpi(outTypes(c), "Kurtosis") + sa.data(j, k, loc) = kurtosis(data); + vargout = 0; + elseif strcmpi(outTypes(c), "Skewness") + sa.data(j, k, loc) = skewness(data); + vargout = 0; + end + end + end + end + outDir = path_check(subdir(inDir, outTypes(c))); + name = split(cases(i).name, '\'); + if length(name) == 1 + name = split(cases(i).name, '\'); + end + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + save(fullfile_check(filename), 'sa'); + catch + end %end try + waitbar((i+(c-1)*length(cases))/... + (length(cases)*length(outTypes)), f) + end + end + close(f) end \ No newline at end of file diff --git a/Measures/time_entropy.m b/Measures/time_entropy.m index a15a045..8c09791 100644 --- a/Measures/time_entropy.m +++ b/Measures/time_entropy.m @@ -1,163 +1,167 @@ -%% time_entropy -% This function computes the entropy related to EEG channels or ROIs. -% -% time_entropy(fs, cf, nEpochs, dt, inDir, outDirs, tStart, outTypes, ... -% filter_name, m, r) -% -% input: -% fs is the sampling frequency -% cf is an array containing the cut frequencies (es, [1 40]) -% nEpochs contains the number of epochs to compute -% dt contains the time (in seconds) of each epoch -% inDir is the directory containing each case -% tStart is the starting time (in seconds) to computate the first sample -% of the first epoch -% outTypes is the list of entropy types to compute (sample_entropy, -% approximate entropy or discretized entropy) -% filter_name is the name of the filtering function (athena_filter as -% default) -% m is the embedding dimension (its maximum value is the number of -% samples of the time series minus 1, 2 by default) -% r is the fraction of standard deviation of the time series, which have -% to be used as tolerance (0.2 by default) - - -function time_entropy(fs, cf, nEpochs, dt, inDir, tStart, outTypes, ... - filter_name, m, r) - - if nargin < 6 - tStart = 0; - end - if nargin < 7 - outTypes = ""; - end - if nargin < 8 - filter_name = 'athena_filter'; - end - if nargin < 9 - m = 2; - end - if nargin < 10 - r = 0.2; - end - - f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); - fchild = allchild(f); - fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) - fchild(1).JavaPeer.setStringPainted(true) - askList = 'Insert the entropy measures separated by a comma'; - - filter_handle = eval(strcat('@', filter_name)); - dt = fs*dt; - tStart = tStart*fs+1; - nBands = length(cf)-1; - - SampEn_names = ["sample_entropy", "SampEn", "SaEn"]; - AppEn_names = ["approximate_entropy", "AppEn", "ApEn"]; - DiscEn_names = ["discretized_entropy", "DiscEn", "DE"]; - - inDir = path_check(inDir); - cases = define_cases(inDir); - outTypes = string(outTypes); - - ctrl = 0; - if length(outTypes) == 1 - if strcmp(outTypes,'') - ctrl = 1; - end - end - - while ctrl == 1 - outTypes = []; - sup = string(inputdlg(askList)); - if strcmp(sup, 'exit') - return - end - sup = split(sup, ','); - for i = 1:length(sup) - if contains(sup(i, 1), SampEn_names) - outTypes = [outTypes, "sample_entropy"]; - elseif contains(sup(i, 1), DiscEn_names) - outTypes = [outTypes, "discretized_entropy"]; - elseif contains(sup(i, 1), AppEn_names) - outTypes = [outTypes, "approximate_entropy"]; - end - end - if length(outDirs) == length(outTypes) - ctrl = 0; - end - end - - for i = 1:length(outTypes) - if contains(outTypes(i), SampEn_names) - outTypes(i) = "sample_entropy"; - elseif contains(outTypes(i), DiscEn_names) - outTypes(i) = "discretized_entropy"; - elseif contains(outTypes(i), AppEn_names) - outTypes(i) = "approximate_entropy"; - end - end - - [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - check = check_filtering(time_series, dt, tStart, fs, cf, ... - filter_handle); - if check - close(f) - return; - end - - for c = 1:length(outTypes) - for i = 1:length(cases) - try - [time_series, fsOld, locations, chanlocs] = ... - load_data(strcat(inDir, cases(i).name), 1); - time_series = resample_signal(time_series, fs, fsOld); - nLoc = size(time_series, 1); - en.data = zeros(nBands, nEpochs, nLoc); - en.locations = locations; - en.chanlocs = chanlocs; - for j = 1:nBands - for k = 1:nEpochs - for loc = 1:nLoc - ti = ((k-1)*dt)+tStart; - tf = k*dt+tStart-1; - data = filter_handle(... - time_series(loc, ti:tf), fs, cf(j), ... - cf(j+1)); - data = data'; - if strcmpi(outTypes(c), "sample_entropy") - en.data(j, k, loc) = ... - sample_entropy(data, m, r); - elseif strcmpi(outTypes(c), ... - "discretized_entropy") - en.data(j, k, loc) = ... - discretized_entropy(data); - elseif strcmpi(outTypes(c), ... - "approximate_entropy") - en.data(j, k, loc) = ... - approximate_entropy(data, m, r); - end - end - end - end - outDir = path_check(subdir(inDir, outTypes(c))); - name = split(cases(i).name, '\'); - if length(name) == 1 - name = split(cases(i).name, '\'); - end - if length(name) > 1 - name = name{2}; - else - name = cases(i).name; - end - filename = strcat(outDir, strtok(name, '.'), '.mat'); - save(fullfile_check(filename), 'en'); - catch - end %end try - waitbar((i+(c-1)*length(cases))/... - (length(cases)*length(outTypes)), f) - end - end - close(f) +%% time_entropy +% This function computes the entropy related to EEG channels or ROIs. +% +% time_entropy(fs, cf, nEpochs, dt, inDir, outDirs, tStart, outTypes, ... +% filter_name, m, r) +% +% input: +% fs is the sampling frequency +% cf is an array containing the cut frequencies (es, [1 40]) +% nEpochs contains the number of epochs to compute +% dt contains the time (in seconds) of each epoch +% inDir is the directory containing each case +% tStart is the starting time (in seconds) to computate the first sample +% of the first epoch +% outTypes is the list of entropy types to compute (sample_entropy, +% approximate entropy or discretized entropy) +% filter_name is the name of the filtering function (athena_filter as +% default) +% m is the embedding dimension (its maximum value is the number of +% samples of the time series minus 1, 2 by default) +% r is the fraction of standard deviation of the time series, which have +% to be used as tolerance (0.2 by default) + + +function vargout = time_entropy(fs, cf, nEpochs, dt, inDir, tStart, outTypes, ... + filter_name, m, r) + + if nargin < 6 + tStart = 0; + end + if nargin < 7 + outTypes = ""; + end + if nargin < 8 + filter_name = 'athena_filter'; + end + if nargin < 9 + m = 2; + end + if nargin < 10 + r = 0.2; + end + + f = waitbar(0, 'Processing your data', 'Color', '[1 1 1]'); + fchild = allchild(f); + fchild(1).JavaPeer.setForeground(fchild(1).JavaPeer.getBackground.BLUE) + fchild(1).JavaPeer.setStringPainted(true) + askList = 'Insert the entropy measures separated by a comma'; + + vargout = -1; + filter_handle = eval(strcat('@', filter_name)); + dt = fs*dt; + tStart = tStart*fs+1; + nBands = length(cf)-1; + + SampEn_names = ["sample_entropy", "SampEn", "SaEn"]; + AppEn_names = ["approximate_entropy", "AppEn", "ApEn"]; + DiscEn_names = ["discretized_entropy", "DiscEn", "DE"]; + + inDir = path_check(inDir); + cases = define_cases(inDir); + outTypes = string(outTypes); + + ctrl = 0; + if length(outTypes) == 1 + if strcmp(outTypes,'') + ctrl = 1; + end + end + + while ctrl == 1 + outTypes = []; + sup = string(inputdlg(askList)); + if strcmp(sup, 'exit') + return + end + sup = split(sup, ','); + for i = 1:length(sup) + if contains(sup(i, 1), SampEn_names) + outTypes = [outTypes, "sample_entropy"]; + elseif contains(sup(i, 1), DiscEn_names) + outTypes = [outTypes, "discretized_entropy"]; + elseif contains(sup(i, 1), AppEn_names) + outTypes = [outTypes, "approximate_entropy"]; + end + end + if length(outDirs) == length(outTypes) + ctrl = 0; + end + end + + for i = 1:length(outTypes) + if contains(outTypes(i), SampEn_names) + outTypes(i) = "sample_entropy"; + elseif contains(outTypes(i), DiscEn_names) + outTypes(i) = "discretized_entropy"; + elseif contains(outTypes(i), AppEn_names) + outTypes(i) = "approximate_entropy"; + end + end + + [time_series, fsOld] = load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + check = check_filtering(time_series, dt, tStart, fs, cf, ... + filter_handle); + if check + close(f) + return; + end + + for c = 1:length(outTypes) + for i = 1:length(cases) + try + [time_series, fsOld, locations, chanlocs] = ... + load_data(strcat(inDir, cases(i).name), 1); + time_series = resample_signal(time_series, fs, fsOld); + nLoc = size(time_series, 1); + en.data = zeros(nBands, nEpochs, nLoc); + en.locations = locations; + en.chanlocs = chanlocs; + for j = 1:nBands + for k = 1:nEpochs + for loc = 1:nLoc + ti = ((k-1)*dt)+tStart; + tf = k*dt+tStart-1; + data = filter_handle(... + time_series(loc, ti:tf), fs, cf(j), ... + cf(j+1)); + data = data'; + if strcmpi(outTypes(c), "sample_entropy") + en.data(j, k, loc) = ... + sample_entropy(data, m, r); + vargout = 0; + elseif strcmpi(outTypes(c), ... + "discretized_entropy") + en.data(j, k, loc) = ... + discretized_entropy(data); + vargout = 0; + elseif strcmpi(outTypes(c), ... + "approximate_entropy") + en.data(j, k, loc) = ... + approximate_entropy(data, m, r); + vargout = 0; + end + end + end + end + outDir = path_check(subdir(inDir, outTypes(c))); + name = split(cases(i).name, '\'); + if length(name) == 1 + name = split(cases(i).name, '\'); + end + if length(name) > 1 + name = name{2}; + else + name = cases(i).name; + end + filename = strcat(outDir, strtok(name, '.'), '.mat'); + save(fullfile_check(filename), 'en'); + catch + end %end try + waitbar((i+(c-1)*length(cases))/... + (length(cases)*length(outTypes)), f) + end + end + close(f) end \ No newline at end of file diff --git a/README.md b/README.md index 1140df9..9ae4012 100644 --- a/README.md +++ b/README.md @@ -1,147 +1,153 @@ -# ATHENA -

- -

- -***ATHENA*** *(Automatic Toolbox for Handling Experimental Neural Analysis)* is a toolbox which allow to -automatically extract commonly analyzed measures, used to study neural time series, such as EEG and MEG. -These measures, which can be connectivity, power or background measures, can be studied through this toolbox, through correlation and -statistical analysis. - -The *pipeline* followed by this toolbox is easy and repetible. -It is possible to choose between a guided or a batch mode: -- The guided mode drives the user throw all the steps which compose the pipeline, allowing him to change parameters during the study -- The batch mode allows the user to compute the entire study automatically, chosing all the parameters before starting it - -

- -

- -The *pipeline* provided by *ATHENA* can be subdivided in 3 main steps: -1) **Feature Engineering**: in this step the user can choose the measure to extract and its parameters as sampling frequency, cut - frequencies to define the studied frequency bands, the number of epochs and the time window of each one and the starting time, extract - it, and finally compute the temporal and spatial average of the measure values and subdivide the studied subjects in their group -2) **Analysis**: in this step the user can choose the analysis to compute, between statistical and visual analysis, and their parameters -3) **Classification**: in this step the user can build a random forest, a decision tree or a multi-layer neural network classifier and its - performance, allowing to change parameters and to study the discriminant capability of the extracted measures - -

- -

- - -It is also possible to **display the signals**, easily switching between them if they are contained in the same directory. -The graphic interface allows: -- To *zoom* or *dezoom* the signal for a better visualization -- To choose the *length of the time window* to show -- To select a *specific time instant* to look at -- To *change the sampling frequency* (if it is not already present inside the signal files) -- To *select the locations* to show -- To *save a piece of signal*, by choosing the *starting time*, the *locations* and the *length of the time window* to save -- To make a *screenshot* of the displayed time window - -Furthermore, it is possible to *display the spectra* of the signals, and to execute a **time-frequency analysis** on them. - -

- -

- -[**Download Athena now to start you analysis!**](https://github.com/smlacava/Athena/archive/master.zip) - - -## Measure extraction -The measures extractably by ATHENA can be divided in: -- *Connectivity measures* - - **PLI** (Phase Lag Index) - - **PLV** (Phase Locking Value) - - **AEC** (Amplitude Envelope Correlation) - - **Corrected AEC** (Amplitude Envelope Correlation) - - **Coherence** - - **ICOH** (Cohenency) - - **Mutual Information** - - **Correlation coefficient** -- *Power measures* - - **relative PSD** (Power Spectral Density) -- *Background measures* - - **Exponent** - - **Offset** -- *Entropy measures* - - **Spectral Entropy** - - **Sample Entropy** - - **Approximate Entropy** - - **Discretized Entropy** -- *Autocorrelation measures* - - **Hurst exponent** -- *Statistical information* - - **Mean** - - **Median** - - **Standard deviation** - - **Skewness** - - **Kurtosis** - - **Variance** - -The *parameters* to choose to extract the measures are: -- **fs**: it is the sampling frequency, it can be automatically chosen if the parameter is present in the time series -- **cf**: it is the list of cut frequencies, they define each studied frequency band which will be define between two of the chosen cut - frequencies -- **epochs number**: it is the number of time epochs to study -- **epochs time**: it is the time window of each epoch -- **starting time**: it is the starting time of the first epoch, it is useful to avoid initial noise or altered signal due to preprocessing, - sources reconstruction, etc. -- **relative band**: it is used to extract the relative PSD as relative band - - -## Averaging and Grouping -In this step a table which contains the subjects names and their group (patients group or healthy controls group) is required. -If it is not present, the user can select the patients between the subject list in a graphic interface that will be create by *ATHENA* -if the button **␚** is pressed and save the resulting table to use it to compute the step. - - -## Analysis -For every analysis a file which contains the locations representing each row of the time series is requested. -This file allow to choose to compute the analysis on the measure pattern of *each location*, on the *asymmetry* between the pattern of the -right hemisphere and the left one, on the *global* pattern or on the *frontal, temporal, central, parietal and occipital areas*. -Furthermore, the user can select the conservativeness level which will define the alpha level for the correlation and statistical -analysis. -In these analysis, the results will be showed to the user through a p-values table and the significant results will be showed through a -table and in a graphical way. -It also possible to investigate the differences intra-subject of the values of a measure in the epochs, or to compute other visual analysis - - -The user can choose between various statistical analysis: -- **U Test**: the user can verify the presence of statistical differences between the patterns of the patients and the patterns of the - healthy controls -- **Index Correlation analysis**: the user can verify the correlation between a group pattern and an index corresponding to each - subject, an external file containing the index for each subject is required -- **Measures Correlation analysis**: the user can verify the presence of a correlation in a group between the patterns of two different - measures -- **Distribution analysis**: the user can compare the distribution of a measure of the group of patients with the one of the group of healty - controls -- **Histogram analysis**: the user can visually analyze the histogram related to the extracted measures -- **Descriptive statistics**: the user can compute some descriptive statistical measures, such as the median, the variance and the kurtosis, - and compare these value between the analyzed groups - - -Furthermore, the user can execute some visual analysis: -- **Scatter Plot analysis**: the user can visually analyze the scatter plot related to two measures, also with different spatial and - frequency parameters, or to different parameters of the same measure -- **Epochs Analysis**: the user can study the variation of a measure through the epochs in every frequency band for a subject -- **Network Measures**: the user can extract some network measures related to the connectivity measures - -

- -

- -## Classification -Finally, the user can train and test a **Decision Tree** classifier, a **Random Forest** classifier or a Multi-Layer **Neural Network** classifier in order -to verify the discriminant capability of the extracted features and to check if it is possible to distinguish between the two groups of subjects. - - -The user can also **merge the data of the significant results** of every previously computed statistical analysis in a csv file which -can be used for a classification or other external analysis. - - -## The wiki -For any doubt, you can consult the toolbox [**wiki**](https://github.com/smlacava/Athena/wiki/Home). - -> The toolbox is still under construction +# ATHENA +

+ +

+ +***ATHENA*** *(Automatic Toolbox for Handling Experimental Neural Analysis)* is a toolbox which allow to +automatically extract commonly analyzed measures, used to study neural time series, such as EEG and MEG. +These measures, which can be connectivity, power or background measures, can be studied through this toolbox, through correlation and +statistical analysis. + +The *pipeline* followed by this toolbox is easy and repetible. +It is possible to choose between a guided or a batch mode: +- The guided mode drives the user throw all the steps which compose the pipeline, allowing him to change parameters during the study +- The batch mode allows the user to compute the entire study automatically, chosing all the parameters before starting it + +

+ +

+ +The *pipeline* provided by *ATHENA* can be subdivided in 3 main steps: +1) **Feature Engineering**: in this step the user can choose the measure to extract and its parameters as sampling frequency, cut + frequencies to define the studied frequency bands, the number of epochs and the time window of each one and the starting time, extract + it, and finally compute the temporal and spatial average of the measure values and subdivide the studied subjects in their group +2) **Analysis**: in this step the user can choose the analysis to compute, between statistical and visual analysis, and their parameters +3) **Classification**: in this step the user can build a random forest, a decision tree or a multi-layer neural network classifier and its + performance, allowing to change parameters and to study the discriminant capability of the extracted measures + +

+ +

+ + +It is also possible to **display the signals**, easily switching between them if they are contained in the same directory. +The graphic interface allows: +- To *zoom* or *dezoom* the signal for a better visualization +- To choose the *length of the time window* to show +- To select a *specific time instant* to look at +- To *change the sampling frequency* (if it is not already present inside the signal files) +- To *select the locations* to show +- To *save a piece of signal*, by choosing the *starting time*, the *locations* and the *length of the time window* to save +- To make a *screenshot* of the displayed time window + +Furthermore, it is possible to *display the spectra* of the signals, and to execute a **time-frequency analysis** on them. + +

+ +

+ +[**Download Athena now to start you analysis!**](https://github.com/smlacava/Athena/archive/master.zip) + + +## Measure extraction +The measures extractably by ATHENA can be divided in: +- *Connectivity measures* + - **PLI** (Phase Lag Index) + - **PLV** (Phase Locking Value) + - **AEC** (Amplitude Envelope Correlation) + - **Corrected AEC** (Amplitude Envelope Correlation) + - **Coherence** + - **ICOH** (Cohenency) + - **Mutual Information** + - **Correlation coefficient** +- *Power measures* + - **relative PSD** (Power Spectral Density) +- *Background measures* + - **Exponent** + - **Offset** +- *Entropy measures* + - **Spectral Entropy** + - **Sample Entropy** + - **Approximate Entropy** + - **Discretized Entropy** +- *Autocorrelation measures* + - **Hurst exponent** +- *Statistical information* + - **Mean** + - **Median** + - **Standard deviation** + - **Skewness** + - **Kurtosis** + - **Variance** + +The *parameters* to choose to extract the measures are: +- **fs**: it is the sampling frequency, it can be automatically chosen if the parameter is present in the time series +- **cf**: it is the list of cut frequencies, they define each studied frequency band which will be define between two of the chosen cut + frequencies +- **epochs number**: it is the number of time epochs to study +- **epochs time**: it is the time window of each epoch +- **starting time**: it is the starting time of the first epoch, it is useful to avoid initial noise or altered signal due to preprocessing, + sources reconstruction, etc. +- **relative band**: it is used to extract the relative PSD as relative band + + +## Averaging and Grouping +In this step a table which contains the subjects names and their group (patients group or healthy controls group) is required. +If it is not present, the user can select the patients between the subject list in a graphic interface that will be create by *ATHENA* +if the button **␚** is pressed and save the resulting table to use it to compute the step. + + +## Analysis +For every analysis a file which contains the locations representing each row of the time series is requested. +This file allow to choose to compute the analysis on the measure pattern of *each location*, on the *asymmetry* between the pattern of the +right hemisphere and the left one, on the *global* pattern or on the *frontal, temporal, central, parietal and occipital areas*. +Furthermore, the user can select the conservativeness level which will define the alpha level for the correlation and statistical +analysis. +In these analysis, the results will be showed to the user through a p-values table and the significant results will be showed through a +table and in a graphical way. +It also possible to investigate the differences intra-subject of the values of a measure in the epochs, or to compute other visual analysis + + +The user can choose between various statistical analysis: +- **U Test**: the user can verify the presence of statistical differences between the patterns of the patients and the patterns of the + healthy controls +- **Index Correlation analysis**: the user can verify the correlation between a group pattern and an index corresponding to each + subject, an external file containing the index for each subject is required +- **Measures Correlation analysis**: the user can verify the presence of a correlation in a group between the patterns of two different + measures +- **Distribution analysis**: the user can compare the distribution of a measure of the group of patients with the one of the group of healty + controls +- **Histogram analysis**: the user can visually analyze the histogram related to the extracted measures +- **Descriptive statistics**: the user can compute some descriptive statistical measures, such as the median, the variance and the kurtosis, + and compare these value between the analyzed groups + + +Furthermore, the user can execute some visual analysis: +- **Scatter Plot analysis**: the user can visually analyze the scatter plot related to two measures, also with different spatial and + frequency parameters, or to different parameters of the same measure +- **Epochs Analysis**: the user can study the variation of a measure through the epochs in every frequency band for a subject +- **Network Measures**: the user can extract some network measures related to the connectivity measures + +

+ +

+ +## Classification +Finally, the user can train and test a **Decision Tree** classifier, a **Random Forest** classifier or a Multi-Layer **Neural Network** classifier in order +to verify the discriminant capability of the extracted features and to check if it is possible to distinguish between the two groups of subjects. + + +The user can also **merge the data of the significant results** of every previously computed statistical analysis in a csv file which +can be used for a classification or other external analysis. + + +## The wiki +For any doubt, you can consult the toolbox [**wiki**](https://github.com/smlacava/Athena/wiki/Home). + +> The toolbox is still under construction + + +## Test +To execute tests, run the command **runtests(directory,'IncludeSubfolders',true)** with the address of the folder **Test** as **directory** or run **runtest('IncludeSubfolders',true)** with the tests directory as **current folder**. + +> The toolbox is still under construction \ No newline at end of file diff --git a/Test/MeasuresTest/Example/signal.mat b/Test/MeasuresTest/Example/signal.mat new file mode 100644 index 0000000..a11234e Binary files /dev/null and b/Test/MeasuresTest/Example/signal.mat differ diff --git a/Test/MeasuresTest/FOOOFer_Test.m b/Test/MeasuresTest/FOOOFer_Test.m new file mode 100644 index 0000000..044765a --- /dev/null +++ b/Test/MeasuresTest/FOOOFer_Test.m @@ -0,0 +1,219 @@ +%This unit test the autocorrelation_measures function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the autocorrelation_measures function for different input +%parameters. Change outTypes for different entropy measures +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%outTyeps type of measure + + +function tests = FOOOFer_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + for x = {'Offset', 'Exponent'} + + try + delete(append(dir, x{1}, '\*')); + rmdir(append(dir, x{1})); + catch + end + + end +end + + +function teardownOnce(testCase) + + global dir + + try + delete(append(dir, '1\*')); + delete(append(dir, 'wrong\*')); + rmdir(append(dir, '1')); + rmdir(append(dir, 'wrong')); + catch + end + + clear + close all + +end + + +function teardown(testCase) + + global dir + + for x = {'Offset', 'Exponent'} + + delete(append(dir, x{1}, '\*')); + + try + rmdir(append(dir, x{1})); + catch + end + + end + +end + + +function testoutTypes_FOOOFer(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = FOOOFer(fs, cf, n_ep, dt, dir, t_start, 1); + verifyEqual(testCase, resp, -1); + bk = loadFromDisk(append(dir, '1\', file_name)); + verifyEmpty(testCase, bk); %verify bad out type, data should be empty + + resp = FOOOFer(fs, cf, n_ep, dt, dir, t_start,'wrong'); + verifyEqual(testCase, resp, -1); + bk = loadFromDisk(append(dir, 'wrong\', file_name)); + verifyEmpty(testCase, bk); %verify wrong out type, data should be empty + +end + + +function testTimeAndEpochs_offset(testCase) + + global fs cf dir file_name + + resp = FOOOFer(fs, cf, 5, 20, dir, 0,'OFF'); + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); %verify too many epochs, should be empty + + resp = FOOOFer(fs, cf, 1, 70, dir, 0,'OFF'); + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); %verify too long windows, should be empty + + resp = FOOOFer(fs, cf, 3, 20, dir, 70, 'OFF'); + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_offset(testCase) + + global cf n_ep dt t_start dir file_name + + resp = FOOOFer(0, cf, n_ep, dt, dir, t_start, 'OFF'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); + + resp = FOOOFer([], cf, n_ep, dt, dir, t_start, 'OFF'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); + +end + + +function testCutFrequencies_offset(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = FOOOFer(fs, [0], n_ep, dt, dir, t_start, 'OFF'); + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); %verify cut frequencies zero, data should be empty + + resp = FOOOFer(fs, [], n_ep, dt, dir, t_start, 'OFF'); + verifyEqual(testCase, resp, -1); + off = loadFromDisk(append(dir, 'Offset\', file_name)); + verifyEmpty(testCase, off); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_exponent(testCase) + + global fs cf dir file_name + + resp = FOOOFer(fs, cf, 5, 20, dir, 0,'EXP'); + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); %verify too many epochs, should be empty + + resp = FOOOFer(fs, cf, 1, 70, dir, 0,'EXP'); + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); %verify too long windows, should be empty + + resp = FOOOFer(fs, cf, 3, 20, dir, 70, 'EXP'); + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_exponent(testCase) + + global cf n_ep dt t_start dir file_name + + resp = FOOOFer(0, cf, n_ep, dt, dir, t_start, 'EXP'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); + + resp = FOOOFer([], cf, n_ep, dt, dir, t_start, 'EXP'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); + +end + + +function testCutFrequencies_exponent(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = FOOOFer(fs, [0], n_ep, dt, dir, t_start, 'EXP'); + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); %verify cut frequencies zero, data should be empty + + resp = FOOOFer(fs, [], n_ep, dt, dir, t_start, 'EXP'); + verifyEqual(testCase, resp, -1); + exp = loadFromDisk(append(dir, 'Exponent\', file_name)); + verifyEmpty(testCase, exp); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end \ No newline at end of file diff --git a/Test/MeasuresTest/PSDr_Test.m b/Test/MeasuresTest/PSDr_Test.m new file mode 100644 index 0000000..41e7436 --- /dev/null +++ b/Test/MeasuresTest/PSDr_Test.m @@ -0,0 +1,140 @@ +%This unit test the PSDr function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the PSDr function for different input +%parameters +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%band frequency bands + +function tests = PSDr_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + try + delete(append(dir, 'PSDr\*')); + rmdir(append(dir, 'PSDr')); + catch + end + +end + + +function teardownOnce(testCase) + + global dir + + try + rmdir(append(dir, 'PSDr')); + catch + end + + clear + close all +end + + +function teardown(testCase) + + global dir + + delete(append(dir, 'PSDr\*')); + +end + + +function testBandPSDr(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = PSDr(fs, cf, n_ep, dt, dir, t_start, []); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify empty band, data should be empty + + resp = PSDr(fs, cf, n_ep, dt, dir, t_start, [0]); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + PSDr(fs, cf, n_ep, dt, dir, t_start, [0 0]); + exp = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEqual(testCase, psdr, exp); %verify zero band, data should be like a single value band + + resp = PSDr(fs, cf, 3, 20, dir, 70, [50]); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + PSDr(fs, cf, n_ep, dt, dir, t_start, [50 50]); + exp = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEqual(testCase, psdr, exp); %verify bad band, data should be like a single value band + +end + +function testTimeAndEpochsPSDr(testCase) + + global fs cf dir band file_name + + resp = PSDr(fs, cf, 5, 20, dir, 0, band); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify too many epochs, data should be empty + + resp = PSDr(fs, cf, 1, 70, dir, 0, band); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify too long windows, should be empty + + resp = PSDr(fs, cf, 3, 20, dir, 70, band); + verifyEqual(testCase, resp, 0); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify too long windows, should be a matrix of zeros +end + + +function testSamplingFrequencyPSDr(testCase) + + global cf n_ep dt t_start dir band file_name + + resp = PSDr(0, cf, n_ep, dt, dir, t_start, band); + verifyEqual(testCase, resp, -1); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify sampling frequency zero, data should be empty + +end + + +function testCutFrequenciesPSDr(testCase) + + global fs n_ep dt t_start dir band series file_name + + resp = PSDr(fs, [0], n_ep, dt, dir, t_start, band); + verifyEqual(testCase, resp, -1); + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify cut frequencies null, data should be empty + + resp = PSDr(fs, [], n_ep, dt, dir, t_start, band); + verifyEqual(testCase, resp, -1) + psdr = loadFromDisk(append(dir, 'PSDr\', file_name)); + verifyEmpty(testCase, psdr); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % %how to test function on single sample + % end + %end +end \ No newline at end of file diff --git a/Test/MeasuresTest/autocorrelation_measures_Test.m b/Test/MeasuresTest/autocorrelation_measures_Test.m new file mode 100644 index 0000000..efb3a97 --- /dev/null +++ b/Test/MeasuresTest/autocorrelation_measures_Test.m @@ -0,0 +1,162 @@ +%This unit test the autocorrelation_measures function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the autocorrelation_measures function for different input +% parameters. Change outTypes for different entropy measures +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%outTyeps type of measure +%filter_name name of the filtering functio + + +function tests = autocorrelation_measures_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + try + delete(append(dir, 'Hurst\*')); + rmdir(append(dir, 'Hurst')); + catch + end + +end + + +function teardownOnce(testCase) + + global dir + + clear + close all + + try + rmdir(append(dir, 'Hurst')); + catch + end + + try + delete(append(dir, '1\*')); + delete(append(dir, 'wrong\*')); + rmdir(append(dir, '1')); + rmdir(append(dir, 'wrong')); + catch + end + +end + + +function teardown(testCase) + + global dir + + delete(append(dir, 'Hurst\*')); + +end + + +function testoutTypes_autocorrelation_measures(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = autocorrelation_measures(fs, cf, n_ep, dt, dir, t_start, 1); + verifyEqual(testCase, resp, -1); + am = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, am); %verify bad out type, data should be empty + + resp = autocorrelation_measures(fs, cf, n_ep, dt, dir, t_start,'wrong'); + verifyEqual(testCase, resp, -1); + am = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, am); %verify wrong out type, data should be empty + +end + + +function testFilterName_autocorrelation_measures(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = autocorrelation_measures(fs, cf, n_ep, dt, dir, t_start, 'Hurst', 'wrong'); + verifyEqual(testCase, resp, -1); + am = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, am); %verify wrong filter name, data should be empty + +end + + +function testTimeAndEpochs_hurst_index(testCase) + + global fs cf dir file_name + + resp = autocorrelation_measures(fs, cf, 5, 20, dir, 0,'Hurst'); + verifyEqual(testCase, resp, 0); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); %verify too many epochs, data should be empty + + resp = autocorrelation_measures(fs, cf, 1, 70, dir, 0,'Hurst'); + verifyEqual(testCase, resp, -1); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifympty(testCase, he); %verify too long windows, should be empty + + resp = autocorrelation_measures(fs, cf, 3, 20, dir, 70, 'Hurst'); + verifyEqual(testCase, resp, 0); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_hurst_index(testCase) + + global cf n_ep dt t_start dir + + resp = autocorrelation_measures(0, cf, n_ep, dt, dir, t_start, 'Hurst'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, 0); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); + + resp = autocorrelation_measures([], cf, n_ep, dt, dir, t_start, 'Hurst'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, 0); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); + +end + + +function testCutFrequencies_hurst_index(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = autocorrelation_measures(fs, [0], n_ep, dt, dir, t_start, 'Hurst'); + verifyEqual(testCase, resp, -1); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); %verify cut frequencies zero, data should be empty + + resp = autocorrelation_measures(fs, [], n_ep, dt, dir, t_start, 'Hurst'); + verifyEqual(testCase, resp, 0); + he = loadFromDisk(append(dir, 'Hurst\', file_name)); + verifyEmpty(testCase, he); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end \ No newline at end of file diff --git a/Test/MeasuresTest/connectivity_Test.m b/Test/MeasuresTest/connectivity_Test.m new file mode 100644 index 0000000..7d90e7d --- /dev/null +++ b/Test/MeasuresTest/connectivity_Test.m @@ -0,0 +1,670 @@ +%This unit test the connectivity function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the connectivity function for different input +%parameters. Change outTypes for different entropy measures +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%outTyeps type of measure +%filter_name name of the filter + +function tests = connectivity_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + for x = {'PLI', 'PLV', 'AEC', 'AECo', 'Coherence', 'ICOH', 'correlation_coefficient', 'wPLI', 'mutual_information'} + + try + delete(append(dir, x{1}, '\*')); + rmdir(append(dir, x{1})); + catch + end + + end +end + + +function teardownOnce(testCase) + + global dir + + try + delete(append(dir, '1\*')); + delete(append(dir, 'wrong\*')); + rmdir(append(dir, '1')); + rmdir(append(dir, 'wrong')); + catch + end + + clear + close all + +end + + +function teardown(testCase) + + global dir + + for x = {'PLI', 'PLV', 'AEC', 'AECo', 'Coherence', 'ICOH', 'correlation_coefficient', 'wPLI', 'mutual_information'} + + delete(append(dir, x{1}, '\*')); + + try + rmdir(append(dir, x{1})); + catch + end + + end + +end + + +function testoutTypes_connectivity(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = connectivity(fs, cf, n_ep, dt, dir, t_start, 1); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, '1\', file_name)); + verifyEqual(testCase, co, zeros(length(co(:, 1)), length(co(1, :)))); %verify bad out type, data should be zero + + resp = connectivity(fs, cf, n_ep, dt, dir, t_start,'wrong'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'wrong\', file_name)); + verifyEqual(testCase, co, zeros(length(co(:, 1)), length(co(1, :)))); %verify wrong out type, data should be zero + +end + +function testFilterName_connectivity(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = connectivity(fs, cf, n_ep, dt, dir, t_start,'PLI', 'wrong'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, co); %verify wrong filter name, data should be empty + +end + +function testTimeAndEpochs_PLI(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'PLI'); + verifyEqual(testCase, resp, 0); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); %verify too many epochs, should be empty + + resp = connectivity(fs, cf, 1, 70, dir, 0,'PLI'); + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); %verify too long windows, should be empty + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'PLI'); + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_PLI(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'PLI'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'PLI'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); + +end + + +function testCutFrequencies_PLI(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'PLI'); + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'PLI'); + verifyEqual(testCase, resp, -1); + pli = loadFromDisk(append(dir, 'PLI\', file_name)); + verifyEmpty(testCase, pli); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_PLV(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'PLV'); + verifyEqual(testCase, resp, 0); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); %verify too many epochs, should be empty + + resp = connectivity(fs, cf, 1, 70, dir, 0,'PLV'); + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); %verify too long windows, should be empty + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'PLV'); + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_PLV(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'PLV'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'PLV'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); + +end + + +function testCutFrequencies_PLV(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'PLV'); + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'PLV'); + verifyEqual(testCase, resp, -1); + plv = loadFromDisk(append(dir, 'PLV\', file_name)); + verifyEmpty(testCase, plv); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_AEC(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'AEC'); + verifyEqual(testCase, resp, 0); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); %verify too many epochs, should be empty + + resp = connectivity(fs, cf, 1, 70, dir, 0,'AEC'); + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); %verify too long windows, should be empty + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'AEC'); + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_AEC(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'AEC'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'AEC'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); + +end + + +function testCutFrequencies_AEC(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'AEC'); + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'AEC'); + verifyEqual(testCase, resp, -1); + aec = loadFromDisk(append(dir, 'AEC\', file_name)); + verifyEmpty(testCase, aec); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_AECo(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'AECo'); + verifyEqual(testCase, resp, 0); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); %verify too many epochs, should be empty + + resp = connectivity(fs, cf, 1, 70, dir, 0,'AECo'); + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); %verify too long windows, should be empty + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'AECo'); + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_AECo(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'AECo'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'AECo'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); + +end + + +function testCutFrequencies_AECo(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'AECo'); + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'AECo'); + verifyEqual(testCase, resp, -1); + aeco = loadFromDisk(append(dir, 'AECo\', file_name)); + verifyEmpty(testCase, aeco); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_coherence(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'Coherence'); + verifyEqual(testCase, resp, 0); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); %verify too many epochs, should be empty + + resp = connectivity(fs, cf, 1, 70, dir, 0,'Coherence'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); %verify too long windows, should be empty + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'Coherence'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); %verify too long windows, should be empty + +end + + +function testSamplingFrequency_coherence(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'Coherence'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'Coherence'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); + +end + + +function testCutFrequencies_coherence(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'Coherence'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'Coherence'); + verifyEqual(testCase, resp, -1); + co = loadFromDisk(append(dir, 'Coherence\', file_name)); + verifyEmpty(testCase, co); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_ICOH(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'ICOH'); + verifyEqual(testCase, resp, 0); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); %verify too many epochs, last column should be zero + + resp = connectivity(fs, cf, 1, 70, dir, 0,'ICOH'); + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); %verify too long windows, single column should be zero + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'ICOH'); + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); %verify too long windows, should be a matrix of zeros + +end + + +function testSamplingFrequency_ICOH(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'ICOH'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'ICOH'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); + +end + + +function testCutFrequencies_ICOH(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'ICOH'); + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'ICOH'); + verifyEqual(testCase, resp, -1); + ico = loadFromDisk(append(dir, 'ICOH\', file_name)); + verifyEmpty(testCase, ico); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_correlation_coefficient(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'correlation_coefficient'); + verifyEqual(testCase, resp, 0); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); %verify too many epochs, last column should be zero + + resp = connectivity(fs, cf, 1, 70, dir, 0,'correlation_coefficient'); + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); %verify too long windows, single column should be zero + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'correlation_coefficient'); + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); %verify too long windows, should be a matrix of zeros + +end + + +function testSamplingFrequency_correlation_coefficient(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'correlation_coefficient'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'correlation_coefficient'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); + +end + + +function testCutFrequencies_correlation_coefficient(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'correlation_coefficient'); + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'correlation_coefficient'); + verifyEqual(testCase, resp, -1); + cc = loadFromDisk(append(dir, 'correlation_coefficient\', file_name)); + verifyEmpty(testCase, cc); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_wPLI(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'wPLI'); + verifyEqual(testCase, resp, 0); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); %verify too many epochs, last column should be zero + + resp = connectivity(fs, cf, 1, 70, dir, 0,'wPLI'); + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); %verify too long windows, single column should be zero + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'wPLI'); + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); %verify too long windows, should be a matrix of zeros + +end + + +function testSamplingFrequency_wPLI(testCase) + + global cf n_ep dt t_start dir file_name + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'wPLI'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'wPLI'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); + +end + + +function testCutFrequencies_wPLI(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'wPLI'); + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'wPLI'); + verifyEqual(testCase, resp, -1); + w = loadFromDisk(append(dir, 'wPLI\', file_name)); + verifyEmpty(testCase, w); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_mutual_information(testCase) + + global fs cf dir file_name + + resp = connectivity(fs, cf, 5, 20, dir, 0,'mutual_information'); + verifyEqual(testCase, resp, 0); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); %verify too many epochs, last column should be zero + + resp = connectivity(fs, cf, 1, 70, dir, 0,'mutual_information'); + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); %verify too long windows, single column should be zero + + resp = connectivity(fs, cf, 3, 20, dir, 70, 'mutual_information'); + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); %verify too long windows, should be a matrix of zeros + +end + + +function testSamplingFrequency_mutual_information(testCase) + + global cf n_ep dt t_start dir + + resp = connectivity(0, cf, n_ep, dt, dir, t_start, 'mutual_information'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); + + resp = connectivity([], cf, n_ep, dt, dir, t_start, 'mutual_information'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); + +end + + +function testCutFrequencies_mutual_information(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = connectivity(fs, [0], n_ep, dt, dir, t_start, 'mutual_information'); + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); %verify cut frequencies zero, data should be empty + + resp = connectivity(fs, [], n_ep, dt, dir, t_start, 'mutual_information'); + verifyEqual(testCase, resp, -1); + mi = loadFromDisk(append(dir, 'mutual_information\', file_name)); + verifyEmpty(testCase, mi); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end \ No newline at end of file diff --git a/Test/MeasuresTest/loadFromDisk.m b/Test/MeasuresTest/loadFromDisk.m new file mode 100644 index 0000000..bd6ce0d --- /dev/null +++ b/Test/MeasuresTest/loadFromDisk.m @@ -0,0 +1,14 @@ +%This function load data from disk +%Many Athena functions use disk instead of cache to save the big amount of data returned + +function data = loadFromDisk(file_dir) + + try + data = load(file_dir); + data = struct2cell(data); + data = data{1}.data; + catch + data = []; + end + +end \ No newline at end of file diff --git a/Test/MeasuresTest/load_test_parameters.m b/Test/MeasuresTest/load_test_parameters.m new file mode 100644 index 0000000..6ee095a --- /dev/null +++ b/Test/MeasuresTest/load_test_parameters.m @@ -0,0 +1,16 @@ +%This function load the parameters used for the test from a csv file +%Edit the csv file and this function to add parameters + +function [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters(file) + + [cf,n_ep,dt,t_start,band,dir,file_name] = readvars(file); + cf = cf'; + n_ep = n_ep(1); + dt = dt(1); + t_start = t_start(1); + band = band'; + band = band(1:2); + dir = dir{1, 1}; + file_name = file_name{1, 1}; + +end \ No newline at end of file diff --git a/Test/MeasuresTest/parameters.csv b/Test/MeasuresTest/parameters.csv new file mode 100644 index 0000000..8b4c3d0 --- /dev/null +++ b/Test/MeasuresTest/parameters.csv @@ -0,0 +1,7 @@ +cf;n_ep;dt;t_start;band;dir;file_name +1;3;20;0;1;C:\Users\Piermario Zoroddu\Desktop\Athena\Test\MeasuresTest\Example\;signal.mat +4;;;;40;; +5;;;;;; +13;;;;;; +30;;;;;; +40;;;;;; diff --git a/Test/MeasuresTest/spectral_entropy_Test.m b/Test/MeasuresTest/spectral_entropy_Test.m new file mode 100644 index 0000000..813c982 --- /dev/null +++ b/Test/MeasuresTest/spectral_entropy_Test.m @@ -0,0 +1,122 @@ +%This unit test the spectral_entropy function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the spectral_entropy function for different input +%parameters +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time + +function tests = spectral_entropy_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + try + delete(append(dir, 'PEntropy\*')); + rmdir(append(dir, 'PEntropy')); + catch + end + +end + + +function teardownOnce(testCase) + + global dir + + try + rmdir(append(dir, 'PEntropy')); + catch + end + + clear + close all + +end + + +function teardown(testCase) + + global dir + + delete(append(dir, 'PEntropy\*')); + +end + + +function testTimeAndEpochs_spectral_entropy(testCase) + + global fs cf dir file_name + + resp = spectral_entropy(fs, cf, 5, 20, dir, 0); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); %verify too many epochs, data should be empty + + resp = spectral_entropy(fs, cf, 1, 70, dir, 0); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); %verify too long windows, should be empty + + resp = spectral_entropy(fs, cf, 3, 20, dir, 70); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); %verify too long windows, should be mepty + +end + + +function testSamplingFrequency_spectral_entropy(testCase) + + global cf n_ep dt t_start dir file_name + + resp = spectral_entropy(0, cf, n_ep, dt, dir, t_start); %verify sampling frequency zero, should be empty + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); + + resp = spectral_entropy([], cf, n_ep, dt, dir, t_start); %verify sampling frequency null, should be empty + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); + +end + + +function testCutFrequencies_spectral_entropy(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = spectral_entropy(fs, [0], n_ep, dt, dir, t_start); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); %verify cut frequencies zero, should be empty + + resp = spectral_entropy(fs, [], n_ep, dt, dir, t_start); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'PEntropy\', file_name)); + verifyEmpty(testCase, se); %verify cut frequencies empty, should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end \ No newline at end of file diff --git a/Test/MeasuresTest/statistical_information_Test.m b/Test/MeasuresTest/statistical_information_Test.m new file mode 100644 index 0000000..a61fbe2 --- /dev/null +++ b/Test/MeasuresTest/statistical_information_Test.m @@ -0,0 +1,479 @@ +%This unit test the statistical_information function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the statistical_information function for different input +%parameters. Change outTypes for different entropy measures +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%outTyeps type of measure +%filter_name name of the filter + +function tests = statistical_information_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + for x = {'Median', 'Mean', 'Standard_deviation', 'Variance', 'Skewness', 'Kurtosis'} + + try + delete(append(dir, x{1}, '\*')); + rmdir(append(dir, x{1})); + catch + end + + end +end + + +function teardownOnce(testCase) + + global dir + + try + delete(append(dir, '1\*')); + delete(append(dir, 'wrong\*')); + rmdir(append(dir, '1')); + rmdir(append(dir, 'wrong')); + catch + end + + clear + close all + +end + + +function teardown(testCase) + + global dir + + for x = {'Median', 'Mean', 'Standard_deviation', 'Variance', 'Skewness', 'Kurtosis'} + + delete(append(dir, x{1}, '\*')); + + try + rmdir(append(dir, x{1})); + catch + end + + end + +end + + +function testoutTypes_statistical_information(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = statistical_information(fs, cf, n_ep, dt, dir, t_start, 1); + verifyEqual(testCase, resp, -1); + sa = loadFromDisk(append(dir, '1\', file_name)); + verifyEqual(testCase, sa, zeros(length(sa(:, 1)), length(sa(1, :)))); %verify bad out type, data should be a matrix of zeros + + resp = statistical_information(fs, cf, n_ep, dt, dir, t_start,'wrong'); + verifyEqual(testCase, resp, -1); + sa = loadFromDisk(append(dir, 'wrong\', file_name)); + verifyEqual(testCase, sa, zeros(length(sa(:, 1)), length(sa(1, :)))); %verify wrong out type, data should be a matrix of zeros + +end + + +function testFilterName_statistical_information(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = statistical_information(fs, cf, n_ep, dt, dir, t_start,'Median', 'wrong'); + verifyEqual(testCase, resp, -1); + sa = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, sa); %verify wrong filter name, data should be empty + +end + + +function testTimeAndEpochs_median(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Median'); + verifyEqual(testCase, resp, 0); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Median'); + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Median'); + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_median(testCase) + + global cf n_ep dt t_start dir file_name + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Median'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Median'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); + +end + + +function testCutFrequencies_median(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Median'); + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Median'); + verifyEqual(testCase, resp, -1); + md = loadFromDisk(append(dir, 'Median\', file_name)); + verifyEmpty(testCase, md); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_mean(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Mean'); + verifyEqual(testCase, resp, 0); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Mean'); + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Mean'); + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_mean(testCase) + + global cf n_ep dt t_start dir file_name + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Mean'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Mean'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); + +end + + +function testCutFrequencies_mean(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Mean'); + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Mean'); + verifyEqual(testCase, resp, -1); + m = loadFromDisk(append(dir, 'Mean\', file_name)); + verifyEmpty(testCase, m); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_standard_deviation(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Standard_deviation'); + verifyEqual(testCase, resp, 0); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Standard_deviation'); + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Standard_deviation'); + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_standard_deviation(testCase) + + global cf n_ep dt t_start dir + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Standard_deviation'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Standard_deviation'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); + +end + + +function testCutFrequencies_standard_deviation(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Standard_deviation'); + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Standard_deviation'); + verifyEqual(testCase, resp, -1); + std = loadFromDisk(append(dir, 'Standard_deviation\', file_name)); + verifyEmpty(testCase, std); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_variance(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Variance'); + verifyEqual(testCase, resp, 0); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Variance'); + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Variance'); + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); %verify too long windows, data should be empty +end + + +function testSamplingFrequency_variance(testCase) + + global cf n_ep dt t_start dir file_name + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Variance'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Variance'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); + +end + + +function testCutFrequencies_variance(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Variance'); + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Variance'); + verifyEqual(testCase, resp, -1); + v = loadFromDisk(append(dir, 'Variance\', file_name)); + verifyEmpty(testCase, v); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end +end + + +function testTimeAndEpochs_skewness(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Skewness'); + verifyEqual(testCase, resp, 0); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Skewness'); + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Skewness'); + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_skewness(testCase) + + global cf n_ep dt t_start dir file_name + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Skewness'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Skewness'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); + +end + + +function testCutFrequencies_skewness(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Skewness'); + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Skewness'); + verifyEqual(testCase, resp, -1); + sk = loadFromDisk(append(dir, 'Skewness\', file_name)); + verifyEmpty(testCase, sk); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end +end + + +function testTimeAndEpochs_kurtosis(testCase) + + global fs cf dir file_name + + resp = statistical_information(fs, cf, 5, 20, dir, 0,'Kurtosis'); + verifyEqual(testCase, resp, 0); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); %verify too many epochs, data should be empty + + resp = statistical_information(fs, cf, 1, 70, dir, 0,'Kurtosis'); + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); %verify too long windows, data should be empty + + resp = statistical_information(fs, cf, 3, 20, dir, 70, 'Kurtosis'); + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_kurtosis(testCase) + + global cf n_ep dt t_start dir file_name + + resp = statistical_information(0, cf, n_ep, dt, dir, t_start, 'Kurtosis'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); + + resp = statistical_information([], cf, n_ep, dt, dir, t_start, 'Kurtosis'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); + +end + + +function testCutFrequencies_kurtosis(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = statistical_information(fs, [0], n_ep, dt, dir, t_start, 'Kurtosis'); + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); %verify cut frequencies zero, data should be empty + + resp = statistical_information(fs, [], n_ep, dt, dir, t_start, 'Kurtosis'); + verifyEqual(testCase, resp, -1); + ku = loadFromDisk(append(dir, 'Kurtosis\', file_name)); + verifyEmpty(testCase, ku); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end +end \ No newline at end of file diff --git a/Test/MeasuresTest/time_entropy_Test.m b/Test/MeasuresTest/time_entropy_Test.m new file mode 100644 index 0000000..d7d39ee --- /dev/null +++ b/Test/MeasuresTest/time_entropy_Test.m @@ -0,0 +1,296 @@ +%This unit test the time_entropy function. Recommended use a .mat signal named +%signal.mat (time serie index must be named data) in Example folder +%These functions test the output of the time_entropy function for different input +%parameters. Change outTypes for different entropy measures +%fs sample frequency +%cf cut frequencies +%n_ep number of epochs +%dt time window +%dir directory +%t_start starting time +%outTyeps type of measure +%filter_name name of the filtering functio +%m embedding dimension +%r fraction of standard deviation + +function tests = time_entropy_Test + + tests = functiontests(localfunctions); + +end + + +function setupOnce(testCase) + + global fs cf n_ep dt t_start band dir file_name series + + [cf,n_ep,dt,t_start,band,dir,file_name] = load_test_parameters("parameters.csv"); + + signal = load(append(dir, file_name)); + series = signal.data.time_series; %time series + fs = signal.data.fs; %sample frequency + + for x = {'discretized_entropy', 'sample_entropy', 'approximate_entropy'} + + try + delete(append(dir, x{1}, '\*')); + rmdir(append(dir, x{1})); + catch + end + + end +end + + +function teardownOnce(testCase) + + global dir + + try + delete(append(dir, '1\*')); + delete(append(dir, 'wrong\*')); + rmdir(append(dir, '1')); + rmdir(append(dir, 'wrong')); + catch + end + + clear + close all + +end + + +function teardown(testCase) + + global dir + + for x = {'discretized_entropy', 'sample_entropy', 'approximate_entropy'} + + delete(append(dir, x{1}, '\*')); + + try + rmdir(append(dir, x{1})); + catch + end + + end + +end + + +function testoutTypes_time_entropy(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = time_entropy(fs, cf, n_ep, dt, dir, t_start, 1); + verifyEqual(testCase, resp, -1); + te = loadFromDisk(append(dir, '1\', file_name)); + verifyEqual(testCase, te(: , end), zeros(length(te(:, end)), 1)); %verify bad out type, data should be a matrix of zeros + + resp = time_entropy(fs, cf, n_ep, dt, dir, t_start,'wrong'); + verifyEqual(testCase, resp, -1); + te = loadFromDisk(append(dir, 'wrong\', file_name)); + verifyEqual(testCase, te(: , end), zeros(length(te(:, end)), 1)); %verify wrong out type, data should be a matrix of zeros + +end + + +function testFilterName_time_entropy(testCase) + + global fs cf n_ep dt t_start dir file_name + + resp = time_entropy(fs, cf, n_ep, dt, dir, t_start,'discretized_entropy', 'wrong'); + verifyEqual(testCase, resp, -1); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify wrong filter name, data should be empty + +end + + +function testTimeAndEpochs_discretized_entropy(testCase) + + global fs cf dir file_name + + resp = time_entropy(fs, cf, 5, 20, dir, 0,'discretized_entropy'); + verifyEqual(testCase, resp, 0); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify too many epochs, last column should be empty + + resp = time_entropy(fs, cf, 1, 70, dir, 0,'discretized_entropy'); + verifyEqual(testCase, resp, -1); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify too long windows, data should be empty + + resp = time_entropy(fs, cf, 3, 20, dir, 70, 'discretized_entropy'); + verifyEqual(testCase, resp, -1); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_discretized_entropy(testCase) + + global cf n_ep dt t_start dir file_name + + resp = time_entropy(0, cf, n_ep, dt, dir, t_start, 'discretized_entropy'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, 0); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); + + resp = time_entropy([], cf, n_ep, dt, dir, t_start, 'discretized_entropy'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, 0); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); + +end + + +function testCutFrequencies_discretized_entropy(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = time_entropy(fs, [0], n_ep, dt, dir, t_start, 'discretized_entropy'); + verifyEqual(testCase, resp, -1); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify cut frequencies zero, data should be empty + + resp = time_entropy(fs, [], n_ep, dt, dir, t_start, 'discretized_entropy'); + verifyEqual(testCase, resp, -1); + de = loadFromDisk(append(dir, 'discretized_entropy\', file_name)); + verifyEmpty(testCase, de); %verify cut frequencies empty, data should be empty + + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_sample_entropy(testCase) + + global fs cf dir file_name + + resp = time_entropy(fs, cf, 5, 20, dir, 0,'sample_entropy'); + verifyEqual(testCase, resp, 0); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEmpty(testCase, se); %verify too many epochs, last column should be empty + + resp = time_entropy(fs, cf, 1, 70, dir, 0,'sample_entropy'); + verifyEqual(testCase, resp, 0); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEqual(testCase, se(: , 1), zeros(length(se(:, 1)), 1)); %verify too long windows, single column should be zero + + resp = time_entropy(fs, cf, 3, 20, dir, 70, 'sample_entropy'); + verifyEqual(testCase, resp, 0); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEqual(testCase, se, zeros(length(se(1, :)), length(se(:, 1)))); %verify too long windows, should be a matrix of zeros + +end + + +function testSamplingFrequency_sample_entropy(testCase) + + global cf n_ep dt t_start dir file_name + + resp = time_entropy(0, cf, n_ep, dt, dir, t_start, 'sample_entropy'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, 0); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEmpty(testCase, se); + + resp = time_entropy([], cf, n_ep, dt, dir, t_start, 'sample_entropy'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, 0); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEmpty(testCase, se); + +end + + +function testCutFrequencies_sample_entropy(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = time_entropy(fs, [0], n_ep, dt, dir, t_start, 'sample_entropy'); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEmpty(testCase, se); %verify cut frequencies zero, data should be empty + + resp = time_entropy(fs, [], n_ep, dt, dir, t_start, 'sample_entropy'); + verifyEqual(testCase, resp, -1); + se = loadFromDisk(append(dir, 'sample_entropy\', file_name)); + verifyEmpty(testCase, se); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + % for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + % + % end + %end + +end + + +function testTimeAndEpochs_approximate_entropy(testCase) + + global fs cf dir file_name + + resp = time_entropy(fs, cf, 5, 20, dir, 0,'approximate_entropy'); + verifyEqual(testCase, resp, 0); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); %verify too many epochs, data should be empty + + resp = time_entropy(fs, cf, 1, 70, dir, 0,'approximate_entropy'); + verifyEqual(testCase, resp, -1); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); %verify too long windows, data should be empty + + resp = time_entropy(fs, cf, 3, 20, dir, 70, 'approximate_entropy'); + verifyEqual(testCase, resp, -1); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); %verify too long windows, data should be empty + +end + + +function testSamplingFrequency_approximate_entropy(testCase) + + global cf n_ep dt t_start dir file_name + + resp = time_entropy(0, cf, n_ep, dt, dir, t_start, 'approximate_entropy'); %verify sampling frequency zero, data should be empty + verifyEqual(testCase, resp, 0); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); + + resp = time_entropy([], cf, n_ep, dt, dir, t_start, 'approximate_entropy'); %verify sampling frequency null, data should be empty + verifyEqual(testCase, resp, 0); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); + +end + + +function testCutFrequencies_approximate_entropy(testCase) + + global fs n_ep dt t_start dir series file_name + + resp = time_entropy(fs, [0], n_ep, dt, dir, t_start, 'approximate_entropy'); + verifyEqual(testCase, resp, -1); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); %verify cut frequencies zero, data should be empty + + resp = time_entropy(fs, [], n_ep, dt, dir, t_start, 'approximate_entropy'); + verifyEqual(testCase, resp, -1); + ae = loadFromDisk(append(dir, 'approximate_entropy\', file_name)); + verifyEmpty(testCase, ae); %verify cut frequencies empty, data should be empty + + %for k = 1:n_ep + %for j = 1:size(series, 1) + % data = squeeze(series(j, dt*(k-1)+1:k*dt)); + + %end + %end +end \ No newline at end of file