Differences

This shows you the differences between two versions of the page.

Link to this comparison view

example:ft_realtime_selectiveaverage [2017/08/17 11:21] (current)
Line 1: Line 1:
 +{{tag>​example realtime}}
 +
 +====== Example real-time selective average ======
 +
 +===== Flowchart =====
 +
 +{{:​example:​realtime:​realtime_selectiveaverage.png?​400}}
 +
 +===== Matlab code =====
 +
 +<code matlab>
 +function ft_realtime_selectiveaverage(cfg)
 +
 +% FT_REALTIME_SELECTIVEAVERAGE is an example realtime application for online
 +% averaging of the data. It should work both for EEG and MEG.
 +%
 +% Use as
 +%   ​ft_realtime_selectiveaverage(cfg)
 +% with the following configuration options
 +%   ​cfg.channel ​   = cell-array, see FT_CHANNELSELECTION (default = '​all'​)
 +%   ​cfg.trialfun ​  = string with the trial function
 +%
 +% The source of the data is configured as
 +%   ​cfg.dataset ​      = string
 +% or alternatively to obtain more low-level control as
 +%   ​cfg.datafile ​     = string
 +%   ​cfg.headerfile ​   = string
 +%   ​cfg.eventfile ​    = string
 +%   ​cfg.dataformat ​   = string, default is determined automatic
 +%   ​cfg.headerformat ​ = string, default is determined automatic
 +%   ​cfg.eventformat ​  = string, default is determined automatic
 +%
 +% To stop the realtime function, you have to press Ctrl-C
 +
 +% Copyright (C) 2008, Robert Oostenveld
 +%
 +% Subversion does not use the Log keyword, use 'svn log <​filename>'​ or 'svn -v log | less' to get detailled information
 +
 +% set the default configuration options
 +if ~isfield(cfg,​ '​dataformat'​), ​    ​cfg.dataformat = [];      end % default is detected automatically
 +if ~isfield(cfg,​ '​headerformat'​), ​  ​cfg.headerformat = [];    end % default is detected automatically
 +if ~isfield(cfg,​ '​eventformat'​), ​   cfg.eventformat = [];     end % default is detected automatically
 +if ~isfield(cfg,​ '​channel'​), ​       cfg.channel = '​all'; ​     end
 +if ~isfield(cfg,​ '​bufferdata'​), ​    ​cfg.bufferdata = '​last'; ​ end % first or last
 +
 +% translate dataset into datafile+headerfile
 +cfg = ft_checkconfig(cfg,​ '​dataset2files',​ '​yes'​);​
 +cfg = ft_checkconfig(cfg,​ '​required',​ {'​datafile'​ '​headerfile'​});​
 +
 +% ensure that the persistent variables related to caching are cleared
 +clear read_header
 +% start by reading the header from the realtime buffer
 +hdr = ft_read_header(cfg.headerfile,​ '​cache',​ true);
 +
 +% define a subset of channels for reading
 +cfg.channel = channelselection(cfg.channel,​ hdr.label);
 +chanindx ​   = match_str(hdr.label,​ cfg.channel);​
 +nchan       = length(chanindx);​
 +if nchan==0
 +  error('​no channels were selected'​);​
 +end
 +
 +prevSample = 0;
 +count      = 0;
 +
 +% initialize the timelock cell-array, each cell will hold the average in one condition
 +timelock = {};
 +
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +% this is the general BCI loop where realtime incoming data is handled
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +while true
 +
 +  % determine latest header and event information
 +  event     = read_event(cfg.dataset,​ '​minsample',​ prevSample+1); ​ % only consider events that are later than the data processed sofar
 +  hdr       = read_header(cfg.dataset,​ '​cache',​ true); ​            % the trialfun might want to use this, but it is not required
 +  cfg.event = event; ​                                              % store it in the configuration,​ so that it can be passed on to the trialfun
 +  cfg.hdr ​  = hdr;                                                 % store it in the configuration,​ so that it can be passed on to the trialfun
 +
 +  % evaluate the trialfun, note that the trialfun should not re-read the events and header
 +  fprintf('​evaluating ''​%s''​ based on %d events\n',​ cfg.trialfun,​ length(event));​
 +  trl = feval(cfg.trialfun,​ cfg);
 +
 +  % the code below assumes that the 4th column of the trl matrix contains the condition index
 +  % set the default condition to one if no condition index was given
 +  if size(trl,​2)<​4
 +    trl(:,4) = 1;
 +  end
 +
 +  fprintf('​processing %d trials\n',​ size(trl,​1));​
 +
 +  for trllop=1:​size(trl,​1)
 +
 +    begsample = trl(trllop,​1);​
 +    endsample = trl(trllop,​2);​
 +    offset ​   = trl(trllop,​3);​
 +    condition = trl(trllop,​4);​
 +
 +    % remember up to where the data was read
 +    prevSample ​ = endsample;
 +    count       = count + 1;
 +    fprintf('​processing segment %d from sample %d to %d, condition = %d\n', count, begsample, endsample, condition);
 +
 +    % read data segment from buffer
 +    dat = ft_read_data(cfg.datafile,​ '​header',​ hdr, '​begsample',​ begsample, '​endsample',​ endsample, '​chanindx',​ chanindx, '​checkboundary',​ false);
 +
 +    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +    % from here onward it is specific to the processing of the data
 +    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
 +    % put the data in a fieldtrip-like raw structure
 +    data.trial{1} = dat;
 +    data.time{1} ​ = offset2time(offset,​ hdr.Fs, endsample-begsample+1);​
 +    data.label ​   = hdr.label(chanindx);​
 +    data.hdr ​     = hdr;
 +    data.fsample ​ = hdr.Fs;
 +
 +    % apply some preprocessing options
 +    data.trial{1} = preproc_baselinecorrect(data.trial{1});​
 +
 +    if length(timelock)<​condition || isempty(timelock{condition})
 +      % this is the first occurence of this condition, initialize an empty timelock structure
 +      timelock{condition}.label ​  = data.label;
 +      timelock{condition}.time ​   = data.time{1};​
 +      timelock{condition}.avg ​    = [];
 +      timelock{condition}.var ​    = [];
 +      timelock{condition}.dimord ​ = '​chan_time';​
 +      nchans ​  = size(data.trial{1},​ 1);
 +      nsamples = size(data.trial{1},​ 2);
 +      % the following elements are for the cumulative computation
 +      timelock{condition}.n ​  = 0;                          % number of trials
 +      timelock{condition}.s ​  = zeros(nchans,​ nsamples); ​   % sum
 +      timelock{condition}.ss ​ = zeros(nchans,​ nsamples); ​   % sum of squares
 +    end
 +
 +    % add the new data to the accumulated data
 +    timelock{condition}.n ​ = timelock{condition}.n ​ + 1;
 +    timelock{condition}.s ​ = timelock{condition}.s ​ + data.trial{1};​
 +    timelock{condition}.ss = timelock{condition}.ss + data.trial{1}.^2;​
 +
 +    % compute the average and variance on the fly
 +    timelock{condition}.avg = timelock{condition}.s ./ timelock{condition}.n;​
 +    timelock{condition}.var = (timelock{condition}.ss - (timelock{condition}.s.^2)./​timelock{condition}.n) ./ (timelock{condition}.n-1);​
 +
 +    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +    % from here onward additional processing of the selective averages could be done
 +    % as an example here the ERP of each condition is plotted in its own figure
 +    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
 +    % compute the t-score versus zero by dividing the average by the standard error of mean
 +    tscore = timelock{condition}.avg ./ (sqrt(timelock{condition}.var)./​(timelock{condition}.n - 1));
 +
 +    figure(condition)
 +    plot(timelock{condition}.time,​ tscore);
 +    title(sprintf('​condition %d, ntrials = %d', condition, timelock{condition}.n));​
 +
 +    % force matlab to redraw the figure
 +    drawnow
 +
 +  end % looping over new trials
 +end % while true
 +</​code>​