diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..34d1f332 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,28 @@ +# Changelog +TAPAS toolbox + + + +## [2.7.3] 2017-08-17 + +## Added +- Added condhalluc_obs and condhalluc_obs2 models. + +## Changed +- Introduced kappa1 in all binary HGF models. + +## [2.7.2] 2017-08-17 +## Added +- New PhysIO CHANGELOG.md specific file. + +## Changed +- Specified in PhysIO/CHANGELOG.md. + +## [2.7.1] 2017-08-17 + +### Added +- Added a CHANGELOG.md file + +### Changed +- Now the README.txt file is in markdown format. +- The documentation is integrated with the wiki in github. diff --git a/HGF/README.md b/HGF/README.md index 77c9366c..8e227531 100644 --- a/HGF/README.md +++ b/HGF/README.md @@ -46,6 +46,10 @@ hgf_demo.pdf. ## Release notes +### v5.1 +- Added condhalluc_obs and condhalluc_obs2 models +- Introduced kappa1 in all binary HGF models + ### v5.0 - Ported interactive demo to Matlab LiveScript - Various additional small improvements diff --git a/HGF/hgf_demo.m b/HGF/hgf_demo.m index a007f230..970e1049 100644 --- a/HGF/hgf_demo.m +++ b/HGF/hgf_demo.m @@ -78,7 +78,7 @@ %% sim = tapas_simModel(u,... 'tapas_hgf_binary',... - [NaN 0 1 NaN 1 1 NaN 0 0 NaN 1 NaN -2.5 -6],... + [NaN 0 1 NaN 1 1 NaN 0 0 1 1 NaN -2.5 -6],... 'tapas_unitsq_sgm',... 5); %% diff --git a/HGF/hgf_demo.mlx b/HGF/hgf_demo.mlx index 257ad307..e3f9cbd6 100644 Binary files a/HGF/hgf_demo.mlx and b/HGF/hgf_demo.mlx differ diff --git a/HGF/tapas_condhalluc_obs.m b/HGF/tapas_condhalluc_obs.m new file mode 100644 index 00000000..cb13301b --- /dev/null +++ b/HGF/tapas_condhalluc_obs.m @@ -0,0 +1,49 @@ +function [logp, yhat, res] = tapas_condhalluc_obs(r, infStates, ptrans) +% Calculates the log-probability of response y=1 under the unit-square sigmoid model +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2015-2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Transform beta to its native space +be = exp(ptrans(1)); + +% Initialize returned log-probabilities as NaNs so that NaN is +% returned for all irregualar trials +n = size(infStates,1); +logp = NaN(n,1); +yhat = NaN(n,1); +res = NaN(n,1); + +% Check input format +if size(r.u,2) ~= 2 + error('tapas:hgf:CondHalluc:InputsIncompatible', 'Inputs incompatible with condhalluc_obs observation model. See tapas_condhalluc_obs_config.m.') +end + +% Get true-positive rate corresponding to stimuli +tp = r.u(:,2); + +% Weed irregular trials out +mu1hat = infStates(:,1,1); +mu1hat(r.irr) = []; +y = r.y(:,1); +y(r.irr) = []; +tp(r.irr) = []; + +% Calculate belief x using Bayes' theorem +x = tp.*mu1hat./(tp.*mu1hat + (1-mu1hat).^2); + +% Belief is mu1hat in trials where there is no tone +x(find(tp==0)) = mu1hat(find(tp==0)); + +% Calculate log-probabilities for non-irregular trials +reg = ~ismember(1:n,r.irr); +logp(reg) = -log(1+exp(-be.*(2.*x-1).*(2.*y-1))); +yhat(reg) = x; +res(reg) = (y-x)./sqrt(x.*(1-x)); + +return; diff --git a/HGF/tapas_condhalluc_obs2.m b/HGF/tapas_condhalluc_obs2.m new file mode 100644 index 00000000..af3c663f --- /dev/null +++ b/HGF/tapas_condhalluc_obs2.m @@ -0,0 +1,48 @@ +function [logp, yhat, res] = tapas_condhalluc_obs2(r, infStates, ptrans) +% Calculates the log-probability of response y=1 under the unit-square sigmoid model +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Transform alpha and beta to their native spaces +be = exp(ptrans(1)); +nu = exp(ptrans(2)); + +% Initialize returned log-probabilities as NaNs so that NaN is +% returned for all irregualar trials +n = size(infStates,1); +logp = NaN(n,1); +yhat = NaN(n,1); +res = NaN(n,1); + +% Check input format +if size(r.u,2) ~= 2 + error('tapas:hgf:CondHalluc:InputsIncompatible', 'Inputs incompatible with condhalluc_obs observation model. See tapas_condhalluc_obs_config.m.') +end + +% Get true-positive rate corresponding to stimuli +tp = r.u(:,2); + +% Weed irregular trials out +mu1hat = infStates(:,1,1); +mu1hat(r.irr) = []; +y = r.y(:,1); +y(r.irr) = []; +tp(r.irr) = []; + +% Update belief using precision-weighted prediction error +% with nu the generalized precision +x = mu1hat + 1/(1 + nu)*(tp - mu1hat); + +% Calculate log-probabilities for non-irregular trials +reg = ~ismember(1:n,r.irr); +logp(reg) = -log(1+exp(-be.*(2.*x-1).*(2.*y-1))); +yhat(reg) = x; +res(reg) = (y-x)./sqrt(x.*(1-x)); + +return; diff --git a/HGF/tapas_condhalluc_obs2_config.m b/HGF/tapas_condhalluc_obs2_config.m new file mode 100644 index 00000000..5fd47ad5 --- /dev/null +++ b/HGF/tapas_condhalluc_obs2_config.m @@ -0,0 +1,55 @@ +function c = tapas_condhalluc_obs2_config +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contains the configuration for the response model used to analyze data from conditioned +% hallucination paradigm by Powers & Corlett +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% The rationale for this model is as follows: +% +% TO BE DESCRIBED... +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Config structure +c = struct; + +% Model name +c.model = 'tapas_condhalluc_obs2'; + +% Sufficient statistics of Gaussian parameter priors + +% Beta +c.logbemu = log(48); +c.logbesa = 1; + +% Nu +c.lognumu = log(1); +c.lognusa = 1; + +% Gather prior settings in vectors +c.priormus = [ + c.logbemu,... + c.lognumu,... + ]; + +c.priorsas = [ + c.logbesa,... + c.lognusa,... + ]; + +% Model filehandle +c.obs_fun = @tapas_condhalluc_obs2; + +% Handle to function that transforms observation parameters to their native space +% from the space they are estimated in +c.transp_obs_fun = @tapas_condhalluc_obs2_transp; + +return; diff --git a/HGF/tapas_condhalluc_obs2_namep.m b/HGF/tapas_condhalluc_obs2_namep.m new file mode 100644 index 00000000..5bb772df --- /dev/null +++ b/HGF/tapas_condhalluc_obs2_namep.m @@ -0,0 +1,15 @@ +function pstruct = tapas_condhalluc_obs2_namep(pvec) +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +pstruct = struct; + +pstruct.be = pvec(1); +pstruct.nu = pvec(2); + +return; diff --git a/HGF/tapas_condhalluc_obs2_sim.m b/HGF/tapas_condhalluc_obs2_sim.m new file mode 100644 index 00000000..42307fa4 --- /dev/null +++ b/HGF/tapas_condhalluc_obs2_sim.m @@ -0,0 +1,35 @@ +function y = tapas_condhalluc_obs2_sim(r, infStates, p) +% Simulates responses according to the condhalluc_obs model +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Get parameters +be = p(1); +nu = p(2); + +% Prediction trajectory +mu1hat = infStates(:,1,1); + +% Get true-positive rate corresponding to stimuli +tp = r.u(:,2); + +% Update belief using precision-weighted prediction error +% with nu the generalized precision +x = mu1hat + 1/(1 + nu)*(tp - mu1hat); + +% Apply the logistic sigmoid to the inferred beliefs +prob = tapas_sgm(be.*(2.*x-1),1); + +% Initialize random number generator +rng('shuffle'); + +% Simulate +y = binornd(1, prob); + +return; diff --git a/HGF/tapas_condhalluc_obs2_transp.m b/HGF/tapas_condhalluc_obs2_transp.m new file mode 100644 index 00000000..c6c393a2 --- /dev/null +++ b/HGF/tapas_condhalluc_obs2_transp.m @@ -0,0 +1,19 @@ +function [pvec, pstruct] = tapas_condhalluc_obs2_transp(r, ptrans) +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +pvec = NaN(1,length(ptrans)); +pstruct = struct; + +pvec(1) = exp(ptrans(1)); % be +pstruct.be = pvec(1); + +pvec(2) = exp(ptrans(2)); % nu +pstruct.nu = pvec(2); + +return; diff --git a/HGF/tapas_condhalluc_obs_config.m b/HGF/tapas_condhalluc_obs_config.m new file mode 100644 index 00000000..0ffd8c2a --- /dev/null +++ b/HGF/tapas_condhalluc_obs_config.m @@ -0,0 +1,94 @@ +function c = tapas_condhalluc_obs_config +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% Contains the configuration for the response model used to analyze data from conditioned +% hallucination paradigm by Powers & Corlett +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% The rationale for this model is as follows: +% +% We apply decision noise (i.e., a logistic sigmoid, see below) to the probability that the subject +% says “yes” on a given trial: +% +% p( yes | belief ) = sigmoid( belief ), +% +% where belief = p( tone | percept, light ), +% +% and where in turn “percept” is the subjective experience of hearing (or not hearing) a tone, while +% “tone” is the objective presentation of a tone, and “light” is the presentation of a light. +% +% In trial where there is a tone, we may use Bayes’ theorem to get the belief: +% +% belief = p( tone | percept, light ) = p( percept | tone )*p( tone | light ) / (p( percept | +% tone )*p( tone | light ) + p( percept | no tone )*p( no tone | light )) +% +% Unpacking the various ingredients, we have +% +% - p( percept | tone ) is given by experimental design: the true positive rate of the tone +% presented without light at each trial - 1/4, 1/2, or 3/4 +% +% - p( tone | light ) is the prior from learning using the HGF: mu1hat +% +% - p( percept | no tone ) is the false positive rate, which we can take to be 1 - mu1hat +% +% - p( no tone | light ) is the other half of the prior: 1 - mu1hat +% +% In trials where there is no tone, the belief is mu1hat. +% +% The logistic sigmoid is +% +% f(x) = 1/(1+exp(-beta*(v1-v0))), +% +% where v1 and v0 are the values of options 1 and 0, respectively, and beta > 0 is a parameter that +% determines the slope of the sigmoid. Beta is sometimes referred to as the (inverse) decision +% temperature. In the formulation above, it represents the probability of choosing option 1. +% Reversing the roles of v1 and v0 yields the probability of choosing option 0. +% +% Beta can be interpreted as inverse decision noise. To have a shrinkage prior on this, choose a +% high value. It is estimated log-space since it has a natural lower bound at zero. In general, v1 +% and v0 can be any real numbers. +% +% This observation model expects the first column of the input matrix to contain the outcomes +% (corresponding in this case to the choices of the subject): 1 or 0 for each trial. The SECOND +% COLUMN of the input matrix is expected to contain the true-positive rate of the stimulus for +% each trial. The response matrix only contains one column consisting of the choices of the +% subject. This means that it will be identical to the first column of the input matrix. +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Config structure +c = struct; + +% Model name +c.model = 'tapas_condhalluc_obs'; + +% Sufficient statistics of Gaussian parameter priors + +% Beta +c.logbemu = log(48); +c.logbesa = 1; + +% Gather prior settings in vectors +c.priormus = [ + c.logbemu,... + ]; + +c.priorsas = [ + c.logbesa,... + ]; + +% Model filehandle +c.obs_fun = @tapas_condhalluc_obs; + +% Handle to function that transforms observation parameters to their native space +% from the space they are estimated in +c.transp_obs_fun = @tapas_condhalluc_obs_transp; + +return; diff --git a/HGF/tapas_condhalluc_obs_namep.m b/HGF/tapas_condhalluc_obs_namep.m new file mode 100644 index 00000000..d62297c8 --- /dev/null +++ b/HGF/tapas_condhalluc_obs_namep.m @@ -0,0 +1,14 @@ +function pstruct = tapas_condhalluc_obs_namep(pvec) +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +pstruct = struct; + +pstruct.be = pvec(1); + +return; diff --git a/HGF/tapas_condhalluc_obs_sim.m b/HGF/tapas_condhalluc_obs_sim.m new file mode 100644 index 00000000..131332a4 --- /dev/null +++ b/HGF/tapas_condhalluc_obs_sim.m @@ -0,0 +1,36 @@ +function y = tapas_condhalluc_obs_sim(r, infStates, p) +% Simulates responses according to the condhalluc_obs model +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2016 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Inverse decision temparature beta is the only parameter +be = p; + +% Prediction trajectory +mu1hat = infStates(:,1,1); + +% Get true-positive rate corresponding to stimuli +tp = r.u(:,2); + +% Calculate belief x using Bayes' theorem +x = tp.*mu1hat./(tp.*mu1hat + (1-mu1hat).^2); + +% Belief is mu1hat in trials where there is no tone +x(find(tp==0)) = mu1hat(find(tp==0)); + +% Apply the logistic sigmoid to the inferred beliefs +prob = tapas_sgm(be.*(2.*x-1),1); + +% Initialize random number generator +rng('shuffle'); + +% Simulate +y = binornd(1, prob); + +return; diff --git a/HGF/tapas_condhalluc_obs_transp.m b/HGF/tapas_condhalluc_obs_transp.m new file mode 100644 index 00000000..2abe4cb1 --- /dev/null +++ b/HGF/tapas_condhalluc_obs_transp.m @@ -0,0 +1,16 @@ +function [pvec, pstruct] = tapas_condhalluc_obs_transp(r, ptrans) +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +pvec = NaN(1,length(ptrans)); +pstruct = struct; + +pvec(1) = exp(ptrans(1)); % be +pstruct.be = pvec(1); + +return; \ No newline at end of file diff --git a/HGF/tapas_hgf_ar1_binary.m b/HGF/tapas_hgf_ar1_binary.m index f46815c1..13b3bfc9 100644 --- a/HGF/tapas_hgf_ar1_binary.m +++ b/HGF/tapas_hgf_ar1_binary.m @@ -14,7 +14,7 @@ % transformed space, and 'trans' is a flag indicating this. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -109,7 +109,7 @@ % 1st level % ~~~~~~~~~ % Prediction - muhat(k,1) = tapas_sgm(muhat(k,2), 1); + muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1); % Precision of prediction pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1))); @@ -129,8 +129,8 @@ pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2))); % Updates - pi(k,2) = pihat(k,2) +1/pihat(k,1); - mu(k,2) = muhat(k,2) +1/pi(k,2) *da(k,1); + pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1); + mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1); % Volatility prediction error da(k,2) = (1/pi(k,2) +(mu(k,2) -muhat(k,2))^2) *pihat(k,2) -1; @@ -203,7 +203,7 @@ end % Implied learning rate at the first level -sgmmu2 = tapas_sgm(mu(:,2), 1); +sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1); dasgmmu2 = u -sgmmu2; lr1 = diff(sgmmu2)./dasgmmu2(2:n,1); lr1(da(2:n,1)==0) = 0; diff --git a/HGF/tapas_hgf_ar1_binary_config.m b/HGF/tapas_hgf_ar1_binary_config.m index d1c9db83..686ac97a 100644 --- a/HGF/tapas_hgf_ar1_binary_config.m +++ b/HGF/tapas_hgf_ar1_binary_config.m @@ -9,6 +9,15 @@ % Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation % for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39. % +% The binary HGF model has since been augmented with a positive factor kappa1 which +% scales the second level with respect to the first, i.e., the relation between the +% first and second level is +% +% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid. +% +% By default, kappa1 is fixed to 1, leading (apart from the AR(1) process) to the +% model introduced in Mathys et al. (2011). +% % This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011)); % for continuous inputs, refer to tapas_hgf_config. % @@ -93,7 +102,7 @@ % the LME increased, so you had a better model. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -148,11 +157,12 @@ % Kappas % Format: row vector of length n_levels-1. -% Undefined (therefore NaN) at the first level. -% This should be fixed (preferably to 1) if the observation model -% does not use mu_i+1 (kappa then determines the scaling of x_i+1). -c.logkamu = [NaN, log(1)]; -c.logkasa = [NaN, 0]; +% Fixing log(kappa1) to log(1) leads to the original HGF model. +% Higher log(kappas) should be fixed (preferably to log(1)) if the +% observation model does not use mu_i+1 (kappa then determines the +% scaling of x_i+1). +c.logkamu = [log(1), log(1)]; +c.logkasa = [ 0, 0]; % Omegas % Format: row vector of length n_levels. diff --git a/HGF/tapas_hgf_binary.m b/HGF/tapas_hgf_binary.m index cb572b5b..c4504327 100644 --- a/HGF/tapas_hgf_binary.m +++ b/HGF/tapas_hgf_binary.m @@ -13,7 +13,7 @@ % transformed space, and 'trans' is a flag indicating this. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -107,7 +107,7 @@ % 1st level % ~~~~~~~~~ % Prediction - muhat(k,1) = tapas_sgm(muhat(k,2), 1); + muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1); % Precision of prediction pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1))); @@ -127,8 +127,8 @@ pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2))); % Updates - pi(k,2) = pihat(k,2) +1/pihat(k,1); - mu(k,2) = muhat(k,2) +1/pi(k,2) *da(k,1); + pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1); + mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1); % Volatility prediction error da(k,2) = (1/pi(k,2) +(mu(k,2) -muhat(k,2))^2) *pihat(k,2) -1; @@ -201,7 +201,7 @@ end % Implied learning rate at the first level -sgmmu2 = tapas_sgm(mu(:,2), 1); +sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1); dasgmmu2 = u -sgmmu2; lr1 = diff(sgmmu2)./dasgmmu2(2:n,1); lr1(da(2:n,1)==0) = 0; diff --git a/HGF/tapas_hgf_binary_condhalluc_plotTraj.m b/HGF/tapas_hgf_binary_condhalluc_plotTraj.m new file mode 100644 index 00000000..0afd8654 --- /dev/null +++ b/HGF/tapas_hgf_binary_condhalluc_plotTraj.m @@ -0,0 +1,106 @@ +function tapas_hgf_binary_condhalluc_plotTraj(r) +% Plots the estimated or generated trajectories for the binary HGF perceptual model +% Usage example: est = tapas_fitModel(responses, inputs); tapas_hgf_binary_plotTraj(est); +% +% -------------------------------------------------------------------------------------------------- +% Copyright (C) 2015 Christoph Mathys, TNU, UZH & ETHZ +% +% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + +% Optional plotting of standard deviations (true or false) +plotsd = true; + +% Optional plotting of responses (true or false) +ploty = true; + +% Set up display +scrsz = get(0,'screenSize'); +outerpos = [0.2*scrsz(3),0.2*scrsz(4),0.8*scrsz(3),0.8*scrsz(4)]; +figure(... + 'OuterPosition', outerpos,... + 'Name', 'HGF trajectories'); + +% Time axis +t = ones(1,size(r.u,1)); + +ts = cumsum(t); +ts = [0, ts]; + +% Number of levels +try + l = r.c_prc.n_levels; +catch + l = (length(r.p_prc.p)+1)/5; +end + +% Upper levels +for j = 1:l-1 + + % Subplots + subplot(l,1,j); + + if plotsd == true + upperprior = r.p_prc.mu_0(l-j+1) +sqrt(r.p_prc.sa_0(l-j+1)); + lowerprior = r.p_prc.mu_0(l-j+1) -sqrt(r.p_prc.sa_0(l-j+1)); + upper = [upperprior; r.traj.mu(:,l-j+1)+sqrt(r.traj.sa(:,l-j+1))]; + lower = [lowerprior; r.traj.mu(:,l-j+1)-sqrt(r.traj.sa(:,l-j+1))]; + + plot(0, upperprior, 'ob', 'LineWidth', 1); + hold all; + plot(0, lowerprior, 'ob', 'LineWidth', 1); + fill([ts, fliplr(ts)], [(upper)', fliplr((lower)')], ... + 'b', 'EdgeAlpha', 0, 'FaceAlpha', 0.15); + end + plot(ts, [r.p_prc.mu_0(l-j+1); r.traj.mu(:,l-j+1)], 'b', 'LineWidth', 2); + hold all; + plot(0, r.p_prc.mu_0(l-j+1), 'ob', 'LineWidth', 2); % prior + xlim([0 ts(end)]); + title(['Posterior expectation of x_' num2str(l-j+1)], 'FontWeight', 'bold'); + ylabel(['\mu_', num2str(l-j+1)]); +end + +% Input level +subplot(l,1,l); + +plot(ts, [tapas_sgm(r.p_prc.mu_0(2), 1); tapas_sgm(r.traj.mu(:,2), 1)], 'r', 'LineWidth', 2); +hold all; +% Get true-positive rate corresponding to stimuli +tp = r.u(:,2); +mu1hat = r.traj.muhat(:,1); +% Calculate belief x using Bayes' theorem +tonebelief = tp.*mu1hat./(tp.*mu1hat + (1-mu1hat).^2); +% Belief is mu1hat in trials where there is no tone +tonebelief(find(tp==0)) = mu1hat(find(tp==0)); +% Plot tone belief +plot(ts, [NaN; tonebelief], 'b', 'LineWidth', 2); +plot(0, tapas_sgm(r.p_prc.mu_0(2), 1), 'or', 'LineWidth', 2); % prior +plot(ts(2:end), r.u(:,2), '.', 'Color', [0 0.6 0]); % inputs +plot(ts(2:end), r.traj.wt(:,1), 'k') % implied learning rate +if (ploty == true) && ~isempty(find(strcmp(fieldnames(r),'y'))) && ~isempty(r.y) + y = r.y(:,1) -0.5; y = 1.16 *y; y = y +0.5; % stretch + if ~isempty(find(strcmp(fieldnames(r),'irr'))) + y(r.irr) = NaN; % weed out irregular responses + plot(ts(r.irr), 1.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses + plot(ts(r.irr), -0.08.*ones([1 length(r.irr)]), 'x', 'Color', [1 0.7 0], 'Markersize', 11, 'LineWidth', 2); % irregular responses + end + plot(ts(2:end), y, '.', 'Color', [1 0.7 0]); % responses + title(['Response y (orange), input u (green), learning rate (fine black), posterior belief (blue), and posterior expectation of input s(\mu_2) ', ... + '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ... + num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ... + 'FontWeight', 'bold'); + ylabel('y, u, s(\mu_2)'); + axis([0 ts(end) -0.15 1.15]); +else + title(['Input u (green), learning rate (fine black), posterior belief (blue), and posterior expectation of input s(\mu_2) ', ... + '(red) for \rho=', num2str(r.p_prc.rho(2:end)), ', \kappa=', ... + num2str(r.p_prc.ka(2:end)), ', \omega=', num2str(r.p_prc.om(2:end))], ... + 'FontWeight', 'bold'); + ylabel('u, s(\mu_2)'); + axis([0 ts(end) -0.1 1.1]); +end +plot(ts(2:end), 0.5, 'k'); +xlabel('Trial number'); +hold off; diff --git a/HGF/tapas_hgf_binary_config.m b/HGF/tapas_hgf_binary_config.m index 82812bb4..7f5d3860 100644 --- a/HGF/tapas_hgf_binary_config.m +++ b/HGF/tapas_hgf_binary_config.m @@ -9,6 +9,15 @@ % Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation % for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39. % +% The binary HGF model has since been augmented with a positive factor kappa1 which +% scales the second level with respect to the first, i.e., the relation between the +% first and second level is +% +% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid. +% +% By default, kappa1 is fixed to 1, leading exactly to the model introduced in +% Mathys et al. (2011). +% % This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011)); % for continuous inputs, refer to tapas_hgf_config. % @@ -86,7 +95,7 @@ % the LME increased, so you had a better model. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2012-2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -133,11 +142,12 @@ % Kappas % Format: row vector of length n_levels-1. -% Undefined (therefore NaN) at the first level. -% This should be fixed (preferably to 1) if the observation model -% does not use mu_i+1 (kappa then determines the scaling of x_i+1). -c.logkamu = [NaN, log(1)]; -c.logkasa = [NaN, 0]; +% Fixing log(kappa1) to log(1) leads to the original HGF model. +% Higher log(kappas) should be fixed (preferably to log(1)) if the +% observation model does not use mu_i+1 (kappa then determines the +% scaling of x_i+1). +c.logkamu = [log(1), log(1)]; +c.logkasa = [ 0, 0]; % Omegas % Format: row vector of length n_levels. diff --git a/HGF/tapas_hgf_binary_mab.m b/HGF/tapas_hgf_binary_mab.m index c4b6a245..73c40883 100644 --- a/HGF/tapas_hgf_binary_mab.m +++ b/HGF/tapas_hgf_binary_mab.m @@ -14,7 +14,7 @@ % transformed space, and 'trans' is a flag indicating this. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -120,7 +120,7 @@ % 1st level % ~~~~~~~~~ % Prediction - muhat(k,1,:) = tapas_sgm(muhat(k,2,:), 1); + muhat(k,1,:) = tapas_sgm(ka(1) *muhat(k,2,:), 1); % Precision of prediction pihat(k,1,:) = 1/(muhat(k,1,:).*(1 -muhat(k,1,:))); @@ -144,10 +144,10 @@ % Updates pi(k,2,:) = pihat(k,2,:); - pi(k,2,y(k)) = pihat(k,2,y(k)) +1/pihat(k,1,y(k)); + pi(k,2,y(k)) = pihat(k,2,y(k)) +ka(1)^2/pihat(k,1,y(k)); mu(k,2,:) = muhat(k,2,:); - mu(k,2,y(k)) = muhat(k,2,y(k)) +1/pi(k,2,y(k)) *da(k,1); + mu(k,2,y(k)) = muhat(k,2,y(k)) +ka(1)/pi(k,2,y(k)) *da(k,1); % Volatility prediction error da(k,2) = (1/pi(k,2,y(k)) +(mu(k,2,y(k)) -muhat(k,2,y(k)))^2) *pihat(k,2,y(k)) -1; @@ -262,7 +262,7 @@ mu2obs = mu2(sub2ind(size(mu2), (1:size(mu2,1))', y)); mu1hat = squeeze(muhat(:,1,:)); mu1hatobs = mu1hat(sub2ind(size(mu1hat), (1:size(mu1hat,1))', y)); -upd1 = tapas_sgm(mu2obs,1) -mu1hatobs; +upd1 = tapas_sgm(ka(1)*mu2obs,1) -mu1hatobs; lr1 = upd1./da(:,1); % Create result data structure diff --git a/HGF/tapas_hgf_binary_mab_config.m b/HGF/tapas_hgf_binary_mab_config.m index d9262a8b..2ee81b5f 100644 --- a/HGF/tapas_hgf_binary_mab_config.m +++ b/HGF/tapas_hgf_binary_mab_config.m @@ -9,6 +9,15 @@ % Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation % for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39. % +% The binary HGF model has since been augmented with a positive factor kappa1 which +% scales the second level with respect to the first, i.e., the relation between the +% first and second level is +% +% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid. +% +% By default, kappa1 is fixed to 1, leading exactly to the model introduced in +% Mathys et al. (2011). +% % This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011)); % for continuous inputs, refer to tapas_hgf_config. % @@ -86,7 +95,7 @@ % the LME increased, so you had a better model. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2013 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2013-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -143,11 +152,12 @@ % Kappas % Format: row vector of length n_levels-1. -% Undefined (therefore NaN) at the first level. -% This should be fixed (preferably to 1) if the observation model -% does not use mu_i+1 (kappa then determines the scaling of x_i+1). -c.logkamu = [NaN, log(1)]; -c.logkasa = [NaN, 0]; +% Fixing log(kappa1) to log(1) leads to the original HGF model. +% Higher log(kappas) should be fixed (preferably to log(1)) if the +% observation model does not use mu_i+1 (kappa then determines the +% scaling of x_i+1). +c.logkamu = [log(1), log(1)]; +c.logkasa = [ 0, 0]; % Omegas % Format: row vector of length n_levels. diff --git a/HGF/tapas_hgf_binary_pu.m b/HGF/tapas_hgf_binary_pu.m index 66729d27..17fea1ad 100644 --- a/HGF/tapas_hgf_binary_pu.m +++ b/HGF/tapas_hgf_binary_pu.m @@ -110,7 +110,7 @@ % 1st level % ~~~~~~~~~ % Prediction - muhat(k,1) = tapas_sgm(muhat(k,2), 1); + muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1); % Precision of prediction pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1))); @@ -132,11 +132,11 @@ pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2))); % Updates - pi(k,2) = pihat(k,2) +1/pihat(k,1); - mu(k,2) = muhat(k,2) +1/pi(k,2) *da(k,1); + pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1); + mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1); % Implied posterior precision at first level - sgmmu2 = tapas_sgm(mu(k,2), 1); + sgmmu2 = tapas_sgm(ka(1) *mu(k,2), 1); pi(k,1) = pi(k,2)/(sgmmu2*(1-sgmmu2)); % Volatility prediction error @@ -210,7 +210,7 @@ end % Implied learning rate at the first level -sgmmu2 = tapas_sgm(mu(:,2), 1); +sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1); lr1 = diff(sgmmu2)./da(2:n,1); lr1(da(2:n,1)==0) = 0; diff --git a/HGF/tapas_hgf_binary_pu_config.m b/HGF/tapas_hgf_binary_pu_config.m index 92342844..a5f2e2e0 100644 --- a/HGF/tapas_hgf_binary_pu_config.m +++ b/HGF/tapas_hgf_binary_pu_config.m @@ -9,6 +9,15 @@ % Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation % for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39. % +% The binary HGF model has since been augmented with a positive factor kappa1 which +% scales the second level with respect to the first, i.e., the relation between the +% first and second level is +% +% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid. +% +% By default, kappa1 is fixed to 1, leading exactly to the model introduced in +% Mathys et al. (2011). +% % This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011)); % for continuous inputs, refer to tapas_hgf_config. % @@ -94,7 +103,7 @@ % the LME increased, so you had a better model. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2012-2015 Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2012-2017 Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -141,11 +150,12 @@ % Kappas % Format: row vector of length n_levels-1. -% Undefined (therefore NaN) at the first level. -% This should be fixed (preferably to 1) if the observation model -% does not use mu_i+1 (kappa then determines the scaling of x_i+1). -c.logkamu = [NaN, log(1)]; -c.logkasa = [NaN, 0]; +% Fixing log(kappa1) to log(1) leads to the original HGF model. +% Higher log(kappas) should be fixed (preferably to log(1)) if the +% observation model does not use mu_i+1 (kappa then determines the +% scaling of x_i+1). +c.logkamu = [log(1), log(1)]; +c.logkasa = [ 0, 0]; % Omegas % Format: row vector of length n_levels. diff --git a/HGF/tapas_hgf_binary_pu_tbt.m b/HGF/tapas_hgf_binary_pu_tbt.m index 0c2d2c2e..d487d170 100644 --- a/HGF/tapas_hgf_binary_pu_tbt.m +++ b/HGF/tapas_hgf_binary_pu_tbt.m @@ -13,7 +13,7 @@ % transformed space, and 'trans' is a flag indicating this. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2016 Rebecca Lawson, Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2016-2017 Rebecca Lawson, Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -110,7 +110,7 @@ % 1st level % ~~~~~~~~~ % Prediction - muhat(k,1) = tapas_sgm(muhat(k,2), 1); + muhat(k,1) = tapas_sgm(ka(1) *muhat(k,2), 1); % Precision of prediction pihat(k,1) = 1/(muhat(k,1)*(1 -muhat(k,1))); @@ -132,11 +132,11 @@ pihat(k,2) = 1/(1/pi(k-1,2) +exp(ka(2) *mu(k-1,3) +om(2))); % Updates - pi(k,2) = pihat(k,2) +1/pihat(k,1); - mu(k,2) = muhat(k,2) +1/pi(k,2) *da(k,1); + pi(k,2) = pihat(k,2) +ka(1)^2/pihat(k,1); + mu(k,2) = muhat(k,2) +ka(1)/pi(k,2) *da(k,1); % Implied posterior precision at first level - sgmmu2 = tapas_sgm(mu(k,2), 1); + sgmmu2 = tapas_sgm(ka(1) *mu(k,2), 1); pi(k,1) = pi(k,2)/(sgmmu2*(1-sgmmu2)); % Volatility prediction error @@ -210,7 +210,7 @@ end % Implied learning rate at the first level -sgmmu2 = tapas_sgm(mu(:,2), 1); +sgmmu2 = tapas_sgm(ka(1) *mu(:,2), 1); lr1 = diff(sgmmu2)./da(2:n,1); lr1(da(2:n,1)==0) = 0; diff --git a/HGF/tapas_hgf_binary_pu_tbt_config.m b/HGF/tapas_hgf_binary_pu_tbt_config.m index 17e33cc9..1a818be6 100644 --- a/HGF/tapas_hgf_binary_pu_tbt_config.m +++ b/HGF/tapas_hgf_binary_pu_tbt_config.m @@ -11,6 +11,14 @@ % Mathys C, Daunizeau J, Friston, KJ, and Stephan KE. (2011). A Bayesian foundation % for individual learning under uncertainty. Frontiers in Human Neuroscience, 5:39. % +% The binary HGF model has since been augmented with a positive factor kappa1 which +% scales the second level with respect to the first, i.e., the relation between the +% first and second level is +% +% p(x1=1|x2) = s(kappa1*x2), where s(.) is the logistic sigmoid. +% +% By default, kappa1 is fixed to 1, leading to the model introduced in Mathys et al. (2011). +% % This file refers to BINARY inputs (Eqs 1-3 in Mathys et al., (2011)); % for continuous inputs, refer to tapas_hgf_config. % @@ -96,7 +104,7 @@ % the LME increased, so you had a better model. % % -------------------------------------------------------------------------------------------------- -% Copyright (C) 2016 Rebecca Lawson, Christoph Mathys, TNU, UZH & ETHZ +% Copyright (C) 2016-2017 Rebecca Lawson, Christoph Mathys, TNU, UZH & ETHZ % % This file is part of the HGF toolbox, which is released under the terms of the GNU General Public % Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL @@ -143,11 +151,12 @@ % Kappas % Format: row vector of length n_levels-1. -% Undefined (therefore NaN) at the first level. -% This should be fixed (preferably to 1) if the observation model -% does not use mu_i+1 (kappa then determines the scaling of x_i+1). -c.logkamu = [NaN, log(1)]; -c.logkasa = [NaN, 0]; +% Fixing log(kappa1) to log(1) leads to the original HGF model. +% Higher log(kappas) should be fixed (preferably to log(1)) if the +% observation model does not use mu_i+1 (kappa then determines the +% scaling of x_i+1). +c.logkamu = [log(1), log(1)]; +c.logkasa = [ 0, 0]; % Omegas % Format: row vector of length n_levels. diff --git a/PhysIO/.gitmodules b/PhysIO/.gitmodules new file mode 100644 index 00000000..e87034b3 --- /dev/null +++ b/PhysIO/.gitmodules @@ -0,0 +1,4 @@ +[submodule "wikidocs"] + path = wikidocs + url = git@tnurepository.ethz.ch:physio/physio-public.wiki.git + branch = master diff --git a/PhysIO/CHANGELOG.md b/PhysIO/CHANGELOG.md new file mode 100644 index 00000000..b63b1e28 --- /dev/null +++ b/PhysIO/CHANGELOG.md @@ -0,0 +1,108 @@ +RELEASE INFORMATION +=============================== + +Current Release +--------------- + +PhysIO_Toolbox_R2017.2 + +August 03, 2017 + +Minor Release Notes (R2017.2) +----------------------------- + +- Included Markdown-based documentation via Wiki (also CITATION, LICENSE, CHANGELOG.md) +- Included FAQ in Wiki +- Split git repositories into public, dev, examples, and added wiki, to disentangle development from deployed toolbox code and data +- Bugfix and Typo correction +- Philips SCANPYHSLOG for their software release 5.1.7. + +Minor Release Notes (R2017.1) +----------------------------- + +- Substantially improved Siemens interface, both for VB/VD and 3T/7T releases + - several bugfixes + - based on extensive user feedback from Berlin and Brisbane +- New functionality tapas_physio_overlay_contrasts.m to display non-physio + contrasts automatically as well + +Major Release Notes (r904 / R2016.1) +------------------------------------ + +- Software version for accepted PhysIO Toolbox Paper: doi:10.1016/j.jneumeth.2016.10.019 +- Tested and expanded versions of examples +- Improved stability by bugfixes and compatibility to Matlab R2016 +- Slice-wise regressor creation +- Detection of constant physiological time series (detachment, clipping) +- Refactoring of report_contrasts and compute_tsnr_gains as standalone functionality +- Improved Read-in capabilities (Siemens respiration data, BioPac .mat) +- Migration from svn (r904) to git (tnurepository) for version control + +Major Release Notes (r835) +-------------------------- + +- Software version for Toolbox Paper submission +- Noise ROIs modeling +- Extended motion models (24 parameters, Volterra expansion) +- HRV/RVT models with optional multiple delay regressors +- Report_contrasts with automatic contrast generation for all regressor groups +- compute_tsnr_gains for individual physiological regressor groups +- consistent module naming (scan_timing, preproc) +- Visualisation improvement (color schemes, legends) + +Minor Release Notes (r666) +-------------------------- + +- Compatibility tested for SPM12, small bugfixes Batch Dependencies +- Cleaner Batch Interface with grouped sub-menus (cfg_choice) +- new model: 'none' to just read out physiological raw data and preprocess, + without noise modelling +- Philips: Scan-timing via gradient log now automatized (gradient_log_auto) +- Siemens: Tics-Logfile read-in (proprietary, needs Siemens-agreement) +- All peak detections (cardiac/respiratory) now via auto_matched algorithm +- Adapt plots/saving for Matlab R2014b + +Major Release Notes (r534) +-------------------------- + +- Read-in of Siemens plain text log files; new example dataset for Siemens +- Speed up and debugging of auto-detection method for noisy cardiac data => new method thresh.cardiac.initial_cpulse_select.method = ???auto_matched??? +- Error handling for temporary breathing belt failures (Eduardo Aponte, TNU Zurich) +- slice-wise regressors can be created by setting sqpar.onset_slice to a index vector of slices + +Major Release Notes (r497) +-------------------------- + +- SPM matlabbatch GUI implemented (Call via Batch -> SPM -> Tools -> TAPAS PhysIO Toolbox) +- improved, automatic heartbeat detection for noisy ECG now standard for ECG and Pulse oximetry (courtesy of Steffen Bollmann) +- QuickStart-Manual and PhysIO-Background presentation expanded/updated +- job .m/.mat-files created for all example datasets +- bugfixes cpulse-initial-select method-handling (auto/manual/load) + +Major Release Notes (r429) +-------------------------- + +- Cardiac and Respiratory response function regressors integrated in workflow (heart rate and breathing volume computation) +- Handling of Cardiac and Respiratory Logfiles only +- expanded documentation (Quickstart.pdf and Handbook.pdf) +- read-in of custom log files, e.g. for BrainVoyager peripheral data +- more informative plots and commenting (especially in tapas_physio_new). + +Minor Release Notes (r354) +-------------------------- + +- computation of heart and breathing rate in Philips/PPU/main_PPU.m +- prefix of functions with tapas_* + +Major Release Notes (r241) +-------------------------- + +- complete modularization of reading/preprocessing/regressor creation for peripheral physiological data +- manual selection of missed heartbeats in ECG/pulse oximetry (courtesy of Jakob Heinzle) +- support for logfiles from GE scanners (courtesy of Steffen Bollmann, KiSpi Zuerich) +- improved detection of pulse oximetry peaks (courtesy of Steffen Bollmann) +- improved documentation +- consistent function names (prefixed by "physio_") + +NOTE: Your main_ECG/PPU.m etc. scripts from previous versions (<=r159) will not work with this one any more. Please adapt one of the example scripts for your needs (~5 min of work). The main benefit of this version is a complete new variable structure that is more sustainable and makes the code more readable. + diff --git a/PhysIO/CITATION b/PhysIO/CITATION new file mode 100644 index 00000000..6add7bf3 --- /dev/null +++ b/PhysIO/CITATION @@ -0,0 +1,10 @@ +Please cite this software as follows: + +Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., +Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., +Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise +in fMRI Data. +Journal of Neuroscience Methods 276, 56–72. +doi:10.1016/j.jneumeth.2016.10.019 + +For a detailed snippet to be included in the methods section of your manuscript, see the FAQ in the documentation. \ No newline at end of file diff --git a/PhysIO/code/COPYING b/PhysIO/LICENSE similarity index 100% rename from PhysIO/code/COPYING rename to PhysIO/LICENSE diff --git a/PhysIO/README.md b/PhysIO/README.md new file mode 100644 index 00000000..a8f27a08 --- /dev/null +++ b/PhysIO/README.md @@ -0,0 +1,228 @@ +TAPAS PhysIO Toolbox Version 2017 +================================= + +> Copyright (C) 2012-2017 Lars Kasper + +> Translational Neuromodeling Unit (TNU) + +> Institute for Biomedical Engineering + +> University of Zurich and ETH Zurich + +Copying +------- + +The PhysIO Toolbox is free software: you can redistribute it and/or +modify it under the terms of the GNU General Public License as +published by the Free Software Foundation, either version 3 of the +License, or (at your option) any later version. + +This program is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program (see the file COPYING). If not, see +. + +Purpose +------- + +The PhysIO Toolbox provides model-based physiological noise correction of +fMRI data using peripheral measures of respiration and cardiac pulsation. +It incorporates noise models of cardiac/respiratory phase (RETROICOR, +Glover et al. 2000), as well as heart rate variability and respiratory +volume per time (cardiac response function, Chang et. al, 2009, respiratory +response function, Birn et al. 2006). The toolbox is usable via the SPM +batch editor, performs automatic pre-processing of noisy peripheral data +and outputs nuisance regressor files directly suitable for SPM +(multiple_regressors.txt). + +Installation +------------ + +### Matlab ### +- Unzip the TAPAS archive +- Add tapas/physio/code to your matlab path + +### SPM ### +- Certain functionality (Batch Editor GUI, pipeline dependencies, model assessment via F-contrasts) require the installation of SPM +- Afterwards, the PhysIO Toolbox has to be registered as an SPM toolbox by copying the `physio/code` folder to `spm/toolbox/physio` + + +Getting Started +--------------- + +Run `example_main_ECG3T.m` in subdirectory `Philips/ECG3T` of the toolbox example repository `physio-examples`. +See subdirectory `physio/docs` and next section of this document. + + +Getting Help/Documentation +-------------------------- + +Several documentation files are provided with this toolbox. They have the extension .md (markdown), i.e. are plain text files, but can be conveniently viewed online as the github/gitlab Wiki. +You can find them in `physio/wikidocs`. + +Alternatively, a pdf and html converted version of the following files is found in `physio/docs/documentation.pdf (or .html)` + +List of Documentation files: +- README.md: this file, purpose, installation, getting started, pointer to more help +- FAQ.md: Frequently asked questions (for users) +- QUICKSTART.md: Example script and how to use on test data, Intro to Batch Editor GUI +- MANUAL.md: Reference Manual (mostly for developers) listing all functions, and rationales of the toolbox, disecting its modular structure + - not provided yet; see the old version in `physio/docs/QuickStart_PhysIO_Toolbox.pdf` for source code documentation +- HOME.md: Landing Page of Wiki. Navigation to all other files and this explanation +- EXAMPLES.md: List and explanation of all examples +- CHANGELOG.md: List of all toolbox versions and the respective release notes, i.e. major changes in functionality, bugfixes etc. + + +Background +---------- + +The PhysIO Toolbox provides physiological noise correction for fMRI-data +from peripheral measures (ECG/pulse oximetry, breathing belt). It is +model-based, i.e. creates nuisance regressors from the physiological +monitoring that can enter a General Linear Model (GLM) analysis, e.g. +SPM8/12. Furthermore, for scanner vendor logfiles (PHILIPS, GE, Siemens), +it provides means to statistically assess peripheral data (e.g. heart rate variability) +and recover imperfect measures (e.g. distorted R-peaks of the ECG). + +Facts about physiological noise in fMRI: +- Physiological noise can explain 20-60 % of variance in fMRI voxel time + series (Birn2006, Hutton2011, Harvey2008) + - Physiological noise affects a lot of brain regions (s. figure, e.g. + brainstem or OFC), especially next to CSF, arteries (Hutton2011). + - If not accounted for, this is a key factor limiting sensitivity for effects of interest. +- Physiological noise contributions increase with field strength; they + become a particular concern at and above 3 Tesla (Kasper2009, Hutton2011). +- In resting state fMRI, disregarding physiological noise leads to wrong + connectivity results (Birn2006). + +Therefore, some kind of physiological noise correction is highly recommended for every statistical fMRI analysis. + +Model-based correction of physiological noise: +- Physiological noise can be decomposed into periodic time series following + heart rate and breathing cycle. +- The Fourier expansion of cardiac and respiratory phases was introduced as + RETROICOR (RETROspective Image CORrection, Glover2000, + see also Josephs1997). +- These Fourier Terms can enter a General Linear Model (GLM) as nuisance + regressors, analogous to movement parameters. +- As the physiological noise regressors augment the GLM and explain + variance in the time series, they increase sensitivity in all contrasts + of interest. + + +Features of this Toolbox +------------------------ + +### Physiological Noise Modeling ### + +- Modeling physiological noise regressors from peripheral data + (breathing belt, ECG, pulse oximeter) + - State of the art RETROICOR cardiac and respiratory phase expansion + - Cardiac response function (Chang et al, 2009) and respiratory response + function (Birn et al. 2006) modelling of heart-rate variability and + respiratory volume per time influence on physiological noise + - Flexible expansion orders to model different contributions of cardiac, + respiratory and interaction terms (see Harvey2008, Hutton2011) +- Data-driven noise regressors + - PCA extraction from nuisance ROIs (CSF, white matter), similar to aCompCor (Behzadi2007) + +### Automatization and Performance Assessment ### +- Automatic creation of nuisance regressors, full integration into standard + GLMs, tested for SPM8/12 ("multiple_regressors.mat") +- Integration in SPM Batch Editor: GUI for parameter input, dependencies to integrate physiological noise correction in preprocessing pipeline +- Performance Assessment: Automatic F-contrast and tSNR Map creation and display for groups of physiological noise regressors, using SPM GLM tools + +### Flexible Read-in ### + +The toolbox is dedicated to seamless integration into a clinical research s +etting and therefore offers correction methods to recover physiological +data from imperfect peripheral measures. + +- General Electric +- Philips SCANPHYSLOG files (all versions from release 2.6 to 5.3) +- Siemens VB (files `.ecg`, `.resp`, `.puls` +- Siemens VD (files (`*_ECG.log`, `*_RESP.log`, `*_PULS.log`) +- Biopac .mat-export + - assuming the following variables (as columns): `data`, `isi`, `isi_units`, `labels`, `start_sample`, `units` + - See `tapas_physio_read_physlogfiles_biopac_mat.m` for details +- Custom logfiles: should contain one amplitude value per line, one logfile per device. Sampling interval(s) are provided as a separate parameter to the toolbox. + + + +Compatibility and Support +------------------------- + +- Matlab Toolbox +- Input: + - Fully integrated to work with physiological logfiles for Philips MR systems (SCANPHYSLOG) + - tested for General Electric (GE) log-files + - implementation for Siemens log-files + - also: interface for 'Custom', i.e. general heart-beat time stamps + & breathing volume time courses from other log formats +- Output: + - Nuisance regressors for mass-univariate statistical analysis with SPM5,8,12 + or as text file for export to any other package + - raw and processed physiological logfile data +- Part of the TAPAS Software Collection of the Translational Neuromodeling Unit (TNU) Zurich:long term support and ongoing development + + +Contributors +------------ + +- Lead Programmer: + - Lars Kasper, TNU & MR-Technology Group, IBT, University of Zurich & ETH Zurich +- Project Team: + - Steffen Bollmann, Centre for Advanced Imaging, University of Queensland, Australia + - Saskia Bollmann, Centre for Advanced Imaging, University of Queensland, Australia +- Contributors: + - Eduardo Aponte, TNU Zurich + - Tobias U. Hauser, FIL London, UK + - Jakob Heinzle, TNU Zurich + - Chloe Hutton, FIL London, UK (previously) + - Miriam Sebold, Charite Berlin, Germany + + +Contact +------- +Send bug reports and suggestions either to +1) our mailing list: tapas@sympa.ethz.ch, or +2) as an issue on our TAPAS github account : https://github.com/translationalneuromodeling/tapas/issues + + +References +---------- + +### Main Toolbox Reference ### +1. Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. Journal of Neuroscience Methods 276, 56–72. doi:10.1016/j.jneumeth.2016.10.019 + +### Related Papers (Implemented noise correction algorithms) ### +2. Glover, G.H., Li, T.Q. & Ress, D. Image‐based method for retrospective correction +of PhysIOlogical motion effects in fMRI: RETROICOR. Magn Reson Med 44, 162-7 (2000). + +3. Hutton, C. et al. The impact of PhysIOlogical noise correction on fMRI at 7 T. +NeuroImage 57, 101‐112 (2011). + +4. Harvey, A.K. et al. Brainstem functional magnetic resonance imaging: +Disentangling signal from PhysIOlogical noise. Journal of Magnetic Resonance +Imaging 28, 1337‐1344 (2008). + +5. Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise +correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, +90–101. doi:10.1016/j.neuroimage.2007.04.042 + +6. Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. The respiration response +function: The temporal dynamics of fMRI s ignal fluctuations related to changes in +respiration. NeuroImage 40, 644–654. doi:10.1016/j.neuroimage.2007.11.059 + +7. Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate on the +BOLD signal: The cardiac response function. NeuroImage 44, 857–869. +doi:10.1016/j.neuroimage.2008.09.029 + +8. Siegel, J.S., Power, J.D., Dubis, J.W., Vogel, A.C., Church, J.A., Schlaggar, B.L., +Petersen, S.E., 2014. Statistical improvements in functional magnetic resonance +imaging analyses produced by censoring high-motion data points. Hum. Brain Mapp. +35, 1981–1996. doi:10.1002/hbm.22307 diff --git a/PhysIO/README.txt b/PhysIO/README.txt deleted file mode 100644 index 901015e1..00000000 --- a/PhysIO/README.txt +++ /dev/null @@ -1,230 +0,0 @@ -TAPAS PhysIO Toolbox Version 2016 - -************************************************************************ -Copyright (C) 2012-2016 Lars Kasper -Translational Neuromodeling Unit (TNU) -Institute for Biomedical Engineering -University of Zurich and ETH Zurich ------------------------------------------------------------------------- - -The PhysIO Toolbox is free software: you can redistribute it and/or -modify it under the terms of the GNU General Public License as -published by the Free Software Foundation, either version 3 of the -License, or (at your option) any later version. - -This program is distributed in the hope that it will be useful, but -WITHOUT ANY WARRANTY; without even the implied warranty of -MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -General Public License for more details. - -You should have received a copy of the GNU General Public License -along with this program (see the file COPYING). If not, see -. -************************************************************************ - -PURPOSE - -The PhysIO Toolbox provides model-based physiological noise correction of -fMRI data using peripheral measures of respiration and cardiac pulsation. -It incorporates noise models of cardiac/respiratory phase (RETROICOR, -Glover et al. 2000), as well as heart rate variability and respiratory -volume per time (cardiac response function, Chang et. al, 2009, respiratory -response function, Birn et al. 2006). The toolbox is usable via the SPM -batch editor, performs automatic pre-processing of noisy peripheral data -and outputs nuisance regressor files directly suitable for SPM -(multiple_regressors.txt). - -BACKGROUND - -The PhysIO Toolbox provides physiological noise correction for fMRI-data -from peripheral measures (ECG/pulse oximetry, breathing belt). It is -model-based, i.e. creates nuisance regressors from the physiological -monitoring that can enter a General Linear Model (GLM) analysis, e.g. -SPM8/12. Furthermore, for scanner vendor logfiles (PHILIPS, GE, Siemens), -it provides means to statistically assess peripheral data (e.g. heart rate variability) -and recover imperfect measures (e.g. distorted R-peaks of the ECG). - -Facts about physiological noise in fMRI: -- Physiological noise can explain 20-60 % of variance in fMRI voxel time - series (Birn2006, Hutton2011, Harvey2008. - - Physiological noise affects a lot of brain regions (s. figure, e.g. - brainstem or OFC), especially next to CSF, arteries (Hutton2011). - - If not accounted for, this is a key factor limiting sensitivity for - effects of interest. -- Physiological noise contributions increase with field strength; they - become a particular concern at and above 3 Tesla (Kasper2009, Hutton2011). -- In resting state fMRI, disregarding physiological noise leads to wrong - connectivity results (Birn2006). - -=> Some kind of physiological noise correction is highly recommended for - every statistical fMRI analysis. - -Model-based correction of physiological noise: -- Physiological noise can be decomposed into periodic time series following - heart rate and breathing cycle. -- The Fourier expansion of cardiac and respiratory phases was introduced as - RETROICOR (RETROspective Image CORrection, Glover2000, - see also Josephs1997). -- These Fourier Terms can enter a General Linear Model (GLM) as nuisance - regressors, analogous to movement parameters. -- As the physiological noise regressors augment the GLM and explain - variance in the time series, they increase sensitivity in all contrasts - of interest. - - -FEATURES OF THIS TOOLBOX - -Physiological Noise Modeling : -- Modeling physiological noise regressors from peripheral data - (breathing belt, ECG, pulse oximeter) -- State of the art RETROICOR cardiac and respiratory phase expansion -- Cardiac response function (Chang et al, 2009) and respiratory response - function (Birn et al. 2006) modelling of heart-rate variability and - respiratory volume per time influence on physiological noise -- Flexible expansion orders to model different contributions of cardiac, - respiratory and interaction terms (see Harvey2008, Hutton2011) -- Automatic creation of nuisance regressors, full integration into standard - GLMs, tested for SPM8/12 ("multiple_regressors.mat") - -Flexible Read-in: -The toolbox is dedicated to seamless integration into a clinical research s -etting and therefore offers correction methods to recover physiological -data from imperfect peripheral measures. -- General Electric logfiles -- Siemens logfiles (synchronization via DICOM time stamp or tics time scale) -- Philips SCANPHYSLOG-file handling, including automatic alignment of scan volume timing and physiological time series - from logged gradient timecourses -- BioPac .mat-files - - -COMPATIBILITY & SUPPORT - -- Matlab Toolbox -- Input: - - Fully integrated to work with physiological logfiles for Philips MR systems (SCANPHYSLOG) - - tested for General Electric (GE) log-files - - implementation for Siemens log-files - - also: interface for 'Custom', i.e. general heart-beat time stamps - & breathing volume time courses from other log formats -- Output: - - Nuisance regressors for mass-univariate statistical analysis with SPM5,8,12 - or as text file for export to any other package - - raw and processed physiological logfile data -- Part of the TNU Software Edition: long term support and ongoing development - - -DOWNLOADS & RELEASE INFORMATION - -- Current Release: - -PhysIO_Toolbox_R2017.1 - -February 19, 2017 - -Minor Release Notes (R2017.1) -* Substantially improved Siemens interface, both for VB/VD and 3T/7T releases -** several bugfixes -** based on extensive user feedback from Berlin and Brisbane -* New functionality tapas_physio_overlay_contrasts.m to display non-physio - contrasts automatically as well - -Major Release Notes (r904 / R2016.1): -- Software version for accepted PhysIO Toolbox Paper: doi:10.1016/j.jneumeth.2016.10.019 -- Tested and expanded versions of examples -- Improved stability by bugfixes and compatibility to Matlab R2016 -- Slice-wise regressor creation -- Detection of constant physiological time series (detachment, clipping) -- Refactoring of report_contrasts and compute_tsnr_gains as standalone functionality -- Improved Read-in capabilities (Siemens respiration data, BioPac .mat) -- Migration from svn (r904) to git (tnurepository) for version control - -Major Release Notes (r835): -- Software version for Toolbox Paper submission -- Noise ROIs modeling -- Extended motion models (24 parameters, Volterra expansion) -- HRV/RVT models with optional multiple delay regressors -- Report_contrasts with automatic contrast generation for all regressor groups -- compute_tsnr_gains for individual physiological regressor groups -- consistent module naming (scan_timing, preproc) -- Visualisation improvement (color schemes, legends) - -Minor Release Notes (r666): -- Compatibility tested for SPM12, small bugfixes Batch Dependencies -- Cleaner Batch Interface with grouped sub-menus (cfg_choice) -- new model: 'none' to just read out physiological raw data and preprocess, - without noise modelling -- Philips: Scan-timing via gradient log now automatized (gradient_log_auto) -- Siemens: Tics-Logfile read-in (proprietary, needs Siemens-agreement) -- All peak detections (cardiac/respiratory) now via auto_matched algorithm -- Adapt plots/saving for Matlab R2014b - -Major Release Notes (r534): -- Read-in of Siemens plain text log files; new example dataset for Siemens -- Speed up and debugging of auto-detection method for noisy cardiac data => new method thresh.cardiac.initial_cpulse_select.method = ???auto_matched??? -- Error handling for temporary breathing belt failures (Eduardo Aponte, TNU Zurich) -- slice-wise regressors can be created by setting sqpar.onset_slice to a index vector of slices - -Major Release Notes (r497): -- SPM matlabbatch GUI implemented (Call via Batch -> SPM -> Tools -> TAPAS PhysIO Toolbox) -- improved, automatic heartbeat detection for noisy ECG now standard for ECG and Pulse oximetry (courtesy of Steffen Bollmann) -- QuickStart-Manual and PhysIO-Background presentation expanded/updated -- job .m/.mat-files created for all example datasets -- bugfixes cpulse-initial-select method-handling (auto/manual/load) - -Major Release Notes (r429): -- Cardiac and Respiratory response function regressors integrated in workflow (heart rate and breathing volume computation) -- Handling of Cardiac and Respiratory Logfiles only -- expanded documentation (Quickstart.pdf and Handbook.pdf) -- read-in of custom log files, e.g. for BrainVoyager peripheral data -- more informative plots and commenting (especially in tapas_physio_new). - -Minor Release Notes (r354): -- computation of heart and breathing rate in Philips/PPU/main_PPU.m -- prefix of functions with tapas_* - -Major Release Notes (r241): -- complete modularization of reading/preprocessing/regressor creation for peripheral physiological data -- manual selection of missed heartbeats in ECG/pulse oximetry (courtesy of Jakob Heinzle) -- support for logfiles from GE scanners (courtesy of Steffen Bollmann, KiSpi Zuerich) -- improved detection of pulse oximetry peaks (courtesy of Steffen Bollmann) -- improved documentation -- consistent function names (prefixed by "physio_") - -NOTE: Your main_ECG/PPU.m etc. scripts from previous versions (<=r159) will not work with this one any more. Please adapt one of the example scripts for your needs (~5 min of work). The main benefit of this version is a complete new variable structure that is more sustainable and makes the code more readable. - - -Lead Programmer: Lars Kasper, TNU & MR-Technology Group, IBT, University & ETH Zurich - -Contributors: -Steffen Bollmann, Children's Hospital Zurich & ETH Zurich -Jakob Heinzle, TNU Zurich -Eduardo Aponte, TNU Zurich - -Send bug reports and suggestions to: kasper@biomed.ee.ethz.ch - - -TUTORIAL - -Run main_ECG3T.m in subdirectory "examples" of the toolbox -See subdirectory "manual" - - -REFERENCES - -Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2016. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. Journal of Neuroscience Methods accepted. doi:10.1016/j.jneumeth.2016.10.019 - -Birn, Rasmus M., Jason B. Diamond, Monica A. Smith, and Peter A. Bandettini. 2006. Separating Respiratory-variation-related Fluctuations from Neuronal-activity-related Fluctuations in fMRI. NeuroImage 31 (4) (July 15): 1536?1548. doi:10.1016/j.neuroimage.2006.02.048. - -Glover, G H, T Q Li, and D Ress. 2000. Image-based Method for Retrospective Correction of Physiological Motion Effects in fMRI: RETROICOR. Magnetic Resonance in Medicine: Official Journal of the Society of Magnetic Resonance in Medicine 44 (1) (July): 162(7). doi:10893535. - -Harvey, Ann K., Kyle T.S. Pattinson, Jonathan C.W. Brooks, Stephen D. Mayhew, Mark Jenkinson, and Richard G. Wise. 2008. Brainstem Functional Magnetic Resonance Imaging: Disentangling Signal from Physiological Noise. Journal of Magnetic Resonance Imaging 28 (6): 1337?1344. doi:10.1002/jmri.21623. - -Hutton, C., O. Josephs, J. Stadler, E. Featherstone, A. Reid, O. Speck, J. Bernarding, and N. Weiskopf. 2011. The Impact of Physiological Noise Correction on fMRI at 7 T. NeuroImage 57 (1) (July 1): 101?112. doi:10.1016/j.neuroimage.2011.04.018. - -Josephs, O., Howseman, A.M., Friston, K., Turner, R., 1997. Physiological noise modelling for multi-slice EPI fMRI using SPM. Proceedings of the 5th Annual Meeting of ISMRM, Vancouver, Canada, p. 1682 - -Kasper, Lars, Sarah Marti, S. Johanna Vannesjo, Chloe Hutton, Ray Dolan, Nikolaus Weiskopf, Klaas Enno Stephan, and Klaas Paul Pruessmann. 2009. Cardiac Artefact Correction for Human Brainstem fMRI at 7 Tesla. In Proc. Org. Hum. Brain Mapping 15, 395. San Francisco. - - -VERSION OF THIS FILE -$Id$ diff --git a/PhysIO/code/tapas_physio_cfg_matlabbatch.m b/PhysIO/code/tapas_physio_cfg_matlabbatch.m index 69e618aa..f77995de 100644 --- a/PhysIO/code/tapas_physio_cfg_matlabbatch.m +++ b/PhysIO/code/tapas_physio_cfg_matlabbatch.m @@ -48,6 +48,7 @@ ' extra acquisition (scan_timing) logfile with' ' time stamps of all volumes and slices' ' ''Biopac_Mat'' - exported mat files from Biopac system' + ' ''BrainProducts'' - .eeg files from BrainProducts EEG system' ' ' ' or' ' ''Custom''' @@ -69,8 +70,8 @@ ' NOTE: the sampling interval has to be specified for these files as' ' well (s.b.)' }; -vendor.labels = {'Philips', 'GE', 'Siemens (VB, *.puls/*.ecg/*.resp)', 'Siemens_Tics (VD: *_PULS.log/*_ECG1.log/*_RESP.log/*_AcquisitionInfo*.log)', 'Biopac_Mat', 'Custom'}; -vendor.values = {'Philips', 'GE', 'Siemens', 'Siemens_Tics', 'Biopac_Mat', 'Custom'}; +vendor.labels = {'Philips', 'GE', 'Siemens (VB, *.puls/*.ecg/*.resp)', 'Siemens_Tics (VD: *_PULS.log/*_ECG1.log/*_RESP.log/*_AcquisitionInfo*.log)', 'Biopac_Mat', 'BrainProducts', 'Custom'}; +vendor.values = {'Philips', 'GE', 'Siemens', 'Siemens_Tics', 'Biopac_Mat', 'BrainProducts', 'Custom'}; vendor.val = {'Philips'}; %-------------------------------------------------------------------------- diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m index f6aeeaa6..be034a20 100644 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m +++ b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_auto_philips.m @@ -23,6 +23,7 @@ % 3. Determine slice events between all detected volumes % (again, creating a slice template and matching it to time series % between consecutive volumes...) +% 4. Correction for incomplete volumes (not enough slices/volume) % % IN % log_files is a structure containing the following filenames (with full @@ -253,10 +254,10 @@ end if debug - % verbose.fig_handles(end+1) = plot_slice_events( LOCS, t, ... - % gradient_choice, templateGradientSlice, secondGuessLOCS); - % - % plot_diff_LOCS(t, LOCS, dt) + verbose.fig_handles(end+1) = plot_slice_events( LOCS, t, ... + gradient_choice, templateGradientSlice, secondGuessLOCS); + + plot_diff_LOCS(t, LOCS, dt) end @@ -299,7 +300,8 @@ % remove erroneous volume events, i.e. those due to slice gaps that are % assumed to be volume end gaps, and create new volumes with not enough % slices -minVolumeDistanceSamplesError = ceil((1-1/sqpar.Nslices) * sqpar.TR/dt); +nDeltaSlicesAllowed = 2; +minVolumeDistanceSamplesError = ceil((1-nDeltaSlicesAllowed/sqpar.Nslices) * sqpar.TR/dt); idxVolError = find(diff(VOLLOCS) < minVolumeDistanceSamplesError); if ~isempty(idxVolError) VOLLOCS(idxVolError(2:2:end)) = []; diff --git a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m index e47145ec..bd0e7900 100644 --- a/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m +++ b/PhysIO/code/tapas_physio_create_scan_timing_from_gradients_philips.m @@ -259,15 +259,15 @@ xlabel('t(s)'); % Plot gradient thresholding for slice timing determination - - if ~isempty(VOLLOCS) - hp(end+1) = stem(t(VOLLOCS), 1.25*ones(size(VOLLOCS))); hold all + heightStem = median(gradient_choice(VOLLOCS)); + hp(end+1) = stem(t(VOLLOCS), 1.25*heightStem*ones(size(VOLLOCS))); hold all lg{end+1} = sprintf('Found volume events (N = %d)', numel(VOLLOCS)); end if ~isempty(LOCS) - hp(end+1) = stem(t(LOCS), ones(size(LOCS))); hold all + heightStem = median(gradient_choice(LOCS)); + hp(end+1) = stem(t(LOCS), heightStem*ones(size(LOCS))); hold all lg{end+1} = sprintf('Found slice events (N = %d)', numel(LOCS)); dLocsSecs = diff(LOCS)*dt*1000; diff --git a/PhysIO/code/tapas_physio_example_philips_ecg3t_v2_spm_job.m b/PhysIO/code/tapas_physio_example_philips_ecg3t_v2_spm_job.m new file mode 100644 index 00000000..bb8f32e1 --- /dev/null +++ b/PhysIO/code/tapas_physio_example_philips_ecg3t_v2_spm_job.m @@ -0,0 +1,43 @@ +%----------------------------------------------------------------------- +% Job saved on 24-Apr-2017 15:26:29 by cfg_util (rev $Rev: 6460 $) +% spm SPM - SPM12 (6906) +% cfg_basicio BasicIO - Unknown +%----------------------------------------------------------------------- +matlabbatch{1}.spm.tools.physio.save_dir = {'/Users/kasperla/Documents/code/matlab/PhysIO/private/examples/Philips/ECG3T_V2'}; +matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Philips'; +matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'/Users/kasperla/Documents/code/matlab/PhysIO/private/examples/Philips/ECG3T_V2/SCANPHYSLOG.log'}; +matlabbatch{1}.spm.tools.physio.log_files.respiration = {'/Users/kasperla/Documents/code/matlab/PhysIO/private/examples/Philips/ECG3T_V2/SCANPHYSLOG.log'}; +matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''}; +matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = []; +matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0; +matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last'; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 32; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = []; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 2.5; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 5; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 305; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 16; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = []; +matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = []; +matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]); +matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU_Wifi'; +matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4; +matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat'; +matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]); +matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt'; +matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat'; +matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none'; +matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3; +matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4; +matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1; +matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]); +matlabbatch{1}.spm.tools.physio.model.hrv.no = struct([]); +matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]); +matlabbatch{1}.spm.tools.physio.model.movement.yes.file_realignment_parameters = {''}; +matlabbatch{1}.spm.tools.physio.model.movement.yes.order = 6; +matlabbatch{1}.spm.tools.physio.model.movement.yes.outlier_translation_mm = Inf; +matlabbatch{1}.spm.tools.physio.model.movement.yes.outlier_rotation_deg = Inf; +matlabbatch{1}.spm.tools.physio.model.other.no = struct([]); +matlabbatch{1}.spm.tools.physio.verbose.level = 2; +matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = ''; +matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false; diff --git a/PhysIO/code/tapas_physio_fill_empty_parameters.m b/PhysIO/code/tapas_physio_fill_empty_parameters.m index aea543b5..05ca6b0f 100644 --- a/PhysIO/code/tapas_physio_fill_empty_parameters.m +++ b/PhysIO/code/tapas_physio_fill_empty_parameters.m @@ -43,6 +43,11 @@ physio.log_files.scan_timing = {''}; end +if strcmpi(physio.preproc.cardiac.initial_cpulse_select.method, 'auto_matched') && ... + isempty(physio.preproc.cardiac.initial_cpulse_select.min) + physio.preproc.cardiac.initial_cpulse_select.min = 0.4; +end + if isempty(physio.log_files.sampling_interval) switch lower(physio.log_files.vendor) case 'philips' @@ -55,7 +60,7 @@ end case 'ge' physio.log_files.sampling_interval = 25e-3; - case {'biopac_mat', 'siemens', 'siemens_tics'} % will be read from file later + case {'biopac_mat', 'siemens', 'siemens_tics', 'brainproducts'} % will be read from file later physio.log_files.sampling_interval = []; otherwise % e.g. custom error('Please specify sampling interval for custom text data'); diff --git a/PhysIO/code/tapas_physio_findpeaks.m b/PhysIO/code/tapas_physio_findpeaks.m index 7b4bb5fa..c8cc5e7c 100644 --- a/PhysIO/code/tapas_physio_findpeaks.m +++ b/PhysIO/code/tapas_physio_findpeaks.m @@ -70,12 +70,24 @@ end %#function dspopts.findpeaks -hopts = tapas_physio_uddpvparse('dspopts.findpeaks',varargin{:}); -Ph = hopts.MinPeakHeight; -Pd = hopts.MinPeakDistance; -Th = hopts.Threshold; -Np = hopts.NPeaks; -Str = hopts.SortStr; +%hopts = tapas_physio_uddpvparse('dspopts.findpeaks',varargin{:}); +if nargin + varargin(1:2:end) = lower(varargin(1:2:end)); +end + +defaults.minpeakheight = -Inf; +defaults.minpeakdistance = []; +defaults.threshold = 0; +defaults.npeaks = []; +defaults.sortstr = 'none'; + +hopts = tapas_physio_propval(varargin, defaults); + +Ph = hopts.minpeakheight; +Pd = hopts.minpeakdistance; +Th = hopts.threshold; +Np = hopts.npeaks; +Str = hopts.sortstr; % Validate MinPeakDistance if ~isempty(Pd) && (~isnumeric(Pd) || ~isscalar(Pd) ||any(rem(Pd,1)) || (Pd < 1)) diff --git a/PhysIO/code/tapas_physio_main_create_regressors.m b/PhysIO/code/tapas_physio_main_create_regressors.m index 79b4e0ca..3a66d866 100644 --- a/PhysIO/code/tapas_physio_main_create_regressors.m +++ b/PhysIO/code/tapas_physio_main_create_regressors.m @@ -200,7 +200,7 @@ minConstantIntervalAlertSamples); end - [ons_secs, scan_timging.sqpar, verbose] = tapas_physio_crop_scanphysevents_to_acq_window(... + [ons_secs, scan_timing.sqpar, verbose] = tapas_physio_crop_scanphysevents_to_acq_window(... ons_secs, scan_timing.sqpar, verbose); sqpar = scan_timing.sqpar; diff --git a/PhysIO/code/tapas_physio_plot_cropped_phys_to_acqwindow.m b/PhysIO/code/tapas_physio_plot_cropped_phys_to_acqwindow.m index dafa0aa4..25a7e443 100644 --- a/PhysIO/code/tapas_physio_plot_cropped_phys_to_acqwindow.m +++ b/PhysIO/code/tapas_physio_plot_cropped_phys_to_acqwindow.m @@ -121,7 +121,7 @@ ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ... ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ... ['cardiac pulse (heartbeat) marker (N = ' int2str(length(cpulse)) ')'], ... - 'cardiac signal', 'respiratory signal'}); + 'cardiac signal (dashed = raw)', 'respiratory signal (dashed = raw)'}); elseif hasCardiacData legend( hs, { 'used cardiac signal', ... @@ -129,14 +129,14 @@ ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ... ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ... ['cardiac pulse (heartbeat) marker (N = ' int2str(length(cpulse)) ')'], ... - 'cardiac signal'}); + 'cardiac signal (dashed = raw)'}); else % only respData legend( hs, { 'used respiratory signal', ... ['dummy scan event marker (N = ' int2str(Ndummies*Nslices) ')'], ... ['volume event marker (N = ' int2str(length(svolpulse)-Ndummies) '), without dummies'], ... ['scan event marker (N = ' int2str(length(spulse)-Ndummies*Nslices) ')'], ... - 'respiratory signal'}); + 'respiratory signal (dashed = raw)'}); end ylim(1.4*maxVal*[-1 1]); diff --git a/PhysIO/code/tapas_physio_plot_gradient.m b/PhysIO/code/tapas_physio_plot_gradient.m index bd9a2489..6302507c 100644 --- a/PhysIO/code/tapas_physio_plot_gradient.m +++ b/PhysIO/code/tapas_physio_plot_gradient.m @@ -46,7 +46,7 @@ 'G_x' 'G_y' 'G_z' - '|G|_2' + '|G|^2' }; for m = 1:4 diff --git a/PhysIO/code/tapas_physio_read_physlogfiles.m b/PhysIO/code/tapas_physio_read_physlogfiles.m index 7e70c4ca..ed389c90 100644 --- a/PhysIO/code/tapas_physio_read_physlogfiles.m +++ b/PhysIO/code/tapas_physio_read_physlogfiles.m @@ -56,13 +56,23 @@ end switch lower(log_files.vendor) - case 'philips' + case 'biopac_mat' [c, r, t, cpulse, acq_codes] = ... - tapas_physio_read_physlogfiles_philips(log_files, cardiac_modality); + tapas_physio_read_physlogfiles_biopac_mat(log_files, cardiac_modality, verbose); + case 'custom' + [c, r, t, cpulse] = ... + tapas_physio_read_physlogfiles_custom(log_files, verbose); + acq_codes = []; + case 'brainproducts' + [c, r, t, cpulse, acq_codes] = ... + tapas_physio_read_physlogfiles_brainproducts(log_files, cardiac_modality, verbose); case 'ge' [c, r, t, cpulse] = ... tapas_physio_read_physlogfiles_GE(log_files, verbose); acq_codes = []; + case 'philips' + [c, r, t, cpulse, acq_codes] = ... + tapas_physio_read_physlogfiles_philips(log_files, cardiac_modality); case 'siemens' [c, r, t, cpulse, verbose] = ... tapas_physio_read_physlogfiles_siemens(log_files, verbose); @@ -71,13 +81,6 @@ [c, r, t, cpulse, verbose] = ... tapas_physio_read_physlogfiles_siemens_tics(log_files, verbose); acq_codes = []; - case 'biopac_mat' - [c, r, t, cpulse, acq_codes] = ... - tapas_physio_read_physlogfiles_biopac_mat(log_files, cardiac_modality, verbose); - case 'custom' - [c, r, t, cpulse] = ... - tapas_physio_read_physlogfiles_custom(log_files, verbose); - acq_codes = []; end % Do not prepend for Siemens Tics, since can be as long as a day diff --git a/PhysIO/code/tapas_physio_read_physlogfiles_brainproducts.m b/PhysIO/code/tapas_physio_read_physlogfiles_brainproducts.m new file mode 100644 index 00000000..cdb275c4 --- /dev/null +++ b/PhysIO/code/tapas_physio_read_physlogfiles_brainproducts.m @@ -0,0 +1,183 @@ +function [c, r, t, cpulse, acq_codes] = tapas_physio_read_physlogfiles_brainproducts(log_files, ... + cardiac_modality, verbose) +% reads out physiological time series (ECG, PMU, resp belt) and timing vector for BrainProducts .eeg file +% +% [cpulse, rpulse, t, c, acq_codes] = tapas_physio_read_physlogfiles_brainproducts(logfiles, ... +% verbose) +% +% +% IN log_files +% .cardiac contains ECG or pulse oximeter time course +% for BrainProducts: usually the same as respiration +% .respiration contains breathing belt amplitude time course +% for BrainProducts: usually the same as cardiac +% .sampling_interval is ignored here, read from logfile +% +% cardiac_modality +% 'ecg1_filtered' filtered 1st ECG channel signal +% (Default) +% 'ecg2_filtered' filteered 2nd ECG channel +% (sometimes less gradient artifacts) +% 'ecg1_raw' raw 1st ECG channel +% +% verbose +% .level debugging plots are created if level >=3 +% .fig_handles appended by handle to output figure +% +% OUT +% r respiratory time series +% c cardiac time series (ECG or pulse oximetry) +% t vector of time points (in seconds) +% cpulse time events of R-wave peak in cardiac time series (seconds) +% for Biopac: usually empty, kept for compatibility +% acq_codes slice/volume start events marked by number <> 0 +% for time points in t +% EXAMPLE +% [ons_secs.cpulse, ons_secs.rpulse, ons_secs.t, ons_secs.c] = +% tapas_physio_read_physlogfiles_GE(logfiles); +% +% See also tapas_physio_main_create_regressors +% +% Author: Lars Kasper +% Created: 2017-01-27 +% Copyright (C) 2017 Institute for Biomedical Engineering, ETH/Uni Zurich. +% +% This file is part of the PhysIO toolbox, which is released under the terms of the GNU General Public +% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL +% (either version 3 or, at your option, any later version). For further details, see the file +% COPYING or . + + +%% user input + +% volume trigger event +event_type = 'Response'; +event_value = 'R128'; + +switch lower(cardiac_modality) + case {'ecg','ecg1_filtered'} + % chose one ecg channel (only for visualisation) + ecg_ch = 1; + case 'ecg2_filtered' + % chose one ecg channel (only for visualisation) + ecg_ch = 2; +end + +% is ecg data flipped? +ecg_is_flipped = 1; + +%% data + +% read data and header info using field trip (included in SPM) +if ~exist('ft_read_header', 'file') + pathSpm = fileparts(which('spm')); + addpath(genpath(fullfile(pathSpm, 'external', 'fieldtrip'))); +end + +hdr = ft_read_header(log_files.cardiac); +data = ft_read_data(log_files.cardiac); + +fs = hdr.Fs; % sampling frequency in Hz +N = hdr.nSamples; % number of samples +dt = 1/fs; % sampling interval in seconds +t = linspace(0,dt*(N-1), N); % time vector in seconds + +fh = []; + +% plot first 10 seconds of the raw data +fh(end+1,1) = tapas_physio_get_default_fig_params(); +plot_end = 10 * fs; +plot(t(1:plot_end), data(:,(1:plot_end))); +xlabel('time in seconds'); +ylabel(hdr.chanunit); +legend(hdr.label); + +% extract ECG data +if ecg_is_flipped + s = -data(ecg_ch,:); +else + s = data(ecg_ch,:); +end + +% plot first 10 seconds again for check +fh(end+1,1) = tapas_physio_get_default_fig_params(); +plot_end = 10 * fs; +plot(t(1:plot_end), s(:,(1:plot_end))); +xlabel('time in seconds'); +ylabel(hdr.chanunit{ecg_ch}); +legend(hdr.label{ecg_ch}); +%% events +% display all events in the data on the command line +cfg = []; +cfg.dataset = log_files.cardiac; +cfg.trialdef.eventtype = '?'; +ft_definetrial(cfg); + +% define volume trigger +cfg.trialdef.eventtype = event_type; +cfg.trialdef.eventvalue = event_value; +cfg = ft_definetrial(cfg); +trigger_pos = cfg.trl(:,1); +n_trigger = length(trigger_pos); + +% plot raw data and events +fh(end+1,1) = tapas_physio_get_default_fig_params(); +max_value = max(s); +plot(t, s); +hold all; +stem(t(trigger_pos), max_value * ones(1,n_trigger)); + +% find segments +% compute the number of samples between each trigger (pos-1) +diff_trigger_pos = diff(trigger_pos); +% find the positons where a change is happening +diff_diff_trigger_pos = find(diff(diff_trigger_pos)); +% add first and last event +start_segment = [1; diff_diff_trigger_pos(2:2:end)+1]; +end_segment = [diff_diff_trigger_pos(1:2:end)+1; n_trigger]; +% number of samples per trials +n_trial_samples = diff_trigger_pos(start_segment); +% position of start segment in samples +sample_start_segment = trigger_pos(start_segment); +% position of end segment in samples (not including the whole trial) +sample_end_segment = trigger_pos(end_segment); +% number of segments +n_segments = length(start_segment); +% number of trials per segment +n_trials = end_segment - start_segment + 1; +% TR per trial +TR_trials = n_trial_samples*dt; + +% plot segments +stem(t(sample_start_segment), max_value * ones(1,n_segments), '--', 'LineWidth', 5); +stem(t(sample_end_segment), max_value * ones(1,n_segments), '--', 'LineWidth', 5); +legend('signal', 'all trigger', 'start segment', 'end segment'); + +%% load only first segment for sanity check +cfg = []; +cfg.dataset = log_files.cardiac; +cfg.trialdef.eventtype = event_type; +cfg.trialdef.eventvalue = event_value; +cfg.trialfun = 'ft_trialfun_segment'; +cfg.trialdef.prestim = 0; +cfg.trialdef.poststim = TR_trials(1); +cfg.trialdef.segment_start = start_segment(1); +cfg.trialdef.segment_end = end_segment(1); +cfg = ft_definetrial(cfg); + +% plot segments +fh(end+1,1) = tapas_physio_get_default_fig_params(); +max_value = max(s); +plot(t, s); +hold all; +stem(t(cfg.trl(:,1)), max_value * ones(1,n_trials(1))); +stem(t(cfg.trl(:,2)), max_value * ones(1,n_trials(1))); +legend('signal', 'start trial', 'end trial'); + +%% save results +[~,name] = fileparts(file_name); +save_name = fullfile(file_path, [name, '_segments.mat']); +save(save_name, 'sample_start_segment', 'sample_end_segment', ... + 'n_segments', 'n_trials', 'TR_trials', 'start_segment', 'end_segment', ... + 'event_type', 'event_value', 'n_trial_samples', 'fs', ... + 'ecg_is_flipped', 'ecg_ch', 'file_path', 'file_name'); \ No newline at end of file diff --git a/PhysIO/code/tapas_physio_read_physlogfiles_philips_matrix.m b/PhysIO/code/tapas_physio_read_physlogfiles_philips_matrix.m index 1e0fe2af..659dddb4 100644 --- a/PhysIO/code/tapas_physio_read_physlogfiles_philips_matrix.m +++ b/PhysIO/code/tapas_physio_read_physlogfiles_philips_matrix.m @@ -1,13 +1,19 @@ function y = tapas_physio_read_physlogfiles_philips_matrix(fileNamePhyslog) -% reads out complete [nSamples,10] numerical matrix from SCANPHYSLOG-file +% reads out complete numerical matrix from SCANPHYSLOG-file +% Physlog file version 1: [nSamples,10] +% Physlog file version 2: [nSamples,11] +% % IN % fileNamePhyslog % string, filename of SCANPHYSLOG-file % OUT -% y [nSamples, 10] SCANPHYSLOG-file matrix with the following +% y [nSamples, 10 or 11] SCANPHYSLOG-file matrix with the following % columns: +% Version 1: % v1raw v2raw v2 v2 PPU Pneumatic_belt Gradient_x Gradient_y % Gradient_z acq_codes(scan/ECG trigger etc.) +% Version 2: +% v1raw v2raw v1 v2 ppu resp gx gy gz mark mark2 % % EXAMPLE % y = tapas_physio_read_physlogfiles_philips_matrix('SCANPHYSLOG.log') @@ -28,19 +34,47 @@ % $Id: tapas_physio_read_physlogfiles_philips_matrix.m 632 2015-01-09 12:36:12Z kasperla $ % use textread as long as it exists, for it is much faster (factor 4) than -% textscan; TODO: use fread ans sscanf to make it even faster... +% textscan; TODO: use fread and sscanf to make it even faster... + +%% check log file version by header line +fid = fopen(fileNamePhyslog, 'r'); +stringFirstLine = fgetl(fid); +fclose(fid); + +if any(strfind(stringFirstLine, 'Physlog file version = 2')) + versionLogfile = 2; +else + versionLogfile = 1; +end + +switch versionLogfile + case 1 + parseString = '%d %d %d %d %d %d %d %d %d %s'; + nColumns = 10; + case 2 + parseString = '%d %d %d %d %d %d %d %d %d %s %s'; + nColumns = 11; +end + + +%% Read decimal value colums if exist('textread') - [z{1:10}] = textread(fileNamePhyslog,'%d %d %d %d %d %d %d %d %d %s','commentstyle', 'shell'); + [z{1:nColumns}] = textread(fileNamePhyslog, parseString,'commentstyle', 'shell'); else fid = fopen(fileNamePhyslog, 'r'); - z = textscan(fid, '%d %d %d %d %d %d %d %d %d %s', 'commentstyle', '#'); + z = textscan(fid, parseString, 'commentstyle', '#'); z(1:9) = cellfun(@double, z(1:9), 'UniformOutput', false); fclose(fid); end -z{10} = hex2dec(z{10}); % hexadecimal acquisition codes converted; -% account for incomplete rows +%% Convert hexadecimal acquisition codes +for iCol = 10:nColumns + z{iCol} = hex2dec(z{iCol}); +end + + +%% Account for incomplete rows nMinRows = min(cellfun(@numel,z)); z = cellfun(@(x) x(1:nMinRows), z, 'UniformOutput', false); y = cell2mat(z); diff --git a/PhysIO/code/tapas_physio_rescale_gradient_gain_fluctuations.m b/PhysIO/code/tapas_physio_rescale_gradient_gain_fluctuations.m index b770195e..019202c2 100644 --- a/PhysIO/code/tapas_physio_rescale_gradient_gain_fluctuations.m +++ b/PhysIO/code/tapas_physio_rescale_gradient_gain_fluctuations.m @@ -13,7 +13,7 @@ % max 1 after rescaling; % if false: gain of last interval chosen for all other % OUT -% G normlized gain gradient time course +% G normalized gain gradient time course % gainArray [nSteps+1,1] vector of detected gains (i.e. maxima in % multiple TR interval) in gradient time course % normFactor max value to which all intervals are scaled in this @@ -42,7 +42,7 @@ defaults.doNormalize = true; % since gain changes happen for whole slice-waveform, look where slice -% waveform ends, i.e. a bit earlier for step-up, and a bit later for +% waveform ends, i.e., a bit earlier for step-up, and a bit later for % step-down events of gain defaults.doExtendStepToSliceRange = true; @@ -51,7 +51,7 @@ minPeakHeight = 1000; % TODO: remove this heuristic! ignoreBoundaryPercent = 30; % for gain estimation in interval margin - % is ignored in case of a slow change + % is ignored in case of a slow change normFactor = 1; % gradients normalized to this value n = minStepDistanceSamples; @@ -69,9 +69,7 @@ warning on tapas_physio_findpeaks:largeMinPeakHeight % plus gains refer to max-changes in the future - - - idxGainPlus = idxGainPlus + n; +idxGainPlus = idxGainPlus + n; % + 1 because of diff idxGainMinus = idxGainMinus + 1; @@ -107,7 +105,8 @@ set(gcf, 'Name', stringTitle); hs(1) = subplot(2,1,1); hp(1) = plot(G); hold all; - hp(2) = plot(mG); + hp(2) = plot(mG, 'k.-'); + hp(2) = plot(-mG, 'k.-'); hp(3) = plot(dmG); hp(4) = stem(idxGainPlus, mG(idxGainPlus)); hp(5) = stem(idxGainMinus, mG(idxGainMinus)); diff --git a/PhysIO/code/tapas_physio_review.m b/PhysIO/code/tapas_physio_review.m index 8c4f5b99..75858b86 100644 --- a/PhysIO/code/tapas_physio_review.m +++ b/PhysIO/code/tapas_physio_review.m @@ -2,7 +2,7 @@ % Reviews performance of PhysIO (recreating output plots and text) % after running tapas_physio_main_create_regressors % -% fhArray = tapas_physio_review(input) +% verbose = tapas_physio_review(physio, newVerboseLevel) % % NOTE: Change physio.verbose.level before running this function to get % additonal output plots not seen during executing of the main-function diff --git a/PhysIO/code/tapas_physio_uddpvparse.m b/PhysIO/code/tapas_physio_uddpvparse.m deleted file mode 100644 index 6b6ff29d..00000000 --- a/PhysIO/code/tapas_physio_uddpvparse.m +++ /dev/null @@ -1,135 +0,0 @@ -function [h, errid, errmsg, msgObj] = tapas_physio_uddpvparse(c, varargin) -%tapas_physio_uddpvparse Parse p/v pairs using a udd object. -% -% NOTE: This copy of function uddpvparse is included in the TAPAS PhysIO -% Toolbox to make the dependency on Matlab's signal processing toolbox. -% explicit. Please do not use this function if you haven't purchased -% the signal processing toolbox. -% -% tapas_physio_uddpvparse(C) Returns a UDD object of class c. -% -% tapas_physio_uddpvparse(C, H) Returns H if it is of class C. Errors if it is not -% the correct class. -% -% tapas_physio_uddpvparse(C, P1, V1, P2, V2, etc.) Returns a UDD object of class C -% with its parameters P1, P2, etc. set to V1, V2, etc. -% -% tapas_physio_uddpvparse(C, CELLC, ...) If the second argument is a cell array, it -% will be used to construct the default object. -% -% [H, ERRID, ERRMSG] = tapas_physio_uddpvparse(...) Returns an error identifier in -% ERRID and an error message in ERRMSG. These will be empty when there -% are no errors. The returned ERRID and ERRMSG can later be passed -% directly to ERROR. This is used to avoid deep stack traces. If one H -% is requested, tapas_physio_uddpvparse will error when ERRMSG is not empty. -% -% [H, ERRID, ERRMSG, MSGOBJ] = tapas_physio_uddpvparse(...) Returns a message object -% that can be used to call error(MSGOBJ). MSGOBJ is empty when no error -% occurred. -% -% This function is meant to be called by functions which can take either -% an options object (C) or P/V pairs for the options object. -% -% Examples: -% -% % Constructs a default object -% tapas_physio_uddpvparse('fdopts.sosscaling') -% -% % Uses the input as the object -% h = fdopts.sosscaling; -% h.MaxNumerator = 10; -% tapas_physio_uddpvparse('fdopts.sosscaling', h) -% -% % Uses the P/V pairs -% tapas_physio_uddpvparse('fdopts.sosscaling', 'MaxNumerator', 10) -% -% % Uses the optional constructor syntax -% Hd = dfilt.df2sos; -% set(Hd, 'Arithmetic', 'Fixed', ... -% 'CoeffAutoScale', false, ... -% 'NumFracLength', 17); -% tapas_physio_uddpvparse('fdopts.sosscaling', {'scaleopts', Hd}, ... -% 'ScaleValueConstraint', 'po2') - -% Author(s): J. Schickler -% Copyright 1988-2008 The MathWorks, Inc. -% $Revision: 1.1.6.6 $ $Date: 2011/05/13 18:14:45 $ - -if verLessThan('matlab','8.4') - error(nargchk(1,11,nargin,'struct')); -else - narginchk(1,11); -end - -identifier = ''; - -h = []; - -if nargin == 1, - h = feval(c); -else - - % If the first extra is a cell array, it must be the "constructor". - if iscell(varargin{1}) - cons = varargin{1}; - varargin(1) = []; - else - cons = {c}; - end - - if length(varargin) == 0 - try - h = feval(cons{:}); - catch ME - throw(ME); - end - else - % Check if the first input is actually the object. Cannot use ISA - % because that returns true for subclasses. - if strcmpi(class(varargin{1}), c) - - % If the object is the input, just return it and ignore the rest of - % the inputs. - h = varargin{1}; - varargin(1) = []; - end - - if ~isempty(varargin) - - if rem(length(varargin), 2) - - % Cannot have an odd # of parameters - identifier = 'signal:tapas_physio_uddpvparse:invalidPVPairs'; - else - - % Assume that the rest of the inputs must be P/V pairs. If they - % are not they will be captured in the 'try/catch'. - if isempty(h) - h = feval(cons{:}); - end - try - set(h, varargin{:}); - catch ME - throw(ME); - end - end - end - end -end - -if ~isempty(identifier) - errid = identifier; - msgObj = message(identifier); - errmsg = getString(msgObj); -else - errid = ''; - errmsg = ''; - msgObj = []; -end - -% Only error if the msg structure is not requested. -if nargout == 1 - if ~isempty(errmsg), error(msgObj); end -end - -% [EOF] diff --git a/PhysIO/manual/QuickStart_PhysIO_Toolbox.pdf b/PhysIO/docs/QuickStart_PhysIO_Toolbox.pdf similarity index 100% rename from PhysIO/manual/QuickStart_PhysIO_Toolbox.pdf rename to PhysIO/docs/QuickStart_PhysIO_Toolbox.pdf diff --git a/PhysIO/docs/documentation.html b/PhysIO/docs/documentation.html new file mode 100644 index 00000000..64fcbb01 --- /dev/null +++ b/PhysIO/docs/documentation.html @@ -0,0 +1,1923 @@ + + + + + + + PhysIO Toolbox Documentation + + + + + + + + + + + + + +
+
+ + + +
+

Home

+

Wiki Main Page

+

Several documentation files are provided with this toolbox, and are linked in this wiki. They have the extension .md (markdown), i.e. are plain text files, but can be conveniently viewed online as the github/gitlab Wiki. +You can find them in physio/docs.

+

List of Documentation files:

+
    +
  • [x] HOME: This page. Landing Page of Wiki. Navigation to all other files and this explanation
  • +
  • [x] FAQ: Frequently asked questions (for users)
  • +
  • [x] QUICKSTART: Example script and how to use on test data, Intro to Batch Editor GUI
  • +
  • [x] EXAMPLES: List and explanation of all examples
  • +
  • [ ] MANUAL: Reference Manual (mostly for developers) listing all functions, and rationales of the toolbox, disecting its modular structure
  • +
  • [x] README: from main physio-public repository, similar info to this page
  • +
  • [x] CHANGELOG: List of all toolbox versions and the respective release notes, i.e. major changes in functionality, bugfixes etc.
  • +
+

Quickstart

+

Quickstart Manual

+

Purpose

+

This page provides simple walk-throughs of the SPM Batch Editor GUI, the scripts to run the main examples, and the most common output plots of the PhysIO Toolbox.

+

One-Page-Quickstart (with SPM)

+

...of SPM Batch Editor GUI for PhysIO Toolbox

+
    +
  1. Copy the PhysIO Toolbox physio/code folder to spm/toolbox
      +
    • (Optional) Rename the folder to something meaningful, e.g. PhysIO (see Figure 1).
    • +
    +
  2. +
  3. (Re-)Start SPM (spm fmri) and the Batch editor.
  4. +
  5. The PhysIO Toolbox should now occur under SPM -> Tools -> TAPAS PhysIO Toolbox
  6. +
  7. Change directory (!) to examples/Philips/ECG3T-folder and load an example spm_job file into the batch editor, e.g : example_spm_job_ECG3T.m
  8. +
  9. Press Play!
  10. +
+

+

Readme

+

TAPAS PhysIO Toolbox Version 2017

+
+

Copyright (C) 2012-2017 Lars Kasper kasper@biomed.ee.ethz.ch

+

Translational Neuromodeling Unit (TNU)

+

Institute for Biomedical Engineering

+

University of Zurich and ETH Zurich

+
+

Copying

+

The PhysIO Toolbox is free software: you can redistribute it and/or +modify it under the terms of the GNU General Public License as +published by the Free Software Foundation, either version 3 of the +License, or (at your option) any later version.

+

This program is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +General Public License for more details.

+

You should have received a copy of the GNU General Public License +along with this program (see the file COPYING). If not, see +http://www.gnu.org/licenses/.

+

Purpose

+

The PhysIO Toolbox provides model-based physiological noise correction of +fMRI data using peripheral measures of respiration and cardiac pulsation. +It incorporates noise models of cardiac/respiratory phase (RETROICOR, +Glover et al. 2000), as well as heart rate variability and respiratory +volume per time (cardiac response function, Chang et. al, 2009, respiratory +response function, Birn et al. 2006). The toolbox is usable via the SPM +batch editor, performs automatic pre-processing of noisy peripheral data +and outputs nuisance regressor files directly suitable for SPM +(multiple_regressors.txt).

+

Installation

+

Matlab

+
    +
  • Unzip the TAPAS archive
  • +
  • Add tapas/physio/code to your matlab path
  • +
+

SPM

+
    +
  • Certain functionality (Batch Editor GUI, pipeline dependencies, model assessment via F-contrasts) require the installation of SPM
  • +
  • Afterwards, the PhysIO Toolbox has to be registered as an SPM toolbox by copying the physio/code folder to spm/toolbox/physio
  • +
+

Getting Started

+

Run example_main_ECG3T.m in subdirectory Philips/ECG3T of the toolbox example repository physio-examples. +See subdirectory physio/docs and next section of this document.

+

Getting Help/Documentation

+

Several documentation files are provided with this toolbox. They have the extension .md (markdown), i.e. are plain text files, but can be conveniently viewed online as the github/gitlab Wiki. +You can find them in physio/wikidocs.

+

Alternatively, a pdf converted version of the following files is found in physio/docs/documentation.pdf

+

List of Documentation files:

+
    +
  • README.md: this file, purpose, installation, getting started, pointer to more help
  • +
  • FAQ.md: Frequently asked questions (for users)
  • +
  • QUICKSTART.md: Example script and how to use on test data, Intro to Batch Editor GUI
  • +
  • MANUAL.md: Reference Manual (mostly for developers) listing all functions, and rationales of the toolbox, disecting its modular structure
  • +
  • WIKIMAIN.md: Landing Page of Wiki. Navigation to all other files and this explanation
  • +
  • EXAMPLES.md: List and explanation of all examples
  • +
  • CHANGELOG.md: List of all toolbox versions and the respective release notes, i.e. major changes in functionality, bugfixes etc.
  • +
+

Background

+

The PhysIO Toolbox provides physiological noise correction for fMRI-data +from peripheral measures (ECG/pulse oximetry, breathing belt). It is +model-based, i.e. creates nuisance regressors from the physiological +monitoring that can enter a General Linear Model (GLM) analysis, e.g. +SPM8/12. Furthermore, for scanner vendor logfiles (PHILIPS, GE, Siemens), +it provides means to statistically assess peripheral data (e.g. heart rate variability) +and recover imperfect measures (e.g. distorted R-peaks of the ECG).

+

Facts about physiological noise in fMRI:

+
    +
  • Physiological noise can explain 20-60 % of variance in fMRI voxel time +series (Birn2006, Hutton2011, Harvey2008)
      +
    • Physiological noise affects a lot of brain regions (s. figure, e.g. +brainstem or OFC), especially next to CSF, arteries (Hutton2011).
    • +
    • If not accounted for, this is a key factor limiting sensitivity for effects of interest.
    • +
    +
  • +
  • Physiological noise contributions increase with field strength; they +become a particular concern at and above 3 Tesla (Kasper2009, Hutton2011).
  • +
  • In resting state fMRI, disregarding physiological noise leads to wrong +connectivity results (Birn2006).
  • +
+

Therefore, some kind of physiological noise correction is highly recommended for every statistical fMRI analysis.

+

Model-based correction of physiological noise:

+
    +
  • Physiological noise can be decomposed into periodic time series following +heart rate and breathing cycle.
  • +
  • The Fourier expansion of cardiac and respiratory phases was introduced as +RETROICOR (RETROspective Image CORrection, Glover2000, +see also Josephs1997).
  • +
  • These Fourier Terms can enter a General Linear Model (GLM) as nuisance +regressors, analogous to movement parameters.
  • +
  • As the physiological noise regressors augment the GLM and explain +variance in the time series, they increase sensitivity in all contrasts +of interest.
  • +
+

Features of this Toolbox

+

Physiological Noise Modeling

+
    +
  • Modeling physiological noise regressors from peripheral data +(breathing belt, ECG, pulse oximeter)
      +
    • State of the art RETROICOR cardiac and respiratory phase expansion
    • +
    • Cardiac response function (Chang et al, 2009) and respiratory response +function (Birn et al. 2006) modelling of heart-rate variability and +respiratory volume per time influence on physiological noise
    • +
    • Flexible expansion orders to model different contributions of cardiac, +respiratory and interaction terms (see Harvey2008, Hutton2011)
    • +
    +
  • +
  • Data-driven noise regressors
      +
    • PCA extraction from nuisance ROIs (CSF, white matter), similar to aCompCor (Behzadi2007)
    • +
    +
  • +
+

Automatization and Performance Assessment

+
    +
  • Automatic creation of nuisance regressors, full integration into standard +GLMs, tested for SPM8/12 ("multiple_regressors.mat")
  • +
  • Integration in SPM Batch Editor: GUI for parameter input, dependencies to integrate physiological noise correction in preprocessing pipeline
  • +
  • Performance Assessment: Automatic F-contrast and tSNR Map creation and display for groups of physiological noise regressors, using SPM GLM tools
  • +
+

Flexible Read-in

+

The toolbox is dedicated to seamless integration into a clinical research s +etting and therefore offers correction methods to recover physiological +data from imperfect peripheral measures.

+
    +
  • General Electric
  • +
  • Philips SCANPHYSLOG files (all versions from release 2.6 to 5.3)
  • +
  • Siemens VB (files .ecg, .resp, .puls
  • +
  • Siemens VD (files (*_ECG.log, *_RESP.log, *_PULS.log)
  • +
  • Biopac .mat-export
      +
    • assuming the following variables (as columns): data, isi, isi_units, labels, start_sample, units
    • +
    • See tapas_physio_read_physlogfiles_biopac_mat.m for details
    • +
    +
  • +
  • Custom logfiles: should contain one amplitude value per line, one logfile per device. Sampling interval(s) are provided as a separate parameter to the toolbox.
  • +
+

Compatibility and Support

+
    +
  • Matlab Toolbox
  • +
  • Input:
      +
    • Fully integrated to work with physiological logfiles for Philips MR systems (SCANPHYSLOG)
    • +
    • tested for General Electric (GE) log-files
    • +
    • implementation for Siemens log-files
    • +
    • also: interface for 'Custom', i.e. general heart-beat time stamps +& breathing volume time courses from other log formats
    • +
    +
  • +
  • Output:
      +
    • Nuisance regressors for mass-univariate statistical analysis with SPM5,8,12 +or as text file for export to any other package
    • +
    • raw and processed physiological logfile data
    • +
    +
  • +
  • Part of the TAPAS Software Collection of the Translational Neuromodeling Unit (TNU) Zurich:long term support and ongoing development
  • +
+

Contributors

+
    +
  • Lead Programmer:
      +
    • Lars Kasper, TNU & MR-Technology Group, IBT, University of Zurich & ETH Zurich
    • +
    +
  • +
  • Project Team:
      +
    • Steffen Bollmann, Centre for Advanced Imaging, University of Queensland, Australia
    • +
    • Saskia Bollmann, Centre for Advanced Imaging, University of Queensland, Australia
    • +
    +
  • +
  • Contributors:
      +
    • Eduardo Aponte, TNU Zurich
    • +
    • Tobias U. Hauser, FIL London, UK
    • +
    • Jakob Heinzle, TNU Zurich
    • +
    • Chloe Hutton, FIL London, UK (previously)
    • +
    • Miriam Sebold, Charite Berlin, Germany
    • +
    +
  • +
+

Contact

+

Send bug reports and suggestions either to +1) our mailing list: tapas@sympa.ethz.ch, or +2) as an issue on our TAPAS github account : https://github.com/translationalneuromodeling/tapas/issues

+

References

+

Main Toolbox Reference

+
    +
  1. Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. Journal of Neuroscience Methods 276, 56–72. doi:10.1016/j.jneumeth.2016.10.019
  2. +
+ +
    +
  1. Glover, G.H., Li, T.Q. & Ress, D. Image‐based method for retrospective correction +of PhysIOlogical motion effects in fMRI: RETROICOR. Magn Reson Med 44, 162-7 (2000).

    +
  2. +
  3. Hutton, C. et al. The impact of PhysIOlogical noise correction on fMRI at 7 T. +NeuroImage 57, 101‐112 (2011).

    +
  4. +
  5. Harvey, A.K. et al. Brainstem functional magnetic resonance imaging: +Disentangling signal from PhysIOlogical noise. Journal of Magnetic Resonance +Imaging 28, 1337‐1344 (2008).

    +
  6. +
  7. Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise +correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, +90–101. doi:10.1016/j.neuroimage.2007.04.042

    +
  8. +
  9. Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. The respiration response +function: The temporal dynamics of fMRI s ignal fluctuations related to changes in +respiration. NeuroImage 40, 644–654. doi:10.1016/j.neuroimage.2007.11.059

    +
  10. +
  11. Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate on the +BOLD signal: The cardiac response function. NeuroImage 44, 857–869. +doi:10.1016/j.neuroimage.2008.09.029

    +
  12. +
  13. Siegel, J.S., Power, J.D., Dubis, J.W., Vogel, A.C., Church, J.A., Schlaggar, B.L., +Petersen, S.E., 2014. Statistical improvements in functional magnetic resonance +imaging analyses produced by censoring high-motion data points. Hum. Brain Mapp. +35, 1981–1996. doi:10.1002/hbm.22307

    +
  14. +
+

FAQ

+

Frequently Asked Questions (FAQ)

+

1. What is the PhysIO Toolbox?

+

PhysIO is a toolbox for model-based physiological noise correction of fMRI data.

+

PhysIO stands for Physiological Input/Output toolbox, which summarizes its core purpose. A quote from our paper:

+
+

In short, the toolbox transforms physiological input, i.e. peripheral recordings, into physiological output, i.e. regressors encoding components of physiological noise [...] A modular Matlab implementation supports command-line operation and is compatible with all major fMRI analysis packages via the export of regressor text-files. For the Statistical Parametric Mapping SPM software package in particular, PhysIO features a full integration as a Batch Editor Tool, which allows user-friendly, GUI-based setup and inclusion into existing preprocessing and modeling pipelines.

+
+

2. How does PhysIO differ from other toolboxes for physiological noise correction for fMRI using peripheral recordings?

+

Citing from the introduction of our paper again

+

>

+
+

Highlights

+
    +
  • A Toolbox to integrate preprocessing of physiological data and fMRI noise modeling.
  • +
  • Robust preprocessing via iterative peak detection, shown for noisy data and patients.
  • +
  • Flexible support of peripheral data formats and noise models (RETROICOR, RVHRCOR).
  • +
  • Fully automated noise correction and performance assessment for group studies.
  • +
  • Integration in fMRI pre-processing pipelines as SPM Toolbox (Batch Editor GUI).
  • +
+
+

3. How do I cite PhysIO?

+

The core reference for PhysIO is: The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data (http://dx.doi.org/10.1016/j.jneumeth.2016.10.019)

+

Please cite this paper if you use PhysIO in your work. Moreover, this paper is also a good source for more information on PhysIO (see next question).

+

A standard snippet to include in your method section could look like the following, assuming you use our specific implementation of RETROICOR, which uses Fourier expansions of different order for the estimated phases of cardiac pulsation (3rd order), respiration (4th order) and cardio-­‐respiratory interactions (1st order) following (Harvey et al., 2008)

+
+

Correction for physiological noise was performed via RETROICOR [1,2] using Fourier +expansions of different order for the estimated phases of cardiac pulsation (3rd order), +respiration (4th order) and cardio-­‐respiratory interactions (1st order) [2]: The +corresponding confound regressors were created using the Matlab PhysIO Toolbox ([4], +open source code available as part of the TAPAS software collection: +https://www.translationalneuromodeling.org/tapas).

+
+
    +
  1. Glover, G.H., Li, T.Q. & Ress, D. Image-­‐based method for retrospective correction +of PhysIOlogical motion effects in fMRI: RETROICOR. Magn Reson Med 44, 162-­‐ +7 (2000).

    +
  2. +
  3. Hutton, C. et al. The impact of PhysIOlogical noise correction on fMRI at 7 T. +NeuroImage 57, 101-­‐112 (2011).

    +
  4. +
  5. Harvey, A.K. et al. Brainstem functional magnetic resonance imaging: +Disentangling signal from PhysIOlogical noise. Journal of Magnetic Resonance +Imaging 28, 1337-­‐1344 (2008).

    +
  6. +
  7. Kasper, L., Bollmann, S., Diaconescu, A.O., Hutton, C., Heinzle, J., Iglesias, S., Hauser, T.U., Sebold, M., Manjaly, Z.-M., Pruessmann, K.P., Stephan, K.E., 2017. The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data. Journal of Neuroscience Methods 276, 56–72. doi:10.1016/j.jneumeth.2016.10.019

    +
  8. +
+

If you use respiratory‐volume-­per time (RVT), heart-‐rate +variability (HRV), noise ROIs or 12/24 regressor motion modeling, also include the +respective references:

+
    +
  1. Behzadi, Y., Restom, K., Liau, J., Liu, T.T., 2007. A component based noise +correction method (CompCor) for BOLD and perfusion based fMRI. NeuroImage 37, +90–101. doi:10.1016/j.neuroimage.2007.04.042

    +
  2. +
  3. Birn, R.M., Smith, M.A., Jones, T.B., Bandettini, P.A., 2008. The respiration response +function: The temporal dynamics of fMRI s ignal fluctuations related to changes in +respiration. NeuroImage 40, 644–654. doi:10.1016/j.neuroimage.2007.11.059 +PhysIO Toolbox | Citing this work 20

    +
  4. +
  5. Chang, C., Cunningham, J.P., Glover, G.H., 2009. Influence of heart rate on the +BOLD signal: The cardiac response function. NeuroImage 44, 857–869. +doi:10.1016/j.neuroimage.2008.09.029

    +
  6. +
  7. Siegel, J.S., Power, J.D., Dubis, J.W., Vogel, A.C., Church, J.A., Schlaggar, B.L., +Petersen, S.E., 2014. Statistical improvements in functional magnetic resonance +imaging analyses produced by censoring high-motion data points. Hum. Brain Mapp. +35, 1981–1996. doi:10.1002/hbm.22307

    +
  8. +
+

4. Where do I find more documentation for PhysIO?

+
    +
  • The paper describing its structure, objective and modules
  • +
  • README.md in the main folder when downloading
      +
    • For help on installation and getting started
    • +
    +
  • +
  • Quickstart
      +
    • PDF (or markdown .md file)
    • +
    • Tutorial matlab-scripts
    • +
    +
  • +
  • Reference Manual (for developers)
  • +
+

5. I am using FSL, AFNI, BrainVoyager, etc., for my fMRI analyses. Do I need SPM for PhysIO to work?

+

No, the basic functionality of PhysIO, i.e. creating nuisance regressors for your GLM analysis, is available in plain Matlab. The following extra functionality related to automatizing and assessing noise correction, require the installation of SPM:

+
    +
  • GUI (SPM Batch Editor)
  • +
  • Pipeline dependencies (automatic input of realignment parameters, feed-in of multiple regressors file to GLM)
  • +
  • Model assessment via F-tests and automatic F-map/tSNR report
  • +
  • Noise-ROIs model (read-in of nifti files via SPM)
  • +
+

6. I am using device X for physiological recordings. Does PhysIO support the physiological logfile format Y?

+

Currently, PhysIO natively supports the following physiological logfile types:

+
    +
  • General Electric
  • +
  • Philips SCANPHYSLOG files (all versions from release 2.6 to 5.3)
  • +
  • Siemens VB (files .ecg, .resp, .puls
  • +
  • Siemens VD (files (*_ECG.log, *_RESP.log, *_PULS.log)
  • +
  • Biopac .mat-export
      +
    • assuming the following variables (as columns): data, isi, isi_units, labels, start_sample, units
    • +
    • See tapas_physio_read_physlogfiles_biopac_mat.m for details
    • +
    +
  • +
+

Furthermore, physiological recordings can be entered via a custom data format, i.e., providing one text file per device. The files should contain one amplitude value per line. The corresponding sampling interval(s) are provided as a separate parameter in the toolbox.

+

If your favourite logfile format is not supported, please contact the developers. We try everything to accomodate the read-in flexibility of the toolbox to your needs.

+

7. I am running the toolbox for a lot of subjects / on a remote server without graphics. Can I somehow reproduce the output figures relevant to assess the data quality?

+

Yes you can, using the toolbox function tapas_physio_review. This function takes the physio-structure as an input argument, which is per default saved as physio.mat in the specified output folder of your batch job.

+

8. How do I interpret the various output plots of the toolbox?

+

Have a look at our publication: The PhysIO Toolbox for Modeling Physiological Noise in fMRI Data (http://dx.doi.org/10.1016/j.jneumeth.2016.10.019)

+

The figures there give a good overview of the toolbox output figures.

+

9. I want to access subject's physiological measures, e.g. heart rate or respiratory volume (per time), before they enter the regressors. Where can I do that?

+

All intermediate data processing steps (e.g. filtering, cropping) of the peripheral data, including the computation of physiologically meaningful time courses, such as heart rate and respiratory volume, are saved in the substructure ons_secs ("onsets in seconds) of the physio-structure mentioned in question 7. This structure is typically saved in a file physio.mat.

+

physio.ons_secs then contains the different time courses, cropped to the acquisition window synchronized to your fMRI scan (the same values before synchronization/cropping, is found in physio.ons_secs.raw). Here are the most important ones:

+
    +
  • ons_secs.t = []; % time vector corresponding to c and r
  • +
  • ons_secs.c = []; % raw cardiac waveform (ECG or PPU)
  • +
  • ons_secs.r = []; % raw respiration amplitude time course
  • +
  • ons_secs.cpulse = []; % onset times of cardiac pulse events (e.g. R-peaks)
  • +
  • ons_secs.fr = []; % filtered respiration amplitude time series
  • +
  • ons_secs.c_sample_phase = []; % phase in heart-cycle when each slice of each volume was acquired
  • +
  • ons_secs.r_sample_phase = []; % phase in respiratory cycle when each slice of each volume was acquired
  • +
  • ons_secs.hr = []; % [nScans,1] estimated heart rate at each scan
  • +
  • ons_secs.rvt = []; % [nScans,1] estimated respiratory volume per time at each scan
  • +
  • ons_secs.c_outliers_high = []; % onset of too long heart beats
  • +
  • ons_secs.c_outliers_low = []; % onsets of too short heart beats
  • +
  • ons_secs.r_hist = []; % histogram of breathing amplitudes
  • +
+

For a detailed list of all properties and their documentation, read the source code of tapas_physio_new.m

+

10. What is the order of the regressor columns in the multiple regressors file?

+

This depends on the physiological models (and their order) specified in the model-submodule of physio (or in the batch editor). The general order is outlined in Fig. 7A of the http://dx.doi.org/10.1016/j.jneumeth.2016.10.019. The []-brackets indicate the number of regressors:

+
    +
  1. RETROICOR cardiac regressors [2 x nOrderCardiac]
  2. +
  3. RETROICOR respiratory regressors [2 x nOrderRespiratory]
  4. +
  5. RETROICOR cardXResp interaction regressors [4 x nOrderCardiacXRespiratory]
  6. +
  7. HRV [nDelaysHRV]
  8. +
  9. RVT [nDelaysRVT]
  10. +
  11. Noise ROIs (PCA signatures and mean of each region) [nNoiseROIs x (nComponents+1)]
  12. +
  13. Other (included other text file) [nColumnsOtherFile]
  14. +
  15. Motion [6 or 12 or 24, depending on motion model]
  16. +
+

If any of the models was not specified, the number of regressors is reduced accordingly.

+

11. I cannot find the answer to my question in the FAQ. Whom do I ask for help?

+

Subscribe to our TAPAS mailing list by clicking Subscribe on the left side of this website

+

Afterwards you can send e-mails with your questions to tapas@sympa.ethz.ch. Both, core developers of PhysIO and experienced users are receiving your e-mail and are eager to help.

+

The mailing list also has a searchable archive (click Archive on the left of above-mentioned website), where you might already find the answer to your question.

+

Examples

+

Example Datasets for PhysIO

+

The following datasets are available and can be downloaded with the toolbox in a specific physio-examples repository, which can be cloned from git@tnurepository.ethz.ch:physio/physio-examples.git.

+

Besides the raw physiological logfiles, each example contains example scripts to run physio as

+
    +
  • SPM job (*spm_job.mat)
  • +
  • editable SPM job (*spm_job.m)
  • +
  • plain matlab script (*matlab_script.m)
  • +
+

Philips

+

ECG 3T

+

Courtesy of Sandra Iglesias, Translational Neuromodeling Unit, ETH & +University of Zurich

+

4-electrode ECG and breathing belt, Philips 3T Achieva scanner

+

Description: Standard example; shows how to use volume counting either +from beginning or end of run to synchronize physiological logfile with +acquisition onsets of fMRI scans.

+

ECG 7T

+

Courtesy of Zina-Mary Manjaly, University Hospital Zurich

+

4-electrode ECG and breathing belt, Philips 7T Achieva scanner

+

Description: The ECG data for ultra-high field data is typically much +noisier than at 3 Tesla. Therefore, R-wave peaks are frequently missed +by prospective trigger detection and not marked correctly in the +logfile. This example shows how to select typical R-wave-peaks manually +and let the algorithm find the heartbeat events.

+

PPU 3T

+

Courtesy of Diana Wotruba, University and University Hospital of Zurich

+

PPU (finger plethysmograph) and breathing belt, Philips 3T Achieva +scanner

+

Description: Similar to ECG3T, but a plethysmograph instead of an ECG +was used to monitor the cardiac pulsation. Example shows how to extract +heart and breathing rate.

+

General Electric

+

PPU 3T

+

Courtesy of Steffen Bollmann, Kinderspital Zurich and ETH Zurich

+

PPU (finger plethysmograph) and breathing belt, General Electric 3T +scanner

+

Description: Similar to PPU, but acquired on a GE system with two +separate output logfiles for pulse oximetry and breathing amplitude, +sampled with 40 Hz. The quality of the signal is particularly +challenging, stemming from a patient population.

+

Siemens

+

ECG 3T

+

Courtesy of Miriam Sebold, Charite Berlin, and Quentin Huys, TNU Zurich

+

4-electrode ECG data, Siemens 3T scanner

+

Description: Similar to ECG 3T, but acquired on a Siemens system with only one logfile for ECG data. The quality of the signal is challenging, stemming from a patient population.

+

Changelog/Versions

+

RELEASE INFORMATION

+

Current Release

+

PhysIO_Toolbox_R2017.1

+

February 19, 2017

+

Minor Release Notes (R2017.1)

+
    +
  • Substantially improved Siemens interface, both for VB/VD and 3T/7T releases
      +
    • several bugfixes
    • +
    • based on extensive user feedback from Berlin and Brisbane
    • +
    +
  • +
  • New functionality tapas_physio_overlay_contrasts.m to display non-physio +contrasts automatically as well
  • +
+

Major Release Notes (r904 / R2016.1)

+
    +
  • Software version for accepted PhysIO Toolbox Paper: doi:10.1016/j.jneumeth.2016.10.019
  • +
  • Tested and expanded versions of examples
  • +
  • Improved stability by bugfixes and compatibility to Matlab R2016
  • +
  • Slice-wise regressor creation
  • +
  • Detection of constant physiological time series (detachment, clipping)
  • +
  • Refactoring of report_contrasts and compute_tsnr_gains as standalone functionality
  • +
  • Improved Read-in capabilities (Siemens respiration data, BioPac .mat)
  • +
  • Migration from svn (r904) to git (tnurepository) for version control
  • +
+

Major Release Notes (r835)

+
    +
  • Software version for Toolbox Paper submission
  • +
  • Noise ROIs modeling
  • +
  • Extended motion models (24 parameters, Volterra expansion)
  • +
  • HRV/RVT models with optional multiple delay regressors
  • +
  • Report_contrasts with automatic contrast generation for all regressor groups
  • +
  • compute_tsnr_gains for individual physiological regressor groups
  • +
  • consistent module naming (scan_timing, preproc)
  • +
  • Visualisation improvement (color schemes, legends)
  • +
+

Minor Release Notes (r666)

+
    +
  • Compatibility tested for SPM12, small bugfixes Batch Dependencies
  • +
  • Cleaner Batch Interface with grouped sub-menus (cfg_choice)
  • +
  • new model: 'none' to just read out physiological raw data and preprocess, +without noise modelling
  • +
  • Philips: Scan-timing via gradient log now automatized (gradient_log_auto)
  • +
  • Siemens: Tics-Logfile read-in (proprietary, needs Siemens-agreement)
  • +
  • All peak detections (cardiac/respiratory) now via auto_matched algorithm
  • +
  • Adapt plots/saving for Matlab R2014b
  • +
+

Major Release Notes (r534)

+
    +
  • Read-in of Siemens plain text log files; new example dataset for Siemens
  • +
  • Speed up and debugging of auto-detection method for noisy cardiac data => new method thresh.cardiac.initial_cpulse_select.method = ???auto_matched???
  • +
  • Error handling for temporary breathing belt failures (Eduardo Aponte, TNU Zurich)
  • +
  • slice-wise regressors can be created by setting sqpar.onset_slice to a index vector of slices
  • +
+

Major Release Notes (r497)

+
    +
  • SPM matlabbatch GUI implemented (Call via Batch -> SPM -> Tools -> TAPAS PhysIO Toolbox)
  • +
  • improved, automatic heartbeat detection for noisy ECG now standard for ECG and Pulse oximetry (courtesy of Steffen Bollmann)
  • +
  • QuickStart-Manual and PhysIO-Background presentation expanded/updated
  • +
  • job .m/.mat-files created for all example datasets
  • +
  • bugfixes cpulse-initial-select method-handling (auto/manual/load)
  • +
+

Major Release Notes (r429)

+
    +
  • Cardiac and Respiratory response function regressors integrated in workflow (heart rate and breathing volume computation)
  • +
  • Handling of Cardiac and Respiratory Logfiles only
  • +
  • expanded documentation (Quickstart.pdf and Handbook.pdf)
  • +
  • read-in of custom log files, e.g. for BrainVoyager peripheral data
  • +
  • more informative plots and commenting (especially in tapas_physio_new).
  • +
+

Minor Release Notes (r354)

+
    +
  • computation of heart and breathing rate in Philips/PPU/main_PPU.m
  • +
  • prefix of functions with tapas_*
  • +
+

Major Release Notes (r241)

+
    +
  • complete modularization of reading/preprocessing/regressor creation for peripheral physiological data
  • +
  • manual selection of missed heartbeats in ECG/pulse oximetry (courtesy of Jakob Heinzle)
  • +
  • support for logfiles from GE scanners (courtesy of Steffen Bollmann, KiSpi Zuerich)
  • +
  • improved detection of pulse oximetry peaks (courtesy of Steffen Bollmann)
  • +
  • improved documentation
  • +
  • consistent function names (prefixed by "physio_")
  • +
+

NOTE: Your main_ECG/PPU.m etc. scripts from previous versions (<=r159) will not work with this one any more. Please adapt one of the example scripts for your needs (~5 min of work). The main benefit of this version is a complete new variable structure that is more sustainable and makes the code more readable.

+

License

+
                GNU GENERAL PUBLIC LICENSE
+                   Version 3, 29 June 2007

Copyright (C) 2007 Free Software Foundation, Inc. http://fsf.org/ + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed.

+
                        Preamble

The GNU General Public License is a free, copyleft license for +software and other kinds of works.

+

The licenses for most software and other practical works are designed +to take away your freedom to share and change the works. By contrast, +the GNU General Public License is intended to guarantee your freedom to +share and change all versions of a program--to make sure it remains free +software for all its users. We, the Free Software Foundation, use the +GNU General Public License for most of our software; it applies also to +any other work released this way by its authors. You can apply it to +your programs, too.

+

When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +them if you wish), that you receive source code or can get it if you +want it, that you can change the software or use pieces of it in new +free programs, and that you know you can do these things.

+

To protect your rights, we need to prevent others from denying you +these rights or asking you to surrender the rights. Therefore, you have +certain responsibilities if you distribute copies of the software, or if +you modify it: responsibilities to respect the freedom of others.

+

For example, if you distribute copies of such a program, whether +gratis or for a fee, you must pass on to the recipients the same +freedoms that you received. You must make sure that they, too, receive +or can get the source code. And you must show them these terms so they +know their rights.

+

Developers that use the GNU GPL protect your rights with two steps: +(1) assert copyright on the software, and (2) offer you this License +giving you legal permission to copy, distribute and/or modify it.

+

For the developers' and authors' protection, the GPL clearly explains +that there is no warranty for this free software. For both users' and +authors' sake, the GPL requires that modified versions be marked as +changed, so that their problems will not be attributed erroneously to +authors of previous versions.

+

Some devices are designed to deny users access to install or run +modified versions of the software inside them, although the manufacturer +can do so. This is fundamentally incompatible with the aim of +protecting users' freedom to change the software. The systematic +pattern of such abuse occurs in the area of products for individuals to +use, which is precisely where it is most unacceptable. Therefore, we +have designed this version of the GPL to prohibit the practice for those +products. If such problems arise substantially in other domains, we +stand ready to extend this provision to those domains in future versions +of the GPL, as needed to protect the freedom of users.

+

Finally, every program is threatened constantly by software patents. +States should not allow patents to restrict development and use of +software on general-purpose computers, but in those that do, we wish to +avoid the special danger that patents applied to a free program could +make it effectively proprietary. To prevent this, the GPL assures that +patents cannot be used to render the program non-free.

+

The precise terms and conditions for copying, distribution and +modification follow.

+
                   TERMS AND CONDITIONS
    +
  1. Definitions.

    +

    "This License" refers to version 3 of the GNU General Public License.

    +

    "Copyright" also means copyright-like laws that apply to other kinds of +works, such as semiconductor masks.

    +

    "The Program" refers to any copyrightable work licensed under this +License. Each licensee is addressed as "you". "Licensees" and +"recipients" may be individuals or organizations.

    +

    To "modify" a work means to copy from or adapt all or part of the work +in a fashion requiring copyright permission, other than the making of an +exact copy. The resulting work is called a "modified version" of the +earlier work or a work "based on" the earlier work.

    +

    A "covered work" means either the unmodified Program or a work based +on the Program.

    +

    To "propagate" a work means to do anything with it that, without +permission, would make you directly or secondarily liable for +infringement under applicable copyright law, except executing it on a +computer or modifying a private copy. Propagation includes copying, +distribution (with or without modification), making available to the +public, and in some countries other activities as well.

    +

    To "convey" a work means any kind of propagation that enables other +parties to make or receive copies. Mere interaction with a user through +a computer network, with no transfer of a copy, is not conveying.

    +

    An interactive user interface displays "Appropriate Legal Notices" +to the extent that it includes a convenient and prominently visible +feature that (1) displays an appropriate copyright notice, and (2) +tells the user that there is no warranty for the work (except to the +extent that warranties are provided), that licensees may convey the +work under this License, and how to view a copy of this License. If +the interface presents a list of user commands or options, such as a +menu, a prominent item in the list meets this criterion.

    +
  2. +
  3. Source Code.

    +

    The "source code" for a work means the preferred form of the work +for making modifications to it. "Object code" means any non-source +form of a work.

    +

    A "Standard Interface" means an interface that either is an official +standard defined by a recognized standards body, or, in the case of +interfaces specified for a particular programming language, one that +is widely used among developers working in that language.

    +

    The "System Libraries" of an executable work include anything, other +than the work as a whole, that (a) is included in the normal form of +packaging a Major Component, but which is not part of that Major +Component, and (b) serves only to enable use of the work with that +Major Component, or to implement a Standard Interface for which an +implementation is available to the public in source code form. A +"Major Component", in this context, means a major essential component +(kernel, window system, and so on) of the specific operating system +(if any) on which the executable work runs, or a compiler used to +produce the work, or an object code interpreter used to run it.

    +

    The "Corresponding Source" for a work in object code form means all +the source code needed to generate, install, and (for an executable +work) run the object code and to modify the work, including scripts to +control those activities. However, it does not include the work's +System Libraries, or general-purpose tools or generally available free +programs which are used unmodified in performing those activities but +which are not part of the work. For example, Corresponding Source +includes interface definition files associated with source files for +the work, and the source code for shared libraries and dynamically +linked subprograms that the work is specifically designed to require, +such as by intimate data communication or control flow between those +subprograms and other parts of the work.

    +

    The Corresponding Source need not include anything that users +can regenerate automatically from other parts of the Corresponding +Source.

    +

    The Corresponding Source for a work in source code form is that +same work.

    +
  4. +
  5. Basic Permissions.

    +

    All rights granted under this License are granted for the term of +copyright on the Program, and are irrevocable provided the stated +conditions are met. This License explicitly affirms your unlimited +permission to run the unmodified Program. The output from running a +covered work is covered by this License only if the output, given its +content, constitutes a covered work. This License acknowledges your +rights of fair use or other equivalent, as provided by copyright law.

    +

    You may make, run and propagate covered works that you do not +convey, without conditions so long as your license otherwise remains +in force. You may convey covered works to others for the sole purpose +of having them make modifications exclusively for you, or provide you +with facilities for running those works, provided that you comply with +the terms of this License in conveying all material for which you do +not control copyright. Those thus making or running the covered works +for you must do so exclusively on your behalf, under your direction +and control, on terms that prohibit them from making any copies of +your copyrighted material outside their relationship with you.

    +

    Conveying under any other circumstances is permitted solely under +the conditions stated below. Sublicensing is not allowed; section 10 +makes it unnecessary.

    +
  6. +
  7. Protecting Users' Legal Rights From Anti-Circumvention Law.

    +

    No covered work shall be deemed part of an effective technological +measure under any applicable law fulfilling obligations under article +11 of the WIPO copyright treaty adopted on 20 December 1996, or +similar laws prohibiting or restricting circumvention of such +measures.

    +

    When you convey a covered work, you waive any legal power to forbid +circumvention of technological measures to the extent such circumvention +is effected by exercising rights under this License with respect to +the covered work, and you disclaim any intention to limit operation or +modification of the work as a means of enforcing, against the work's +users, your or third parties' legal rights to forbid circumvention of +technological measures.

    +
  8. +
  9. Conveying Verbatim Copies.

    +

    You may convey verbatim copies of the Program's source code as you +receive it, in any medium, provided that you conspicuously and +appropriately publish on each copy an appropriate copyright notice; +keep intact all notices stating that this License and any +non-permissive terms added in accord with section 7 apply to the code; +keep intact all notices of the absence of any warranty; and give all +recipients a copy of this License along with the Program.

    +

    You may charge any price or no price for each copy that you convey, +and you may offer support or warranty protection for a fee.

    +
  10. +
  11. Conveying Modified Source Versions.

    +

    You may convey a work based on the Program, or the modifications to +produce it from the Program, in the form of source code under the +terms of section 4, provided that you also meet all of these conditions:

    +

    a) The work must carry prominent notices stating that you modified +it, and giving a relevant date.

    +

    b) The work must carry prominent notices stating that it is +released under this License and any conditions added under section

    +
      +
    1. This requirement modifies the requirement in section 4 to +"keep intact all notices".
    2. +
    +

    c) You must license the entire work, as a whole, under this +License to anyone who comes into possession of a copy. This +License will therefore apply, along with any applicable section 7 +additional terms, to the whole of the work, and all its parts, +regardless of how they are packaged. This License gives no +permission to license the work in any other way, but it does not +invalidate such permission if you have separately received it.

    +

    d) If the work has interactive user interfaces, each must display +Appropriate Legal Notices; however, if the Program has interactive +interfaces that do not display Appropriate Legal Notices, your +work need not make them do so.

    +

    A compilation of a covered work with other separate and independent +works, which are not by their nature extensions of the covered work, +and which are not combined with it such as to form a larger program, +in or on a volume of a storage or distribution medium, is called an +"aggregate" if the compilation and its resulting copyright are not +used to limit the access or legal rights of the compilation's users +beyond what the individual works permit. Inclusion of a covered work +in an aggregate does not cause this License to apply to the other +parts of the aggregate.

    +
  12. +
  13. Conveying Non-Source Forms.

    +

    You may convey a covered work in object code form under the terms +of sections 4 and 5, provided that you also convey the +machine-readable Corresponding Source under the terms of this License, +in one of these ways:

    +

    a) Convey the object code in, or embodied in, a physical product +(including a physical distribution medium), accompanied by the +Corresponding Source fixed on a durable physical medium +customarily used for software interchange.

    +

    b) Convey the object code in, or embodied in, a physical product +(including a physical distribution medium), accompanied by a +written offer, valid for at least three years and valid for as +long as you offer spare parts or customer support for that product +model, to give anyone who possesses the object code either (1) a +copy of the Corresponding Source for all the software in the +product that is covered by this License, on a durable physical +medium customarily used for software interchange, for a price no +more than your reasonable cost of physically performing this +conveying of source, or (2) access to copy the +Corresponding Source from a network server at no charge.

    +

    c) Convey individual copies of the object code with a copy of the +written offer to provide the Corresponding Source. This +alternative is allowed only occasionally and noncommercially, and +only if you received the object code with such an offer, in accord +with subsection 6b.

    +

    d) Convey the object code by offering access from a designated +place (gratis or for a charge), and offer equivalent access to the +Corresponding Source in the same way through the same place at no +further charge. You need not require recipients to copy the +Corresponding Source along with the object code. If the place to +copy the object code is a network server, the Corresponding Source +may be on a different server (operated by you or a third party) +that supports equivalent copying facilities, provided you maintain +clear directions next to the object code saying where to find the +Corresponding Source. Regardless of what server hosts the +Corresponding Source, you remain obligated to ensure that it is +available for as long as needed to satisfy these requirements.

    +

    e) Convey the object code using peer-to-peer transmission, provided +you inform other peers where the object code and Corresponding +Source of the work are being offered to the general public at no +charge under subsection 6d.

    +

    A separable portion of the object code, whose source code is excluded +from the Corresponding Source as a System Library, need not be +included in conveying the object code work.

    +

    A "User Product" is either (1) a "consumer product", which means any +tangible personal property which is normally used for personal, family, +or household purposes, or (2) anything designed or sold for incorporation +into a dwelling. In determining whether a product is a consumer product, +doubtful cases shall be resolved in favor of coverage. For a particular +product received by a particular user, "normally used" refers to a +typical or common use of that class of product, regardless of the status +of the particular user or of the way in which the particular user +actually uses, or expects or is expected to use, the product. A product +is a consumer product regardless of whether the product has substantial +commercial, industrial or non-consumer uses, unless such uses represent +the only significant mode of use of the product.

    +

    "Installation Information" for a User Product means any methods, +procedures, authorization keys, or other information required to install +and execute modified versions of a covered work in that User Product from +a modified version of its Corresponding Source. The information must +suffice to ensure that the continued functioning of the modified object +code is in no case prevented or interfered with solely because +modification has been made.

    +

    If you convey an object code work under this section in, or with, or +specifically for use in, a User Product, and the conveying occurs as +part of a transaction in which the right of possession and use of the +User Product is transferred to the recipient in perpetuity or for a +fixed term (regardless of how the transaction is characterized), the +Corresponding Source conveyed under this section must be accompanied +by the Installation Information. But this requirement does not apply +if neither you nor any third party retains the ability to install +modified object code on the User Product (for example, the work has +been installed in ROM).

    +

    The requirement to provide Installation Information does not include a +requirement to continue to provide support service, warranty, or updates +for a work that has been modified or installed by the recipient, or for +the User Product in which it has been modified or installed. Access to a +network may be denied when the modification itself materially and +adversely affects the operation of the network or violates the rules and +protocols for communication across the network.

    +

    Corresponding Source conveyed, and Installation Information provided, +in accord with this section must be in a format that is publicly +documented (and with an implementation available to the public in +source code form), and must require no special password or key for +unpacking, reading or copying.

    +
  14. +
  15. Additional Terms.

    +

    "Additional permissions" are terms that supplement the terms of this +License by making exceptions from one or more of its conditions. +Additional permissions that are applicable to the entire Program shall +be treated as though they were included in this License, to the extent +that they are valid under applicable law. If additional permissions +apply only to part of the Program, that part may be used separately +under those permissions, but the entire Program remains governed by +this License without regard to the additional permissions.

    +

    When you convey a copy of a covered work, you may at your option +remove any additional permissions from that copy, or from any part of +it. (Additional permissions may be written to require their own +removal in certain cases when you modify the work.) You may place +additional permissions on material, added by you to a covered work, +for which you have or can give appropriate copyright permission.

    +

    Notwithstanding any other provision of this License, for material you +add to a covered work, you may (if authorized by the copyright holders of +that material) supplement the terms of this License with terms:

    +

    a) Disclaiming warranty or limiting liability differently from the +terms of sections 15 and 16 of this License; or

    +

    b) Requiring preservation of specified reasonable legal notices or +author attributions in that material or in the Appropriate Legal +Notices displayed by works containing it; or

    +

    c) Prohibiting misrepresentation of the origin of that material, or +requiring that modified versions of such material be marked in +reasonable ways as different from the original version; or

    +

    d) Limiting the use for publicity purposes of names of licensors or +authors of the material; or

    +

    e) Declining to grant rights under trademark law for use of some +trade names, trademarks, or service marks; or

    +

    f) Requiring indemnification of licensors and authors of that +material by anyone who conveys the material (or modified versions of +it) with contractual assumptions of liability to the recipient, for +any liability that these contractual assumptions directly impose on +those licensors and authors.

    +

    All other non-permissive additional terms are considered "further +restrictions" within the meaning of section 10. If the Program as you +received it, or any part of it, contains a notice stating that it is +governed by this License along with a term that is a further +restriction, you may remove that term. If a license document contains +a further restriction but permits relicensing or conveying under this +License, you may add to a covered work material governed by the terms +of that license document, provided that the further restriction does +not survive such relicensing or conveying.

    +

    If you add terms to a covered work in accord with this section, you +must place, in the relevant source files, a statement of the +additional terms that apply to those files, or a notice indicating +where to find the applicable terms.

    +

    Additional terms, permissive or non-permissive, may be stated in the +form of a separately written license, or stated as exceptions; +the above requirements apply either way.

    +
  16. +
  17. Termination.

    +

    You may not propagate or modify a covered work except as expressly +provided under this License. Any attempt otherwise to propagate or +modify it is void, and will automatically terminate your rights under +this License (including any patent licenses granted under the third +paragraph of section 11).

    +

    However, if you cease all violation of this License, then your +license from a particular copyright holder is reinstated (a) +provisionally, unless and until the copyright holder explicitly and +finally terminates your license, and (b) permanently, if the copyright +holder fails to notify you of the violation by some reasonable means +prior to 60 days after the cessation.

    +

    Moreover, your license from a particular copyright holder is +reinstated permanently if the copyright holder notifies you of the +violation by some reasonable means, this is the first time you have +received notice of violation of this License (for any work) from that +copyright holder, and you cure the violation prior to 30 days after +your receipt of the notice.

    +

    Termination of your rights under this section does not terminate the +licenses of parties who have received copies or rights from you under +this License. If your rights have been terminated and not permanently +reinstated, you do not qualify to receive new licenses for the same +material under section 10.

    +
  18. +
  19. Acceptance Not Required for Having Copies.

    +

    You are not required to accept this License in order to receive or +run a copy of the Program. Ancillary propagation of a covered work +occurring solely as a consequence of using peer-to-peer transmission +to receive a copy likewise does not require acceptance. However, +nothing other than this License grants you permission to propagate or +modify any covered work. These actions infringe copyright if you do +not accept this License. Therefore, by modifying or propagating a +covered work, you indicate your acceptance of this License to do so.

    +
  20. +
  21. Automatic Licensing of Downstream Recipients.

    +

    Each time you convey a covered work, the recipient automatically +receives a license from the original licensors, to run, modify and +propagate that work, subject to this License. You are not responsible +for enforcing compliance by third parties with this License.

    +

    An "entity transaction" is a transaction transferring control of an +organization, or substantially all assets of one, or subdividing an +organization, or merging organizations. If propagation of a covered +work results from an entity transaction, each party to that +transaction who receives a copy of the work also receives whatever +licenses to the work the party's predecessor in interest had or could +give under the previous paragraph, plus a right to possession of the +Corresponding Source of the work from the predecessor in interest, if +the predecessor has it or can get it with reasonable efforts.

    +

    You may not impose any further restrictions on the exercise of the +rights granted or affirmed under this License. For example, you may +not impose a license fee, royalty, or other charge for exercise of +rights granted under this License, and you may not initiate litigation +(including a cross-claim or counterclaim in a lawsuit) alleging that +any patent claim is infringed by making, using, selling, offering for +sale, or importing the Program or any portion of it.

    +
  22. +
  23. Patents.

    +

    A "contributor" is a copyright holder who authorizes use under this +License of the Program or a work on which the Program is based. The +work thus licensed is called the contributor's "contributor version".

    +

    A contributor's "essential patent claims" are all patent claims +owned or controlled by the contributor, whether already acquired or +hereafter acquired, that would be infringed by some manner, permitted +by this License, of making, using, or selling its contributor version, +but do not include claims that would be infringed only as a +consequence of further modification of the contributor version. For +purposes of this definition, "control" includes the right to grant +patent sublicenses in a manner consistent with the requirements of +this License.

    +

    Each contributor grants you a non-exclusive, worldwide, royalty-free +patent license under the contributor's essential patent claims, to +make, use, sell, offer for sale, import and otherwise run, modify and +propagate the contents of its contributor version.

    +

    In the following three paragraphs, a "patent license" is any express +agreement or commitment, however denominated, not to enforce a patent +(such as an express permission to practice a patent or covenant not to +sue for patent infringement). To "grant" such a patent license to a +party means to make such an agreement or commitment not to enforce a +patent against the party.

    +

    If you convey a covered work, knowingly relying on a patent license, +and the Corresponding Source of the work is not available for anyone +to copy, free of charge and under the terms of this License, through a +publicly available network server or other readily accessible means, +then you must either (1) cause the Corresponding Source to be so +available, or (2) arrange to deprive yourself of the benefit of the +patent license for this particular work, or (3) arrange, in a manner +consistent with the requirements of this License, to extend the patent +license to downstream recipients. "Knowingly relying" means you have +actual knowledge that, but for the patent license, your conveying the +covered work in a country, or your recipient's use of the covered work +in a country, would infringe one or more identifiable patents in that +country that you have reason to believe are valid.

    +

    If, pursuant to or in connection with a single transaction or +arrangement, you convey, or propagate by procuring conveyance of, a +covered work, and grant a patent license to some of the parties +receiving the covered work authorizing them to use, propagate, modify +or convey a specific copy of the covered work, then the patent license +you grant is automatically extended to all recipients of the covered +work and works based on it.

    +

    A patent license is "discriminatory" if it does not include within +the scope of its coverage, prohibits the exercise of, or is +conditioned on the non-exercise of one or more of the rights that are +specifically granted under this License. You may not convey a covered +work if you are a party to an arrangement with a third party that is +in the business of distributing software, under which you make payment +to the third party based on the extent of your activity of conveying +the work, and under which the third party grants, to any of the +parties who would receive the covered work from you, a discriminatory +patent license (a) in connection with copies of the covered work +conveyed by you (or copies made from those copies), or (b) primarily +for and in connection with specific products or compilations that +contain the covered work, unless you entered into that arrangement, +or that patent license was granted, prior to 28 March 2007.

    +

    Nothing in this License shall be construed as excluding or limiting +any implied license or other defenses to infringement that may +otherwise be available to you under applicable patent law.

    +
  24. +
  25. No Surrender of Others' Freedom.

    +

    If conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot convey a +covered work so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you may +not convey it at all. For example, if you agree to terms that obligate you +to collect a royalty for further conveying from those to whom you convey +the Program, the only way you could satisfy both those terms and this +License would be to refrain entirely from conveying the Program.

    +
  26. +
  27. Use with the GNU Affero General Public License.

    +

    Notwithstanding any other provision of this License, you have +permission to link or combine any covered work with a work licensed +under version 3 of the GNU Affero General Public License into a single +combined work, and to convey the resulting work. The terms of this +License will continue to apply to the part which is the covered work, +but the special requirements of the GNU Affero General Public License, +section 13, concerning interaction through a network will apply to the +combination as such.

    +
  28. +
  29. Revised Versions of this License.

    +

    The Free Software Foundation may publish revised and/or new versions of +the GNU General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns.

    +

    Each version is given a distinguishing version number. If the +Program specifies that a certain numbered version of the GNU General +Public License "or any later version" applies to it, you have the +option of following the terms and conditions either of that numbered +version or of any later version published by the Free Software +Foundation. If the Program does not specify a version number of the +GNU General Public License, you may choose any version ever published +by the Free Software Foundation.

    +

    If the Program specifies that a proxy can decide which future +versions of the GNU General Public License can be used, that proxy's +public statement of acceptance of a version permanently authorizes you +to choose that version for the Program.

    +

    Later license versions may give you additional or different +permissions. However, no additional obligations are imposed on any +author or copyright holder as a result of your choosing to follow a +later version.

    +
  30. +
  31. Disclaimer of Warranty.

    +

    THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY +APPLICABLE LAW. EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT +HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY +OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, +THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM +IS WITH YOU. SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF +ALL NECESSARY SERVICING, REPAIR OR CORRECTION.

    +
  32. +
  33. Limitation of Liability.

    +

    IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS +THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY +GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE +USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF +DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD +PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), +EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF +SUCH DAMAGES.

    +
  34. +
  35. Interpretation of Sections 15 and 16.

    +

    If the disclaimer of warranty and limitation of liability provided +above cannot be given local legal effect according to their terms, +reviewing courts shall apply local law that most closely approximates +an absolute waiver of all civil liability in connection with the +Program, unless a warranty or assumption of liability accompanies a +copy of the Program in return for a fee.

    +
               END OF TERMS AND CONDITIONS
    +
    +  How to Apply These Terms to Your New Programs

    If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms.

    +

    To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +state the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found.

    +

    <one line to give the program's name and a brief idea of what it does.> +Copyright (C)

    +

    This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version.

    +

    This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details.

    +

    You should have received a copy of the GNU General Public License +along with this program. If not, see http://www.gnu.org/licenses/.

    +
  36. +
+

Also add information on how to contact you by electronic and paper mail.

+

If the program does terminal interaction, make it output a short +notice like this when it starts in an interactive mode:

+
<program>  Copyright (C) <year>  <name of author>
+This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+This is free software, and you are welcome to redistribute it
+under certain conditions; type `show c' for details.

The hypothetical commands show w' andshow c' should show the appropriate +parts of the General Public License. Of course, your program's commands +might be different; for a GUI interface, you would use an "about box".

+

You should also get your employer (if you work as a programmer) or school, +if any, to sign a "copyright disclaimer" for the program, if necessary. +For more information on this, and how to apply and follow the GNU GPL, see +http://www.gnu.org/licenses/.

+

The GNU General Public License does not permit incorporating your program +into proprietary programs. If your program is a subroutine library, you +may consider it more useful to permit linking proprietary applications with +the library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. But first, please read +http://www.gnu.org/philosophy/why-not-lgpl.html.

+ +
+
+
+ + + \ No newline at end of file diff --git a/PhysIO/docs/documentation.pdf b/PhysIO/docs/documentation.pdf new file mode 100644 index 00000000..18952419 Binary files /dev/null and b/PhysIO/docs/documentation.pdf differ diff --git a/PhysIO/wikidocs b/PhysIO/wikidocs new file mode 160000 index 00000000..39f3c031 --- /dev/null +++ b/PhysIO/wikidocs @@ -0,0 +1 @@ +Subproject commit 39f3c0313ebc2bc915506374284987226c41b1ed diff --git a/README.txt b/README.md similarity index 70% rename from README.txt rename to README.md index 3eda0eaa..8167d68a 100644 --- a/README.txt +++ b/README.md @@ -1,10 +1,7 @@ -Version 2.7.0.0 - -________________________________________________________________________ - - T A P A S - TNU Algorithms for Psychiatry-Advancing Science. -________________________________________________________________________ +*Version 2.7.3.0* +T A P A S - TNU Algorithms for Psychiatry-Advancing Science. +======================================================================== This file describes the installation and usage of the software. Full details can be found on the TNU website: @@ -21,57 +18,68 @@ particularly concerning the application of neuroimaging and computational modeling to research questions in psychiatry and neurology. Problems that can be addressed by tools in TAPAS presently include: - - Correction of physiological noise in fMRI data. - - Bayesian inference on computational processes from observed behaviour. - - Bayesian mixed-effects inference for classification studies. - - Variational-Bayes Linear Regression. - - Efficient integration of the DCMs using massive parallelization. - - Modeling of eye movements. +- Correction of physiological noise in fMRI data. +- Bayesian inference on computational processes from observed behaviour. +- Bayesian mixed-effects inference for classification studies. +- Variational-Bayes Linear Regression. +- Efficient integration of the DCMs using massive parallelization. +- Modeling of eye movements. TAPAS is written in MATLAB and distributed as open source code under the GNU General Public License (GPL, Version 3). ------------ -INSTRUCTIONS +INSTALLATION ------------ -TAPAS is a collection of toolboxes written in MATLAB (Version R2012b). The key +TAPAS is a collection of toolboxes written in MATLAB (Version R2016b). The key requirement is the installation of MATLAB (produced by The MathWorks, Inc. Natick, MA, USA. http://www.mathworks.com/). -To add the TAPAS directory to the MATLAB path, run the script tapas_init.m in +To add the TAPAS directory to the MATLAB path, run the script `tapas_init.m` in the directory where tapas is installed/extracted. +For the individual toolboxes included in TAPAS, please refer to their +documentation (s.b.) for specific installation instructions. + ------------- DOCUMENTATION ------------- -The TNU website contains all the necessary documentation, user manuals, -publications relevant to the software. - http://www.translationalneuromodeling.org/tapas-documentation/ +For now, the documentation of TAPAS can be found on +- GitHub (latest TAPAS versions) +- the TNU website (previous TAPAS versions) -For queries and discussions please join the mailing list at -http://sympa.ethz.ch/sympa/info/tapas -Bugs in the software can be emailed directly to: tapas-bugs@biomed.ee.ethz.ch +### GitHub ### +- The [Wiki](/../wikis/home) on the + [TAPAS GitHub page](https://github.com/translationalneuromodeling/tapas) contains + the latest documentation about TAPAS. +- If you have a GitHub account, please report bugs or feature requests via the + [Issues](https://github.com/translationalneuromodeling/tapas/issues) page there. ------------------------- -Changes in this version ------------------------- +### TNU-Website ### +- The TNU website contains all the necessary documentation, user manuals, + publications relevant to the previous versions of this software (<= 2.7.0.0) + http://www.translationalneuromodeling.org/tapas-documentation/ -Major Release Notes: +- For queries and discussions please join the mailing list at + http://sympa.ethz.ch/sympa/info/tapas +- Bugs in the software can also be emailed directly to: [tapas-bugs@biomed.ee.ethz.ch] -This release includes the SERIA model as described in +### Cite me ### -http://www.biorxiv.org/content/early/2017/06/08/109090 +Information about citations and current version can be printed from matlab with +the command: -Highlights: +~~~ +tapas_version(1); +~~~ - - The SERIA late race model is implemented. - - The underlying algorithm for computing the likelihood has been updated. - - Includes examples to run the model. +### Current release ### -Detailed release notes can be found in sem/README.md +Information about changes in the current release can be found in the CHANGELOG.md +file. ------- LICENSE diff --git a/misc/tapas_version.m b/misc/tapas_version.m index 4148fc50..fe525d58 100644 --- a/misc/tapas_version.m +++ b/misc/tapas_version.m @@ -14,7 +14,7 @@ % copyright (C) 2017 % -version = {2, 7, 0, 0}; +version = {2, 7, 3, 0}; hash = ''; % In a future implementation if nargin < 1 @@ -26,11 +26,11 @@ fprintf(1, '\n\nVersion %d.%d.%d.%d\n', version{:}); fprintf(1, 'In your citation please include the current version.\n'); fprintf(1, 'Please cite the corresponding masucript according to:\n') - fprintf(1, 'Physio: https://www.ncbi.nlm.nih.gov/pubmed/27832957\n') + fprintf(1, 'PhysIO: https://www.ncbi.nlm.nih.gov/pubmed/27832957\n') fprintf(1, 'HGF: https://www.ncbi.nlm.nih.gov/pubmed/21629826\n') - fprintf(1, ' https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4237059\n'); + fprintf(1, ' https://www.ncbi.nlm.nih.gov/pubmed/25477800\n'); fprintf(1, 'MPDCM: https://www.ncbi.nlm.nih.gov/pubmed/26384541\n'); - fprintf(1, 'SERIA: doi: https://doi.org/10.1101/109090\n'); + fprintf(1, 'SERIA: https://www.ncbi.nlm.nih.gov/pubmed/28767650\n'); end end