function bl_ds=cosmo_meeg_baseline_correct(ds, reference, method)
% correct baseline of MEEG dataset
%
% bl_ds=cosmo_meeg_baseline_correct(ds, reference, method)
%
% Inputs:
% ds MEEG dataset struct with 'time' feature dimension
% reference Either:
% - interval [start, stop], with start and stop in
% seconds (in this case, parameters are estimated in
% this time interval)
% - MEEG dataset struct with 'time' feature dimension,
% with features (in this case, parameters are
% estimated using this dataset)
% method One of:
% - 'absolute' : bl=samples-mu
% - 'relative' : bl=samples/mu
% - 'relchange' : bl=(samples-mu)/mu
% - 'db' : bl=10*log10(samples/mu)
%
% Output:
% bl_ds MEEG dataset struct with the same shape as ds,
% where samples are baseline corrected using reference.
% This is done separately for each combination of values
% in feature dimensions different from 'time', e.g.
% for 'chan' for a timelock dataset and for 'chan' and
% 'freq' combinations for a timefreq dataset
%
%
% Examples:
% % illustrate 'relative' baseline correction
% ds=cosmo_synthetic_dataset('type','timelock','size','small');
% ds=cosmo_slice(ds,1:2); % take first two samples
% ds_rel=cosmo_meeg_baseline_correct(ds,[-.3,-.18],'relative');
% cosmo_disp(ds_rel.samples);
% %|| [ 1 0.572 1 -1.3 1 1.56
% %|| 1 -1.45 1 1.89 1 -0.171 ]
%
% % illustrate 'absolute' baseline correction
% ds=cosmo_synthetic_dataset('type','timelock','size','small');
% ds=cosmo_slice(ds,1:2); % take first two samples
% ds_abs=cosmo_meeg_baseline_correct(ds,[-.3,-.18],'absolute');
% cosmo_disp(ds_abs.samples);
% %|| [ 0 -0.869 0 2.05 0 -0.465
% %|| 0 -1.43 0 1.65 0 -1.36 ]
%
% % illustrate use of another dataset as reference
% ds=cosmo_synthetic_dataset('type','timelock','size','small');
% ds=cosmo_slice(ds,1:2); % take first two samples
% ref=cosmo_synthetic_dataset('type','timelock','size','small');
% ref=cosmo_slice(ref,1:2);
% ds_ref_relch=cosmo_meeg_baseline_correct(ds,ref,'relchange');
% cosmo_disp(ds_ref_relch.samples);
% %|| [ 0.272 -0.272 -7.72 7.72 -0.22 0.22
% %|| -5.41 5.41 -0.309 0.309 1.41 -1.41 ]
%
% # For CoSMoMVPA's copyright information and license terms, #
% # see the COPYING file distributed with CoSMoMVPA. #
if isstruct(reference)
f=@baseline_correct_ds;
elseif isnumeric(reference)
f=@baseline_correct_interval;
else
error('illegal reference: expected dataset struct or vector');
end
baseline_label='time';
check_dataset(ds, baseline_label);
bl_ds=f(ds, reference, baseline_label, method);
function check_dataset(ds, baseline_label)
cosmo_check_dataset(ds);
[dim, index]=cosmo_dim_find(ds,baseline_label,true);
if dim~=2
error(['''%s'' must be feature dimension, found as '...
'sample dimension'], baseline_label);
end
function bl=baseline_correct(samples, mu, method)
assert(size(samples,1)==size(mu,1));
assert(isvector(mu));
assert(numel(size(samples))==2);
switch method
case 'absolute'
bl=bsxfun(@minus,samples,mu);
case 'relchange'
% use absolute and relative, i.e. (samples-mu)/mu
bl=baseline_correct(baseline_correct(samples,mu,'absolute'),...
mu,'relative');
case 'relative'
bl=bsxfun(@rdivide,samples,mu);
case 'vssum'
bl=bsxfun(@rdivide,baseline_correct(samples,mu,'absolute'),...
bsxfun(@plus,samples,mu));
case 'db'
bl=10*log10(baseline_correct(samples,mu,'relative'));
otherwise
error('illegal baseline correction method ''%s''',method);
end
function bl_ds=baseline_correct_interval(ds,interval,baseline_label,method)
if numel(interval)~=2
error('interval must have two values');
end
matcher=@(x) interval(1) <= x & x <= interval(2);
msk=cosmo_dim_match(ds,baseline_label,matcher);
reference=cosmo_slice(ds,msk,2);
bl_ds=baseline_correct_ds(ds,reference,baseline_label,method);
function bl_ds=baseline_correct_ds(ds,reference,baseline_label,method)
check_dataset(reference,baseline_label);
% ensure compatible at sample dimension
[ds1,ds_sa]=first_feature(ds);
[ref1,ref_sa]=first_feature(reference);
if ~isequal(ref_sa, ds_sa)
error('.sa mismatch for input and reference dataset');
end
% split both
ds_split=split_by_other(ds, baseline_label);
reference_split=split_by_other(reference, baseline_label);
n=numel(ds_split);
nref=numel(reference_split);
if n~=nref
error(['Input dataset has different number of dimension '...
'combinations (%d) than reference dataset (%d)'],...
n,nref);
end
for k=1:n
part_ds=ds_split{k};
part_reference=reference_split{k};
fa_ds=cosmo_slice(part_ds.fa,1,2,'struct');
fa_reference=cosmo_slice(part_reference.fa,1,2,'struct');
if ~isequal(rmfield(fa_ds,baseline_label),...
rmfield(fa_reference,baseline_label))
error(['Feature mismatch between '...
'input dataset and reference dataset']);
end
mu=mean(part_reference.samples,2);
ds_split{k}.samples=baseline_correct(part_ds.samples,mu,method);
end
bl_ds=cosmo_stack(ds_split,2);
function ds_split=split_by_other(ds, baseline_label)
other_dims=setdiff(ds.a.fdim.labels,{baseline_label});
ds_split=cosmo_split(ds,other_dims,2);
function [ds1, sa]=first_feature(ds, remove_fields)
ds1=cosmo_slice(ds,1,2);
sa=struct();
if isfield(ds,'sa')
sa=ds.sa;
end
