diff --git a/src/backroom/pspm_find_data_epochs.m b/src/backroom/pspm_find_data_epochs.m index 2ef36973a..8d7cd0e93 100644 --- a/src/backroom/pspm_find_data_epochs.m +++ b/src/backroom/pspm_find_data_epochs.m @@ -3,25 +3,25 @@ % of a channel. It then writes it to a file depending on the input and % returns the location of that file. % -% FORMAT: +% FORMAT: % [sts, epochfile] = pspm_find_data_epochs(datafile, chan_id, options) -% +% % ARGUMENTS: % datafile: File which contains the corresponding channel to % look for epochs. -% chan_id: Number or channeltype of the the channel. If -% multiple channels match the specified channeltype +% chan_id: Number or chantype of the the channel. If +% multiple channels match the specified chantype % only the first found channel will be used. -% options: +% options: % output_file: Name and path to the output file. Default is the % same as datafile but prepended with an 'e'. % overwrite: Overwrite output file if it already exists. Default % is 0. -% +% % OUTPUT: % sts: Return status of the function. If equals 1 no % error or warnings were produced. -% epochfile: File where the epochs have been saved to. +% epochfile: File where the epochs have been saved to. %__________________________________________________________________________ % PsPM 3.1 % (C) 2016 Tobias Moser (University of Zurich) @@ -90,7 +90,7 @@ if (isempty(stop) && ~isempty(start)) || (start(end) > stop(end)) stop = [stop; numel(logi)]; end; - + if (isempty(start) && ~isempty(stop)) || (start(1) > stop(1)) start = [1; start]; end; @@ -102,7 +102,7 @@ file_exist = exist(options.output_file, 'file'); write_ok = false; if file_exist - if options.overwrite + if options.overwrite write_ok = true; elseif ~options.dont_ask_overwrite if feature('ShowFigureWindows') diff --git a/src/pspm_cfg/pspm_cfg_data_convert.m b/src/pspm_cfg/pspm_cfg_data_convert.m index 428c79fea..62b8d11cd 100644 --- a/src/pspm_cfg/pspm_cfg_data_convert.m +++ b/src/pspm_cfg/pspm_cfg_data_convert.m @@ -117,7 +117,7 @@ visangle2sps.help = {['Takes paires of channels with gaze data in spherical',... ' coordinates (i.e. visual angle) and computes scanpath speed',... ' (i.e. scalar distance per second). Saves result into a ',... - 'new channel with channeltype ''sps'' (scanpath speed)']}; + 'new channel with chantype ''sps'' (scanpath speed)']}; %% Mode mode = cfg_choice; diff --git a/src/pspm_cfg/pspm_cfg_run_sf.m b/src/pspm_cfg/pspm_cfg_run_sf.m index 2da224e68..9638d4f9c 100644 --- a/src/pspm_cfg/pspm_cfg_run_sf.m +++ b/src/pspm_cfg/pspm_cfg_run_sf.m @@ -77,7 +77,11 @@ if ~isempty(job.fresp) options.fresp = job.fresp; end - +if ~isempty(job.missing) + if ischar(job.missing.missingdata) + options.missing = job.missing; + end +end options.dispwin = job.dispwin; options.dispsmallwin = job.dispsmallwin; diff --git a/src/pspm_cfg/pspm_cfg_sf.m b/src/pspm_cfg/pspm_cfg_sf.m index fa6411f93..f6b89ddf6 100644 --- a/src/pspm_cfg/pspm_cfg_sf.m +++ b/src/pspm_cfg/pspm_cfg_sf.m @@ -143,7 +143,7 @@ direction.labels = {'Unidirectional', 'Bidirectional'}; direction.values = {'uni', 'bi'}; direction.help = {['A unidirectional filter is applied twice in the forward direction. ' ... - 'A �bidirectional� filter is applied once in the forward direction and once in the ' ... + 'A "bidirectional" filter is applied once in the forward direction and once in the ' ... 'backward direction to correct the temporal shift due to filtering in forward direction.']}; filter_edit = cfg_branch; @@ -165,6 +165,8 @@ filter.values = {filter_def, filter_edit}; filter.help = {'Specify how you want filter the SCR data.'}; + + %% Epochs epochfile = cfg_files; epochfile.name = 'Epoch File'; @@ -191,11 +193,11 @@ mrk_chan.name = 'Marker Channel'; mrk_chan.tag = 'mrk_chan'; mrk_chan.strtype = 'i'; -mrk_chan.val = {0}; +mrk_chan.val = {1}; mrk_chan.num = [1 1]; mrk_chan.help = {['Indicate the marker channel. By default the first marker channel is ' ... 'assumed to contain the relevant markers.'], ['Markers are only used if you have ' ... - 'specified the time units as �markers�.']}; + 'specified the time units as "markers".']}; %% Timeunits @@ -228,8 +230,8 @@ timeunits.tag = 'timeunits'; timeunits.values = {seconds, samples, markers, whole}; timeunits.help = {['Indicate the time units on which the specification of the conditions will be based. ' ... - 'Time units can be specified in �seconds�, number of �markers�, or number of data �samples� . Time units ' ... - 'refer to the beginning of the data file and not to the beginning of the original recordings e. g. if ' ... + 'Time units can be specified in "seconds", number of "markers", or number of data "samples". Time units ' ... + 'refer to the beginning of the data file and not to the beginning of the original recordings e.g. if ' ... 'data were trimmed.']}; %% Channel nr @@ -292,6 +294,29 @@ fresp.help = {'Frequency of responses to model.'}; fresp.hidden = true; + + +missingepoch_none = cfg_const; +missingepoch_none.name = 'Not added'; +missingepoch_none.tag = 'missingdata'; +missingepoch_none.val = {0}; +missingepoch_none.help = {'Do not add missing epochs.'}; + +missingepoch_file = cfg_files; +missingepoch_file.name = 'Missing epoch file'; +missingepoch_file.tag = 'missingdata'; +missingepoch_file.num = [1 1]; +missingepoch_file.filter = '.*\.(mat|MAT)$'; +missingepoch_file.help = {['Missing (e.g. artefact) epochs in the data file, where ',... + 'data must always be specified in seconds.']}; + +missing = cfg_choice; +missing.name = 'Missing Epoch Settings'; +missing.tag = 'missing'; +missing.val = {missingepoch_none}; +missing.values = {missingepoch_none, missingepoch_file}; +missing.help = {'Specify whether you would like to include missing epochs.'}; + % Show figures dispwin = cfg_menu; dispwin.name = 'Display Progress Window'; @@ -314,7 +339,7 @@ sf = cfg_exbranch; sf.name = 'SF'; sf.tag = 'sf'; -sf.val = {datafile, modelfile, outdir, method, timeunits, filter, chan, overwrite, threshold, theta, fresp, dispwin, dispsmallwin}; +sf.val = {datafile, modelfile, outdir, method, timeunits, filter, chan, overwrite, threshold, missing, theta, fresp, dispwin, dispsmallwin}; sf.prog = @pspm_cfg_run_sf; sf.vout = @pspm_cfg_vout_sf; sf.help = {['This suite of models is designed for analysing spontaneous fluctuations (SF) in skin conductance ' ... diff --git a/src/pspm_dcm.m b/src/pspm_dcm.m index d668beb8e..732368585 100644 --- a/src/pspm_dcm.m +++ b/src/pspm_dcm.m @@ -30,7 +30,7 @@ % │ trial. If this is not the case, include `dummy` events with % │ negative onsets. % │ ▶︎ Optional -% ├───.missing: Allows to specify missing (e. g. artefact) epochs in the +% ├───.missing: Allows to specify missing (e.g. artefact) epochs in the % │ data file. See pspm_get_timing for epoch definition; specify % │ a cell array for multiple input files. This must always be % │ specified in SECONDS. diff --git a/src/pspm_dcm_inv.m b/src/pspm_dcm_inv.m index a3bdf7908..5ae46e9b1 100644 --- a/src/pspm_dcm_inv.m +++ b/src/pspm_dcm_inv.m @@ -34,6 +34,9 @@ % │ transformed back into raw data units. % ├───.flexevents: flexible events to adjust amplitude priors % ├────.fixevents: fixed events to adjust amplitude priors +% ├─.missing_data: missing epoch data, originally loaded as model.missing +% │ from pspm_dcm, but calculated into .missing_data (created +% │ in pspm_dcm and then transferred to pspm_dcm_inv. % └──.constrained: [optional] % constrained model for flexible responses which have % have fixed dispersion (0.3 s SD) but flexible latency @@ -64,7 +67,6 @@ % │ [numeric, default: 0.1, unit: second] % │ minimum dispersion (standard deviation) for flexible % │ responses. -% ├──────.missing: data points to be disregarded by inversion. % │ ▶︎ display % ├──────.dispwin: [bool, default as 1] % │ display progress window. @@ -126,9 +128,9 @@ try model.meanSCR; catch, model.meanSCR = 0; end % These parameters were set with default fallback values but will be % determined later by processing (same to pspm_dcm) -%try model.fixevents; catch, warning('model.fixevents not defined.'); end -%try model.flexevents; catch, warning('model.flexevents not defined.'); end -%try model.missing_data; catch, warning('model.missing_data not defined.'); end +% try model.fixevents; catch, warning('model.fixevents not defined.'); end +% try model.flexevents; catch, warning('model.flexevents not defined.'); end +% try model.missing_data; catch, warning('model.missing_data not defined.'); end % These parameters do not need to have a default value and will be % determined later (same to pspm_dcm) @@ -551,9 +553,10 @@ priors.SigmaX0 = zeros(7); priors.muX0 = zeros(7, 1); if trl == 1 - priors.muX0(1) = mean(y(1:3));%) - min(y); - priors.muX0(2) = mean(diff(y(1:3))); - priors.muX0(3) = diff(diff(y(1:3))); + y_non_nan = y(~isnan(y)); + priors.muX0(1) = mean(y_non_nan(1:3));%) - min(y); + priors.muX0(2) = mean(diff(y_non_nan(1:3))); + priors.muX0(3) = diff(diff(y_non_nan(1:3))); priors.muX0(7) = 0;%min(y); for n = [1:3 7] priors.SigmaX0(n, n) = 1e-2; diff --git a/src/pspm_guide.m b/src/pspm_guide.m index da1a514a7..8c95bba8f 100644 --- a/src/pspm_guide.m +++ b/src/pspm_guide.m @@ -1,7 +1,7 @@ function varargout = pspm_guide(varargin) -% ? Description +% ● Description % pspm_guide handles the main GUI for PsPM -% ? Developer's Guide +% ● Developer's Guide % * Template % function varargout = FunctionName(hObject, eventdata, handles, varargin) % varargout cell array for returning output args (see VARARGOUT); @@ -15,7 +15,7 @@ % contents{get(hObject,'Value')} returns selected item from the list % * Hint % popupmenu controls usually have a white background on Windows. -% ? Copyright +% ● Copyright % Introduced in PsPM 1.0 and terminated in PsPM 6.1. % Written by Dominik R Bach (Wellcome Centre for Human Neuroimaging) in 2008-2021 % Lastly updated in PsPM 6.1 by Teddy Chao (UCL) in 2022 diff --git a/src/pspm_interp1.m b/src/pspm_interp1.m new file mode 100644 index 000000000..86b1e5aae --- /dev/null +++ b/src/pspm_interp1.m @@ -0,0 +1,111 @@ +function Y = pspm_interp1(varargin) +% ● Description +% pspm_interp1 is a shared PsPM function for interpolating data with NaNs +% based on the reference of missing epochs and first order interpolation. +% ● Format +% Y = pspm_interp1(varargin) +% ● Arguments +% X: data that contains NaNs to be interpolated +% index_missing: index of missing epochs with the same size of X in binary +% values. 1 if NaNs, 0 if non-NaNs. +% Y: processed data +% ● History +% Introduced in PsPM 6.1 +% Written in 2023 by Teddy Chao (UCL) + +%% 1 Load inputs +switch nargin + case 1 + X = varargin{1}; + index_missing = zeros(size(X)); + case 2 + X = varargin{1}; + index_missing = varargin{2}; + otherwise + warning('ID:invalid_input','pspm_interp1 accepts up to two arguments'); +end +%% 2 Check inputs +switch sum(~isnan(X)) + case 0 + % if there are no non-nans, do not process any interpolation, give a + % warning and return + warning('ID:invalid_input',... + 'Input data contains only NaNs thus cannot be interpolated.') + Y = X; + return + case 1 + % if there are only 1 non-nan, do not process any interpolation, + % give a warning and explain the reason + warning('ID:invalid_input',... + 'Input data contains only 1 non-NaN thus cannot be interpolated.') + Y = X; + return + otherwise + % if there are less than 10^ non-nan, still perform interpolation, + % however give a warning and explain the reason + non_nan_percentage = sum(~isnan(X))/length(X); + if non_nan_percentage<0.1 + warning('ID:invalid_input',... + 'Input data contains less than 10% non-NaN. Interpolation can ',... + 'still be performed but results could be inaccurate.') + end +end +%% 3 find nan head and tail +X_nan_head = 0; +X_nan_tail = 0; +X_nan_head_range = []; +X_nan_tail_range = []; +X_nan_head_interp = []; +X_nan_tail_interp = []; +index_non_nan_full = 1:length(X); +index_non_nan_full = index_non_nan_full(~isnan(X)); +if index_non_nan_full(1) > 1 + X_nan_head = 1; + X_nan_head_range = 1:(index_non_nan_full(1)-1); +end +if index_non_nan_full(end) < length(X) + X_nan_tail = 1; + X_nan_tail_range = (index_non_nan_full(end)+1):length(X); +end +X_body = X(index_non_nan_full(1):index_non_nan_full(end)); +% processing body +index = 1:length(X_body); +index_nan = index(isnan(X_body) | index_missing(index_non_nan_full(1):index_non_nan_full(end))) ; +index_non_nan = 1 - index_nan; +if ~isempty(index_nan) + X_body_interp = interp1(index_non_nan,X_body(index_non_nan),index_nan); +else + X_body_interp = X_body; +end +% interpolate head +if X_nan_head + X_nan_head_interp = interp1(... + (1:length(X_body_interp))+length(X_nan_head_range),... + X_body_interp,... + X_nan_head_range,'linear','extrap'); +end +% interpolate tail +if X_nan_tail + X_nan_tail_interp = interp1(... + (1:length(X_body_interp))+length(X_nan_head_range),... + X_body_interp,... + X_nan_tail_range,'linear','extrap'); +end +if iscolumn(X_body_interp) + if ~iscolumn(X_nan_head_interp) + X_nan_head_interp = transpose(X_nan_head_interp); + end + if ~iscolumn(X_nan_tail_interp) + X_nan_tail_interp = transpose(X_nan_tail_interp); + end + Y = [X_nan_head_interp; X_body_interp; X_nan_tail_interp]; +else + if iscolumn(X_nan_head_interp) + X_nan_head_interp = transpose(X_nan_head_interp); + end + if iscolumn(X_nan_tail_interp) + X_nan_tail_interp = transpose(X_nan_tail_interp); + end + Y = [X_nan_head_interp, X_body_interp, X_nan_tail_interp]; +end +return diff --git a/src/pspm_load_data.m b/src/pspm_load_data.m index 1ffbc574f..c0d58822e 100644 --- a/src/pspm_load_data.m +++ b/src/pspm_load_data.m @@ -98,7 +98,6 @@ warning('ID:invalid_input', 'Too many inputs specified.'); return end - %% 3 Check fn % fn has to be a file or a struct switch class(fn) @@ -153,7 +152,6 @@ warning('ID:invalid_input', 'fn needs to be an existing file or a struct.'); return end - %% 4 Check channel switch class(channel) case 'double' @@ -190,7 +188,6 @@ otherwise warning('ID:invalid_input', 'Unknown channel option.'); end - %% 5 Check infos if isstruct(channel) infos = channel.infos; @@ -217,8 +214,7 @@ warning('ID:invalid_data_structure', 'Input data does not have sufficient infos'); return end - -%% 6 Check data +%% 6 Load data if isstruct(channel) data = channel.data; else @@ -232,18 +228,24 @@ clear loaded_data end end -% initialise error flags +%% 7 Check data +% 7.1 initialise error flags -- vflag = zeros(numel(data), 1); % records data structure, valid if 0 wflag = zeros(numel(data), 1); % records whether data is out of range, valid if 0 nflag = zeros(numel(data), 1); zflag = zeros(numel(data), 1); % records whether data is empty % loop through channels for k = 1:numel(data) - % check header + % 7.2 Check header -- if ~isfield(data{k}, 'header') vflag(k) = 1; else - if (~isfield(data{k}.header, 'channeltype') && ~isfield(data{k}.header, 'chantype')) || ... + % 7.2.1 Convert header channeltype into chantype if there are -- + if isfield(data{k}.header, 'channeltype') + data{k}.header.chantype = data{k}.header.channeltype; + data{k}.header = rmfield(data{k}.header, 'channeltype'); + end + if ~isfield(data{k}.header, 'chantype') || ... ~isfield(data{k}.header, 'sr') || ... ~isfield(data{k}.header, 'units') vflag(k) = 1; @@ -259,7 +261,7 @@ end end end - % check data + % 7.3 Check data -- if vflag(k)==0 && nflag(k)==0 && flag_infos==0 % required information is available and valid in header and infos if ~isfield(data{k}, 'data') @@ -307,18 +309,6 @@ warning('ID:missing_data', 'Channel %01.0f is empty.', find(zflag,1)); % if there is empty data, give a warning but do not suspend end - - -%% 7 Autofill information in header -% some other optional fields which can be autofilled with default values -% should be added here. -for k = 1:numel(data) - if isfield(data{k}.header, 'channeltype') - data{k}.header.chantype = data{k}.header.channeltype; - data{k}.header = rmfield( data{k}.header, 'channeltype'); - end -end - %% 8 Analyse file structure filestruct.numofwavechan = 0; filestruct.numofeventchan = 0; @@ -339,7 +329,6 @@ elseif numel(filestruct.posofmarker) > 1 filestruct.posofmarker = filestruct.posofmarker(1); % first marker channel end - %% 9 Return channels, or save file if isstruct(channel) infos = channel.infos; @@ -391,7 +380,6 @@ else filestruct.posofchannels = []; end - sts = 1; return @@ -438,7 +426,6 @@ combined_channels = combined_channels & pupil_channels; besteye_channels = besteye_channels & pupil_channels & ~preprocessed_channels; % best eye will not select preprocessed eyes - if any(combined_channels) flag = combined_channels; elseif any(preprocessed_channels) && ~prefer_unprocessed @@ -450,7 +437,6 @@ flag = besteye_channels; end - function flag = get_chans_to_load_for_sps(data, best_eye) % 16-06-21 This is a tempory patch for loading sps data, copied from % pupil data @@ -494,7 +480,6 @@ preprocessed_channels = preprocessed_channels & sps_channels; combined_channels = combined_channels & sps_channels; besteye_channels = besteye_channels & sps_channels; - if any(combined_channels) flag = combined_channels; elseif any(preprocessed_channels) @@ -504,4 +489,4 @@ end else flag = besteye_channels; -end +end \ No newline at end of file diff --git a/src/pspm_options.m b/src/pspm_options.m index 09ee2086d..4b982bd52 100644 --- a/src/pspm_options.m +++ b/src/pspm_options.m @@ -411,8 +411,10 @@ options = autofill(options,'dispwin', 1, 0 ); options = autofill(options,'fresp', 0.5, '>=', 0 ); options = autofill(options,'marker_chan_num', 1, '*Int*Char' ); + options = autofill(options,'missingthresh', 2, '>', 0 ); options = autofill(options,'overwrite', 1, 0 ); options = autofill(options,'threshold', 0.1, '>', 0 ); + options = autofill(options,'missingthresh', 2, '>', 0 ); options = autofill(options,'theta', [0.923581, ... 3.921034, ... 2.159389, ... @@ -424,7 +426,6 @@ options = autofill(options,'dispwin', 1, 0 ); options = autofill(options,'dispsmallwin', 0, 1 ); options = autofill(options,'fresp', 0.5, '>', 0 ); - options = autofill(options,'missingthresh', 2, '>', 0 ); options = autofill(options,'threshold', 0.1, '>', 0 ); options = autofill(options,'theta', [0.923581, ... 3.921034, ... @@ -579,7 +580,7 @@ if contains(optional_value, '*Int') flag_is_allowed_value = flag_is_allowed_value || ... all([isnumeric(options.(field_name)), ... - options.(field_name)>0, ... + options.(field_name)>=0, ... mod(options.(field_name), 1)==0]); end if contains(optional_value, '*Struct') diff --git a/src/pspm_overwrite.m b/src/pspm_overwrite.m index f6c446550..f8aa1e642 100644 --- a/src/pspm_overwrite.m +++ b/src/pspm_overwrite.m @@ -96,5 +96,7 @@ case 2 varargout{1} = sts; varargout{2} = overwrite_final; + otherwise + varargout{1} = overwrite_final; end return diff --git a/src/pspm_prepdata.m b/src/pspm_prepdata.m index a2838a471..d9b2a4f0b 100644 --- a/src/pspm_prepdata.m +++ b/src/pspm_prepdata.m @@ -197,63 +197,3 @@ varargout{3} = newsr; end return -function data_interp = pspm_interp1(data) -% find nan head and tail -data_nan_head = 0; -data_nan_tail = 0; -data_nan_head_range = []; -data_nan_tail_range = []; -data_nan_head_interp = []; -data_nan_tail_interp = []; -index_non_nan_full = 1:length(data); -index_non_nan_full = index_non_nan_full(~isnan(data)); -if index_non_nan_full(1) > 1 - data_nan_head = 1; - data_nan_head_range = 1:(index_non_nan_full(1)-1); -end -if index_non_nan_full(end) < length(data) - data_nan_tail = 1; - data_nan_tail_range = (index_non_nan_full(end)+1):length(data); -end -data_body = data(index_non_nan_full(1):index_non_nan_full(end)); -% processing body -index = 1:length(data_body); -index_nan = index(isnan(data_body)); -index_non_nan = index(~isnan(data_body)); -if ~isempty(index_nan) - data_body_interp = interp1(index_non_nan,data_body(index_non_nan),index_nan); -else - data_body_interp = data_body; -end -% interpolate head -if data_nan_head - data_nan_head_interp = interp1(... - (1:length(data_body_interp))+length(data_nan_head_range),... - data_body_interp,... - data_nan_head_range,'linear','extrap'); -end -% interpolate tail -if data_nan_tail - data_nan_tail_interp = interp1(... - (1:length(data_body_interp))+length(data_nan_head_range),... - data_body_interp,... - data_nan_tail_range,'linear','extrap'); -end -if iscolumn(data_body_interp) - if ~iscolumn(data_nan_head_interp) - data_nan_head_interp = transpose(data_nan_head_interp); - end - if ~iscolumn(data_nan_tail_interp) - data_nan_tail_interp = transpose(data_nan_tail_interp); - end - data_interp = [data_nan_head_interp; data_body_interp; data_nan_tail_interp]; -else - if iscolumn(data_nan_head_interp) - data_nan_head_interp = transpose(data_nan_head_interp); - end - if iscolumn(data_nan_tail_interp) - data_nan_tail_interp = transpose(data_nan_tail_interp); - end - data_interp = [data_nan_head_interp, data_body_interp, data_nan_tail_interp]; -end -return diff --git a/src/pspm_resp_pp.m b/src/pspm_resp_pp.m index ed3f954cc..96fbb561d 100644 --- a/src/pspm_resp_pp.m +++ b/src/pspm_resp_pp.m @@ -85,7 +85,7 @@ filt.hporder = 1; filt.direction = 'bi'; filt.down = 10; -[sts, newresp, newsr] = pspm_prepdata(resp - mean(resp), filt); +[sts, newresp, newsr] = pspm_prepdata(resp - mean(resp,"omitnan"), filt); % Median filter newresp = medfilt1(newresp, ceil(newsr) + 1); %% detect breathing cycles diff --git a/src/pspm_sf.m b/src/pspm_sf.m index 93fdaadd5..48c7d426f 100644 --- a/src/pspm_sf.m +++ b/src/pspm_sf.m @@ -25,21 +25,33 @@ % │ as 'dcm' % │ [cell_array] a cell array of methods mentioned above. % ├─────────.filter: filter settings; modality specific default -% └────────.channel: channel number; default: first SCR channel +% ├────────.missing: [string/cell_array] [default: no missing values] +% │ Allows to specify missing (e.g. artefact) epochs in the +% │ data file. See pspm_get_timing for epoch definition; specify +% │ a cell array for multiple input files. This must always be +% │ specified in SECONDS. +% └────────.channel: [integer] [default: first SCR channel] +% channel number. % ┌─────────options -% ├──────.overwrite: [logical] (0 or 1) +% ├──────.overwrite: [logical] [default: determined by pspm_overwrite] % │ Define whether to overwrite existing output files or not. -% │ Default value: determined by pspm_overwrite. -% ├.marker_chan_num: marker channel number +% ├.marker_chan_num: [integer] +% │ marker channel number % │ if undefined or 0, first marker channel is used. -% │ additional options for individual methods: +% │ * Additional options for individual methods: % │ dcm related options -% ├──────.threshold: threshold for SN detection (default 0.1 mcS) -% ├──────────.theta: a (1 x 5) vector of theta values for f_SF -% │ (default: read from pspm_sf_theta) -% ├──────────.fresp: frequency of responses to model (default 0.5 Hz) -% ├────────.dispwin: display progress window (default 1) -% └───.dispsmallwin: display intermediate windows (default 0); +% ├──────.threshold: [numeric] [default: 0.1] [unit: mcS] +% │ threshold for SN detection (default 0.1 mcS) +% ├──────────.theta: [vector] [default: read from pspm_sf_theta] +% │ A (1 x 5) vector of theta values for f_SF. +% ├──────────.fresp: [numeric] [unit: Hz] [default: 0.5] +% │ frequency of responses to model. +% ├────────.dispwin: [logical] [default: 1] +% │ display progress window. +% ├───.dispsmallwin: [logical] [default: 0] +% │ display intermediate windows. +% └──.missingthresh: [numeric] [default: 2] [unit: second] +% threshold value for controlling missing epochs. % ● References % 1.[DCM for SF] % Bach DR, Daunizeau J, Kuelzow N, Friston KJ, Dolan RJ (2010). Dynamic @@ -65,7 +77,7 @@ outfile = []; sts = -1; %% 2 Check input -% 2.1 Check missing input +% 2.1 Check missing input -- if nargin<1 warning('ID:invalid_input', 'Nothing to do.'); return; elseif nargin<2 @@ -80,7 +92,7 @@ elseif ~isfield(model, 'timing') && ~strcmpi(model.timeunits, 'file') warning('ID:invalid_input', 'No epochs specified.'); return; end -% 2.2 Check faulty input +% 2.2 Check faulty input -- if ~ischar(model.datafile) && ~iscell(model.datafile) warning('ID:invalid_input', 'Input data must be a cell or string.'); return; elseif ~ischar(model.modelfile) && ~iscell(model.modelfile) @@ -103,7 +115,7 @@ if ischar(model.modelfile) model.modelfile = {model.modelfile}; end -% 2.4 Check number of files +% 2.4 Check number of files -- if ~strcmpi(model.timeunits, 'whole') && numel(model.datafile) ~= numel(model.timing) warning('ID:number_of_elements_dont_match',... 'Number of data files and epoch definitions does not match.'); return; @@ -111,7 +123,7 @@ warning('ID:number_of_elements_dont_match',... 'Number of data files and model files does not match.'); return; end -% 2.5 check methods +% 2.5 check methods -- if ~isfield(model, 'method') model.method = {'dcm'}; elseif ischar(model.method) @@ -150,39 +162,50 @@ end end end -% 2.6 Check timing +% 2.6 Check timing -- if strcmpi(model.timeunits, 'whole') epochs = repmat({[1 1]}, numel(model.datafile), 1); else - for iSn = 1:numel(model.datafile) - [sts_get_timing, epochs{iSn}] = pspm_get_timing('epochs', model.timing{iSn}, model.timeunits); + for iFile = 1:numel(model.datafile) + [sts_get_timing, epochs{iFile}] = pspm_get_timing('epochs', model.timing{iFile}, model.timeunits); if sts_get_timing == -1 warning('ID:invalid_input', 'Call of pspm_get_timing failed.'); return; end end end -% 2.7 Check filter +% 2.7 Check filter -- if ~isfield(model, 'filter') model.filter = settings.dcm{2}.filter; elseif ~isfield(model.filter, 'down') || ~isnumeric(model.filter.down) warning('ID:invalid_input', 'Filter structure needs a numeric ''down'' field.'); return; end -% 2.8 Set options +% 2.8 Set options -- try model.channel; catch, model.channel = 'scr'; end options = pspm_options(options, 'sf'); if options.invalid return end +% 2.9 Set missing epochs -- +if ~isfield(model, 'missing') + model.missing = cell(numel(model.datafile), 1); +elseif ischar(model.missing) || isnumeric(model.missing) + model.missing = {model.missing}; +elseif ~iscell(model.missing) + warning('ID:invalid_input',... + 'Missing values must be a filename, matrix, or cell array of these.'); + return +end %% 3 Get data +missing = cell(size(model.missing)); for iFile = 1:numel(model.datafile) - % 3.1 User output + % 3.1 User output -- fprintf('SF analysis: %s ...', model.datafile{iFile}); - % 3.2 Check whether model file exists + % 3.2 Check whether model file exists -- if ~pspm_overwrite(model.modelfile, options) return end - % 3.3 get and filter data + % 3.3 get and filter data -- [sts_load_data, ~, data] = pspm_load_data(model.datafile{iFile}, model.channel); if sts_load_data < 0, return; end Y{1} = data{1}.data; sr(1) = data{1}.header.sr; @@ -192,40 +215,42 @@ warning('ID:invalid_input', 'Call of pspm_prepdata failed.'); return; end - % 3.4 Check data units + % 3.4 Check data units -- if ~strcmpi(data{1}.header.units, 'uS') && any(strcmpi('dcm', method)) fprintf(['\nYour data units are stored as %s, ',... 'and the method will apply an amplitude threshold in uS. ',... 'Please check your results.\n'], ... data{1}.header.units); end - % 3.5 Get marker data - if any(strcmp(model.timeunits, {'marker', 'markers'})) - if options.marker_chan_num - [nsts, ~, ndata] = pspm_load_data(model.datafile, options.marker_chan_num); - if nsts == -1 - warning('ID:invalid_input', 'Could not load data'); - return; - end - if ~strcmp(ndata{1}.header.chantype, 'marker') - warning('ID:invalid_option', ... - ['Channel %i is no marker channel. ',... - 'The first marker channel in the file is used instead'],... - options.marker_chan_num); - [nsts, ~, ~] = pspm_load_data(model.datafile, 'marker'); - if nsts == -1 - warning('ID:invalid_input', 'Could not load data'); - return; + % 3.5 Get missing epochs -- + % 3.5.1 Load missing epochs -- + if ~isempty(model.missing{iFile}) + [~, missing{iFile}] = pspm_get_timing('epochs', model.missing{iFile}, 'seconds'); + % 3.5.2 sort missing epochs -- + if size(missing{iFile}, 1) > 0 + [~, sortindx] = sort(missing{iFile}(:, 1)); + missing{iFile} = missing{iFile}(sortindx,:); + % check for overlap and merge + for k = 2:size(missing{iFile}, 1) + if missing{iFile}(k, 1) <= missing{iFile}(k - 1, 2) + missing{iFile}(k, 1) = missing{iFile}(k - 1, 1); + missing{iFile}(k - 1, :) = []; end end + end + else + missing{iFile} = []; + end + % 3.6 Get marker data -- + if any(strcmp(model.timeunits, {'marker', 'markers'})) + if options.marker_chan_num + [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, options.marker_chan_num); + if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end else - [nsts, ~, ~] = pspm_load_data(model.datafile, 'marker'); - if nsts == -1 - warning('ID:invalid_input', 'Could not load data'); - return; - end + [nsts, ~, ndata] = pspm_load_data(model.datafile{iFile}, 'marker'); + if nsts == -1; warning('ID:invalid_input', 'Could not load data'); return; end end - events = data{1}.data; + events = ndata{1}.data; end for iEpoch = 1:size(epochs{iFile}, 1) if iEpoch > 1, fprintf('\n\t\t\t'); end @@ -238,7 +263,7 @@ case 'samples' win = round(epochs{iFile}(iEpoch, :) * sr(datatype(k)) / sr(1)); case 'markers' - win = round(events(epochs{1}(iEpoch, :)) * sr(datatype(k))); + win = round(events(epochs{iFile}(iEpoch, :)) * sr(datatype(k))); case 'whole' win = [1 numel(Y{datatype(k)})]; end @@ -249,7 +274,7 @@ win(1) = max(win(1), 1); win(2) = min(win(2), numel(Y{datatype(k)})); end - % 3.6 collect information + % 3.6.1 collect information -- sf.model{k}(iEpoch).modeltype = method{k}; sf.model{k}(iEpoch).boundaries = squeeze(epochs{iFile}(iEpoch, :)); sf.model{k}(iEpoch).timeunits = model.timeunits; @@ -258,8 +283,18 @@ % escr = Y{datatype(k)}(win(1):win(end)); sf.model{k}(iEpoch).data = escr; - % 3.7 do the analysis and collect results - invrs = fhandle{k}(escr, sr(datatype(k)), options); + if any(missing{iFile}) + model.missing_data = zeros(size(escr)); + missing_index = pspm_time2index(missing, sr(datatype(k))); + model.missing_data((missing_index{iFile}(:,1)+1):(missing_index{iFile}(:,2)+1)) = 1; + end + % 3.6.2 do the analysis and collect results -- + if any(missing{iFile}) + model_analysis = struct('scr', escr, 'sr', sr(datatype(k)), 'missing_data', model.missing_data); + else + model_analysis = struct('scr', escr, 'sr', sr(datatype(k))); + end + invrs = fhandle{k}(model_analysis, options); if any(strcmpi(method{k}, {'dcm', 'mp'})) sf.model{k}(iEpoch).inv = invrs; sf.stats(iEpoch, k) = invrs.f; diff --git a/src/pspm_sf_auc.m b/src/pspm_sf_auc.m index bc2b2a0de..96bdc378a 100644 --- a/src/pspm_sf_auc.m +++ b/src/pspm_sf_auc.m @@ -1,4 +1,4 @@ -function varargout = pspm_sf_auc(scr, sr, options) +function varargout = pspm_sf_auc(model, options) % ● Description % pspm_sf_auc returns the integral/area under the curve of an SCR time series % ● Format @@ -28,6 +28,10 @@ if nargin < 1 warning('No data specified'); return; end; +try model.scr; catch, warning('Input data is not defined.'); return; end +try model.sr; catch, warning('Sample rate is not defined.'); return; end +scr = model.scr; +sr = model.sr; scr = scr - min(scr); auc = mean(scr); sts = 1; diff --git a/src/pspm_sf_dcm.m b/src/pspm_sf_dcm.m index 494a854f1..2df113ac7 100644 --- a/src/pspm_sf_dcm.m +++ b/src/pspm_sf_dcm.m @@ -1,37 +1,37 @@ -function varargout = pspm_sf_dcm(scr, sr, options) +function varargout = pspm_sf_dcm(model, options) % ● Description % pspm_sf_dcm does dynamic causal modelling for SF of the skin conductance % uses f_SF and g_Id % the input data is assumed to be in mcS, and sampling rate in Hz % ● Format -% function out = pspm_sf_dcm(scr, sr, options) -% ● Output -% out: output -% .n: number of responses above threshold -% .f: frequency of responses above threshold in Hz -% .ma: mean amplitude of responses above threshold -% .t: timing of (all) responses -% .a: amplitude of (all) responses -% .theta: parameters used for f_SF -% .threshold: threshold -% .if: initial frequency for f_SF -% .yhat: fitted time series -% .model: information about the DCM inversion +% function out = pspm_sf_dcm(model, options) % ● Arguments -% scr: skin conductance epoch (maximum size depends on computing -% power, a sensible size is 60 s at 10 Hz) -% sr: sampling rate in Hz -% options: options structure -% .threshold: threshold for SN detection (default 0.1 mcS) -% .theta: a (1 x 5) vector of theta values for f_SF -% (default: read from pspm_sf_theta) -% .fresp: frequency of responses to model (default 0.5 Hz) -% .dispwin: display progress window (default 1) -% .dispsmallwin: -% display intermediate windows (default 0); -% .missing: index of missing values to ignore -% .missingthresh: -% threshold value for controlling missing epochs (default 2s). +% ┌──────model +% │ ▶︎ Mandatory +% ├───────.scr: skin conductance epoch (maximum size depends on computing +% │ power, a sensible size is 60 s at 10 Hz) +% ├────────.sr: [numeric] [unit: Hz] +% │ sampling rate. +% │ ▶︎ Optional +% └.missing_data:missing epoch data, originally loaded as model.missing +% from pspm_sf, but calculated into .missing_data (created +% in pspm_sf and then transferred to pspm_sf_dcm. +% +% ┌────options +% ├─.threshold: [numeric] [default: 0.1] [unit: mcS] +% │ threshold for SN detection (default 0.1 mcS) +% ├─────.theta: [vector] [default: read from pspm_sf_theta] +% │ a (1 x 5) vector of theta values for f_SF +% ├─────.fresp: [numeric] [unit: Hz] [default: 0.5] +% │ frequency of responses to model +% ├───.dispwin: [logical] [default: 1] +% │ display progress window. +% ├.dispsmallwin:[logical] [default: 0] +% │ display intermediate windows. +% └.missingthresh: +% [numeric] [default: 2] [unit: second] +% threshold value for controlling missing epochs, +% which is originally inherited from SF % ● References % Bach DR, Daunizeau J, Kuelzow N, Friston KJ, & Dolan RJ (2011). Dynamic % causal modelling of spontaneous fluctuations in skin conductance. @@ -48,18 +48,22 @@ sts = -1; tstart = tic; %% 2 Check input arguments -if nargin < 2 || ~isnumeric(sr) || numel(sr) > 1 +% 2.1 set model --- +try model.scr; catch, warning('Input data is not defined.'); return; end +try model.sr; catch, warning('Sample rate is not defined.'); return; end +% 2.2 Validate parameters --- +if ~isnumeric(model.sr) || numel(model.sr) > 1 errmsg = sprintf('No valid sample rate given.'); -elseif (sr < 1) || (sr > 1e5) +elseif (model.sr < 1) || (model.sr > 1e5) errmsg = sprintf('Sample rate out of range.'); -elseif exist('osr', 'var') && osr ~= sr +elseif exist('osr', 'var') && osr ~= model.sr errmsg = sprintf('Sample rate of theta file is different from sample rate of data.'); -elseif nargin < 1 || ~isnumeric(scr) +elseif nargin < 1 || ~isnumeric(model.scr) errmsg = 'No data.'; -elseif ~any(size(scr) == 1) +elseif ~any(size(model.scr) == 1) errmsg = 'Input SCR is not a vector'; else - scr = scr(:); + model.scr = model.scr(:); end if exist('errmsg', 'var') == 1 warning(errmsg); @@ -71,8 +75,8 @@ if options.invalid return end -options.DisplayWin = options.dispwin; -options.GnFigs = options.dispsmallwin; +% options.DisplayWin = options.dispwin; +% options.GnFigs = options.dispsmallwin; fresp = options.fresp; threshold = options.threshold; try @@ -93,35 +97,36 @@ priors.a_alpha = Inf; priors.b_alpha = 0; % 4.2 initialise priors in correct dimensions -priors.iQy = cell(numel(scr), 1); -priors.iQx = cell(numel(scr), 1); -for k = 1:numel(scr) % default priors on noise covariance +priors.iQy = cell(numel(model.scr), 1); +priors.iQx = cell(numel(model.scr), 1); +for k = 1:numel(model.scr) % default priors on noise covariance priors.iQy{k} = 1; priors.iQx{k} = eye(dim.n); end options.inG.ind = 1; -options.inF.dt = 1/sr; +options.inF.dt = 1/model.sr; % 4.3 prepare data -y = scr; +y = model.scr; y = y - min(y); % 4.4 determine initial conditions -x0 = y(1:3); +y_non_nan = y(~isnan(y)); +x0 = y_non_nan(1:3); X0(1, 1) = mean(x0); X0(2, 1) = mean(diff(x0)); X0(3, 1) = diff(diff(x0)); priors.muX0 = X0; -nresp = floor(fresp * numel(y)/sr) + 1; +nresp = floor(fresp * numel(y)/model.sr) + 1; u = []; -u(1, :) = (1:numel(y))/sr; +u(1, :) = (1:numel(y))/model.sr; u(2, :) = nresp; -priors.muTheta = theta(1:3)'; +priors.muTheta = transpose(theta(1:3)); priors.muTheta(4:2:(2 * nresp + 3)) = 1/fresp * (0:(nresp-1)); priors.muTheta(5:2:(2 * nresp + 4)) = -10; dim.n_theta = numel(priors.muTheta); % nb of evolution parameters priors.SigmaTheta = zeros(dim.n_theta); for k = (4:2:(2 * nresp + 3)), priors.SigmaTheta(k, k) = 1e-2;end for k = (5:2:(2 * nresp + 4)), priors.SigmaTheta(k, k) = 1e2; end -priors.muPhi = phi'; +priors.muPhi = transpose(phi); priors.SigmaPhi = zeros(dim.n_phi); priors.SigmaX0 = 1e-8*eye(dim.n); options.priors = priors; @@ -129,8 +134,8 @@ c = clock; fprintf(['\n\nEstimating model parameters for f_SF ... \t%02.0f:%02.0f:%02.0f', ... '\n=========================================================\n'], c(4:6)); -if isfield(options, 'missing') - ymissing = options.missing; +if isfield(model, 'missing_data') + ymissing = model.missing_data; else ymissing = isnan(y); end @@ -143,22 +148,30 @@ end miss_epoch = [ymissing_start(:),ymissing_end(:)]; flag_missing_too_long = 0; -if any(diff(miss_epoch, 1, 2)/sr > options.missingthresh) - warning_message = ['Imported data includes too long miss epoches (over ',... - num2str(options.missingthresh), 's), thus estimation has been skipped.']; +if any(diff(miss_epoch, 1, 2)/model.sr > 0) + if any(diff(miss_epoch, 1, 2)/model.sr > options.missingthresh) + warning_message = ['Imported data includes too long miss epoches (over ',... + num2str(options.missingthresh), 's), thus estimation has been skipped.']; + flag_missing_too_long = 1; + else + warning_message = ['Imported data includes miss epoches (over ',... + num2str(options.missingthresh), 's), but the trial has been allowed. ',... + 'Please adjust options.missingthresh to skip if you wish.']; + end warning('ID:missing_data', warning_message); - flag_missing_too_long = 1; end options.isYout = ymissing(:)'; +% 4.6 interpolate data body to fill NaNs +y_interpolated = pspm_interp1(y, ymissing); %% 5 Extract parameters if ~flag_missing_too_long - [posterior, output] = VBA_NLStateSpaceModel(y(:)',u,f_fname,g_fname,dim,options); + [posterior, output] = VBA_NLStateSpaceModel(y_interpolated(:)',u,f_fname,g_fname,dim,options); for i = 1:length(output) output(i).options = rmfield(output(i).options, 'hf'); end t = posterior.muTheta(4:2:end); a = exp(posterior.muTheta(5:2:end) - theta(5)); % rescale - ex = find(t < -2 | t > (numel(scr)/sr - 1)); % find SA responses the SCR peak of which is outside episode + ex = find(t < -2 | t > (numel(model.scr)/model.sr - 1)); % find SA responses the SCR peak of which is outside episode t(ex) = []; a(ex) = []; end @@ -167,7 +180,7 @@ out.t = t - theta(4); % subtract conduction delay out.a = a; out.n = numel(find(a > threshold)); - out.f = out.n/(numel(scr)/sr); + out.f = out.n/(numel(model.scr)/model.sr); out.ma = mean(a(a > threshold)); out.theta = theta; out.if = fresp; @@ -176,7 +189,7 @@ out.model.posterior = posterior; out.model.output = output; out.model.u = u; - out.model.y = y(:)'; + out.model.y = y_interpolated(:)'; out.time = toc(tstart); else out.t = NaN; diff --git a/src/pspm_sf_mp.m b/src/pspm_sf_mp.m index 944ae5379..6a9f8f0b2 100644 --- a/src/pspm_sf_mp.m +++ b/src/pspm_sf_mp.m @@ -1,4 +1,4 @@ -function varargout = pspm_sf_mp(scr, sr, options) +function varargout = pspm_sf_mp(model, options) % ● Description % pspm_sf_mp does the inversion of a DCM for SF of the skin conductance, using % a matching pursuit algorithm, and f_SF for the forward model @@ -47,9 +47,14 @@ sts = -1; tstart = tic; +try model.scr; catch, warning('Input data is not defined.'); return; end +try model.sr; catch, warning('Sample rate is not defined.'); return; end +scr = model.scr; +sr = model.sr; + % check input arguments % ------------------------------------------------------------------------ -if nargin < 2 || ~isnumeric(sr) || numel(sr) > 1 +if ~isnumeric(sr) || numel(sr) > 1 errmsg = sprintf('No valid sample rate given.'); elseif (sr < 1) || (sr > 1e5) errmsg = sprintf('Sample rate out of range.'); @@ -73,6 +78,7 @@ S.dt = 1/sr; % sampling interval of the data S.n = numel(scr); % n: number of samples in the data segment S.sfduration = 30; % duration of the modelled SF +S.fresp = options.fresp; S.sftail = 10; % model tail of SF in previous seconds S.ntail = S.sftail/S.dt; % number of samples in SF tail to model S.nsf = S.sfduration/S.dt; % number of samples in modelled SF @@ -82,6 +88,7 @@ S.theta = pspm_sf_theta; % get SF CRF S.maxres = 0.001 * S.n; % residual threshold per sample S.options = options; +S.threshold = options.threshold; % generate over-complete dictionary D.D % ------------------------------------------------------------------------- @@ -157,12 +164,15 @@ S.nosf = 0; k = 1; S.Yres = y; -ind = []; % initialise diagnostics S.diagnostics.neg = 0; % stop because of zero or negative amplitude S.diagnostics.num = 0; % number of iterations S.diagnostics.error = NaN; % error +a = []; +asf = []; +ind = []; + % run algorithm --- while S.cont % remove used atoms from dictionary @@ -252,9 +262,9 @@ sts = 1; switch nargout case 1 - varargout{1} = mp; + varargout{1} = out; case 2 varargout{1} = sts; - varargout{2} = mp; + varargout{2} = out; end return diff --git a/src/pspm_sf_scl.m b/src/pspm_sf_scl.m index e82c7c5d4..f8e1d9ee7 100644 --- a/src/pspm_sf_scl.m +++ b/src/pspm_sf_scl.m @@ -1,4 +1,4 @@ -function varargout = pspm_sf_scl(scr, sr, options) +function varargout = pspm_sf_scl(model, options) % ● Description % pspm_sf_scl returns the mean skin conductance level for an epoch % ● Format @@ -25,6 +25,8 @@ if nargin < 1 warning('No data specified'); return; end; +try model.scr; catch, warning('Input data is not defined.'); return; end +try model.sr; catch, warning('Sample rate is not defined.'); return; end scl = mean(scr); sts = 1; switch nargout diff --git a/src/pspm_split_sessions.m b/src/pspm_split_sessions.m index b0c255ffc..4c9a756aa 100644 --- a/src/pspm_split_sessions.m +++ b/src/pspm_split_sessions.m @@ -205,7 +205,7 @@ newepochfile{sn} = fullfile(p_epochs, sprintf('%s_sn%02.0f%s', f_epochs, sn, ex_epochs)); end % 2.4.2 Split data - trimoptions = struct('drop_offset_markers', 1); + trimoptions = struct('drop_offset_markers', 1, 'marker_chan_num', markerchannel); newdata = pspm_trim(struct('data', {indata}, 'infos', ininfos), ... options.prefix, suffix(sn), trimpoint(sn, 1:2), trimoptions); options.overwrite = pspm_overwrite(newdatafile{sn}, options); @@ -213,7 +213,7 @@ pspm_load_data(newdatafile{sn}, newdata); % 2.4.5 Split Epochs if ~isempty(options.missing) - dummydata{1,1}.header = struct('channeltype', 'custom', ... + dummydata{1,1}.header = struct('chantype', 'custom', ... 'sr', missingsr, ... 'units', 'unknown'); dummydata{1,1}.data = dp_epochs;