/trunk/examples/toolboxes/fieldtrip/ft_combineplanar.m
MATLAB | 278 lines | 174 code | 29 blank | 75 comment | 42 complexity | a94584ac2e2dcb5727da8363b6138423 MD5 | raw file
Possible License(s): BSD-3-Clause, AGPL-1.0, LGPL-2.0, GPL-3.0, GPL-2.0, LGPL-2.1, LGPL-3.0, BSD-2-Clause
- function [data] = ft_combineplanar(cfg, data)
- % FT_COMBINEPLANAR computes the planar gradient magnitude over both directions
- % combining the two gradients at each sensor to a single positive-valued number.
- % This can be done for averaged ERFs or TFRs (i.e. powerspectra).
- %
- % Use as
- % [data] = ft_combineplanar(cfg, data)
- % where data contains an averaged planar gradient (either ERF or TFR).
- %
- % In the case of ERFs, the configuration can contain
- % cfg.blc = 'yes' or 'no' (default)
- % cfg.blcwindow = [begin end]
- %
- % After combining the planar data, the planar gradiometer definition does not
- % match the data any more and therefore it is removed from the data. With
- % cfg.combinegrad = 'yes'
- % the function will try to reconstruct the axial gradiometer definition.
- %
- % See also FT_MEGPLANAR
- % Undocumented local options:
- % cfg.baseline
- % cfg.combinemethod
- % cfg.foilim
- % cfg.trials
- % cfg.inputfile = one can specifiy preanalysed saved data as input
- % cfg.outputfile = one can specify output as file to save to disk
- % Copyright (C) 2004, Ole Jensen, Robert Oostenveld
- %
- % This file is part of FieldTrip, see http://www.ru.nl/neuroimaging/fieldtrip
- % for the documentation and details.
- %
- % FieldTrip 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.
- %
- % FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
- %
- % $Id: ft_combineplanar.m 1261 2010-06-22 15:09:23Z roboos $
- fieldtripdefs
- % set the defaults
- if ~isfield(cfg, 'blc'), cfg.blc = 'no'; end
- if ~isfield(cfg, 'foilim'), cfg.foilim = [-inf inf]; end
- if ~isfield(cfg, 'blcwindow'), cfg.blcwindow = [-inf inf]; end
- if ~isfield(cfg, 'combinemethod'), cfg.combinemethod = 'sum'; end
- if ~isfield(cfg, 'trials'), cfg.trials = 'all'; end
- if ~isfield(cfg, 'feedback'), cfg.feedback = 'none'; end
- if ~isfield(cfg, 'inputfile'), cfg.inputfile = []; end
- if ~isfield(cfg, 'outputfile'), cfg.outputfile = []; end
- hasdata = (nargin>1);
- if ~isempty(cfg.inputfile)
- % the input data should be read from file
- if hasdata
- error('cfg.inputfile should not be used in conjunction with giving input data to this function');
- else
- data = loadvar(cfg.inputfile, 'data');
- end
- end
- % check if the input data is valid for this function
- data = checkdata(data, 'datatype', {'raw', 'freq', 'timelock'}, 'feedback', 'yes', 'senstype', {'ctf151_planar', 'ctf275_planar', 'neuromag122', 'neuromag306', 'bti248_planar', 'bti148_planar', 'itab153_planar'});
- cfg = checkconfig(cfg, 'trackconfig', 'on');
- cfg = checkconfig(cfg, 'forbidden', {'combinegrad'});
- cfg = checkconfig(cfg, 'deprecated', {'baseline'});
- if isfield(cfg, 'baseline')
- warning('only supporting cfg.baseline for backwards compatibility, please update your cfg');
- cfg.blc = 'yes';
- cfg.blcwindow = cfg.baseline;
- end
- israw = datatype(data, 'raw');
- isfreq = datatype(data, 'freq');
- istimelock = datatype(data, 'timelock');
- try, dimord = data.dimord; end
- % select trials of interest
- if ~strcmp(cfg.trials, 'all')
- error('trial selection has not been implemented yet') % first fix checkdata (see above)
- end
- % find the combination of horizontal and vertical channels that should be combined
- planar = planarchannelset(data);
- sel_dH = match_str(data.label, planar(:,1)); % indices of the horizontal channels
- sel_dV = match_str(data.label, planar(:,2)); % indices of the vertical channels
- lab_dH = data.label(sel_dH);
- lab_dV = data.label(sel_dV);
- if length(sel_dH)~=length(sel_dV)
- error('not all planar channel combinations are complete')
- end
- % find the other channels that are present in the data
- sel_other = setdiff(1:length(data.label), [sel_dH(:)' sel_dV(:)']);
- lab_other = data.label(sel_other);
- % define the channel names after combining the planar combinations
- % they should be sorted according to the order of the planar channels in the data
- [dum, sel_planar] = match_str(data.label, planar(:,1));
- lab_comb = planar(sel_planar,3);
- % perform baseline correction
- if strcmp(cfg.blc, 'yes')
- if ~(istimelock || israw)
- error('baseline correction is only supported for timelocked or raw input data')
- end
- if ischar(cfg.blcwindow) && strcmp(cfg.blcwindow, 'all')
- cfg.blcwindow = [-inf inf];
- end
- % find the timebins corresponding to the baseline interval
- tbeg = nearest(data.time, cfg.blcwindow(1));
- tend = nearest(data.time, cfg.blcwindow(2));
- cfg.blcwindow(1) = data.time(tbeg);
- cfg.blcwindow(2) = data.time(tend);
- data.avg = blc(data.avg, tbeg, tend);
- end
- if isfreq
-
- switch cfg.combinemethod
- case 'sum'
- if isfield(data, 'powspctrm'),
- % compute the power of each planar channel, by summing the horizontal and vertical gradients
- dimtok = tokenize(dimord,'_');
- catdim = strmatch('chan',dimtok);
- if catdim==1,
- combined = data.powspctrm(sel_dH,:,:,:) + data.powspctrm(sel_dV,:,:,:);
- other = data.powspctrm(sel_other,:,:,:);
- elseif catdim==2,
- combined = data.powspctrm(:,sel_dH,:,:,:) + data.powspctrm(:,sel_dV,:,:,:);
- other = data.powspctrm(:,sel_other,:,:,:);
- else
- error('unsupported dimension order of frequency data');
- end
- data.powspctrm = cat(catdim, combined, other);
- data.label = cat(1, lab_comb(:), lab_other(:));
- else
- error('cfg.combinemethod = ''%s'' only works for frequency data with powspctrm', cfg.combinemethod);
- end
- case 'svd'
- if isfield(data, 'fourierspctrm'),
- fbin = nearest(data.freq, cfg.foilim(1)):nearest(data.freq, cfg.foilim(2));
- Nrpt = size(data.fourierspctrm,1);
- Nsgn = length(sel_dH);
- Nfrq = length(fbin);
- Ntim = size(data.fourierspctrm,4);
- %fourier= complex(zeros(Nrpt,Nsgn,Nfrq,Ntim),zeros(Nrpt,Nsgn,Nfrq,Ntim));
- fourier= zeros(Nrpt,Nsgn,Nfrq,Ntim)+nan;
- progress('init', cfg.feedback, 'computing the svd');
- for j = 1:Nsgn
- progress(j/Nsgn, 'computing the svd of signal %d/%d\n', j, Nsgn);
- for k = 1:Nfrq
- dum = reshape(data.fourierspctrm(:,[sel_dH(j) sel_dV(j)],fbin(k),:), [Nrpt 2 Ntim]);
- dum = permute(dum, [2 3 1]);
- dum = reshape(dum, [2 Ntim*Nrpt]);
- timbin = ~isnan(dum(1,:));
- dum2 = svdfft(dum(:,timbin),1,data.cumtapcnt);
- dum(1,timbin) = dum2;
- dum = reshape(dum(1,:),[Ntim Nrpt]);
- fourier(:,j,k,:) = transpose(dum);
-
- %for m = 1:Ntim
- % dum = data.fourierspctrm(:,[sel_dH(j) sel_dV(j)],fbin(k),m);
- % timbin = find(~isnan(dum(:,1)));
- % [fourier(timbin,j,k,m)] = svdfft(transpose(dum(timbin,:)),1);
- %end
- end
- end
- progress('close');
- other = data.fourierspctrm(:,sel_other,fbin,:);
- data = rmfield(data,'fourierspctrm');
- data.fourierspctrm = cat(2, fourier, other);
- data.label = cat(1, lab_comb(:), lab_other(:));
- data.freq = data.freq(fbin);
- else
- error('cfg.combinemethod = ''%s'' only works for frequency data with fourierspctrm', cfg.combinemethod);
- end
- otherwise
- error('cfg.combinemethod = ''%s'' is not supported for frequency data', cfg.combinemethod);
- end
-
- elseif (israw || istimelock)
- if istimelock,
- % convert timelock to raw
- data = checkdata(data, 'datatype', 'raw', 'feedback', 'yes');
- end
-
- switch cfg.combinemethod
- case 'sum'
- Nrpt = length(data.trial);
- for k = 1:Nrpt
- combined = sqrt(data.trial{k}(sel_dH,:).^2 + data.trial{k}(sel_dV,:).^2);
- other = data.trial{k}(sel_other,:);
- data.trial{k} = [combined; other];
- end
- data.label = cat(1, lab_comb(:), lab_other(:));
- case 'svd'
- Nrpt = length(data.trial);
- Nsgn = length(sel_dH);
- Nsmp = cellfun('size', data.trial, 2);
- Csmp = cumsum([0 Nsmp]);
- % do a 'fixed orientation' across all trials approach here
- % this is different from the frequency case FIXME
- tmpdat = zeros(2, sum(Nsmp));
- for k = 1:Nsgn
- for m = 1:Nrpt
- tmpdat(:, (Csmp(m)+1):Csmp(m+1)) = data.trial{m}([sel_dH(k) sel_dV(k)],:);
- end
- tmpdat2 = abs(svdfft(tmpdat,1));
- tmpdat2 = mat2cell(tmpdat2, 1, Nsmp);
- for m = 1:Nrpt
- if k==1, trial{m} = zeros(Nsgn, Nsmp(m)); end
- trial{m}(k,:) = tmpdat2{m};
- end
- end
- for m = 1:Nrpt
- other = data.trial{m}(sel_other,:);
- trial{m} = [trial{m}; other];
- end
- data.trial = trial;
- data.label = cat(1, lab_comb(:), lab_other(:));
- otherwise
- error('cfg.combinemethod = ''%s'' is not supported for timelocked or raw data', cfg.combinemethod);
- end
-
- if istimelock,
- % convert raw to timelock
- data = checkdata(data, 'datatype', 'timelock', 'feedback', 'yes');
- end
-
- else
- error('unsupported input data');
- end % which datatype
- % remove the fields for which the planar gradient could not be combined
- try, data = rmfield(data, 'crsspctrm'); end
- try, data = rmfield(data, 'labelcmb'); end
- % accessing this field here is needed for the configuration tracking
- % by accessing it once, it will not be removed from the output cfg
- cfg.outputfile;
- % get the output cfg
- cfg = checkconfig(cfg, 'trackconfig', 'off', 'checksize', 'yes');
- % store the configuration of this function call, including that of the previous function call
- try
- % get the full name of the function
- cfg.version.name = mfilename('fullpath');
- catch
- % required for compatibility with Matlab versions prior to release 13 (6.5)
- [st, i] = dbstack;
- cfg.version.name = st(i);
- end
- cfg.version.id = '$Id: ft_combineplanar.m 1261 2010-06-22 15:09:23Z roboos $';
- % remember the configuration details of the input data
- try, cfg.previous = data.cfg; end
- % remember the exact configuration details in the output
- data.cfg = cfg;
- % the output data should be saved to a MATLAB file
- if ~isempty(cfg.outputfile)
- savevar(cfg.outputfile, 'data', data); % use the variable name "data" in the output file
- end