ft_realtime_baseline.m

function ft_realtime_baseline()

% FT_REALTIME_BASELINE is a neurofeedback application based on Hilbert phase estimation.
%
% Use as
%   ft_realtime_baseline()
% with the following configuration options that are coded inside the function
%   cfg.channel    = cell-array, see FT_CHANNELSELECTION (default = 'gui')
%   cfg.foilim     = [Flow Fhigh] (default = [0 120])
%   cfg.blocksize  = number, size of the blocks/chuncks that are processed (default = 1 second)
%   cfg.bufferdata = whether to start on the 'first or 'last' data that is available (default = 'first')
%
% The source of the data is configured as
%   cfg.dataset       = string
% or alternatively to obtain more low-level control as
%   cfg.datafile      = string ( es:  cfg.datafile='buffer://localhost:1972';)
%   cfg.headerfile    = string ( es:  cfg.headerfile='buffer://localhost:1972';)
%   cfg.filename      = string ( es:  cfg.filename = 'buffer://localhost:1972';)
%   cfg.eventfile     = string ( es:  cfg.eventfile = 'buffer://localhost:1972';)
%   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

%Take off the warning message to avoid problems with ATAN2 and hilbert.
warning off all;

cfg = [];
cfg.datafile = 'buffer://localhost:1972';
cfg.headerfile = 'buffer://localhost:1972';
cfg.filename = 'buffer://localhost:1972';

% 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, 'blocksize'),      cfg.blocksize = 1;        end % in seconds
if ~isfield(cfg, 'overlap'),        cfg.overlap = 0.5;        end % in seconds
if ~isfield(cfg, 'channel'),        cfg.channel = {'all', '-Fp1', '-Fp2',...
                                    '-F7','-Fz','-F8','-FT7', '-FC3', '-FCz',...
                                    '-FC4', '-FT8', '-T3', '-C3', '-Cz', '-C4','-T4','-TP7',...
                                    '-CP3', '-CPz', '-CP4','-TP8','-A1', '-T5', '-Pz',...
                                    '-T6','-A2', '-O1', '-Oz', '-O2'};
end                                                               % This will select F3 F4 P3 P4 for EEG
                                                                  % X1 X2 X3 X4 for EMG
                                                                  % X5 X6 X7 X8 for EOG
if ~isfield(cfg, 'foilim'),         cfg.foilim = [0.1 100];   end
if ~isfield(cfg, 'bufferdata'),     cfg.bufferdata='last';   end

% translate dataset into datafile+headerfile
cfg = checkconfig(cfg, 'dataset2files', 'yes');
cfg = checkconfig(cfg, 'required', {'datafile' 'headerfile'});

% ensure that the persistent variables related to caching are cleared
clear ft_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
lab=['X1 '; 'X2 '; 'Fp1'; 'Fp2'; 'X3 '; 'X4 '; 'F7 '; 'F3 '; 'Fz '; 'F4 ';...
    'F8 ';'FT7'; 'FC3'; 'FCz'; 'FC4'; 'FT8'; 'T3 '; 'C3 '; 'Cz '; 'C4 ';...
    'T4 ';'TP7'; 'CP3'; 'CPz'; 'CP4'; 'TP8'; 'A1 '; 'T5 '; 'P3 '; 'Pz ';...
    'P4 '; 'T6 ';'A2 '; 'O1 '; 'Oz '; 'O2 '; 'X5 '; 'X6 '; 'X7 '; 'X8 '];
label=cellstr(lab);
cfg.channel = ft_channelselection(cfg.channel, label);
chanindx    = match_str(label, cfg.channel);
nchan       = length(chanindx);
if nchan==0
  error('no channels were selected');
end

% determine the size of blocks to process
blocksize = cfg.blocksize * hdr.Fs;
overlap   = cfg.overlap * hdr.Fs;
prevSample  = 0;
count       = 0;

% Create arrays that contains the rhos found
nummaxarray = 240; % the step is 500ms so 240 are 2 minutes
i=1;
vlmean=zeros(1,nummaxarray);
vrmean=zeros(1,nummaxarray);
frontmeans=zeros(1,nummaxarray);
neckmeans=zeros(1,nummaxarray);

% This is used to plot the screen for the subject
fullscreen = get(0,'ScreenSize');
fig1 = figure('NumberTitle','off', ...
    'MenuBar','none', ...
    'Units','pixels', ...
    'Position',[0 0 fullscreen(3) fullscreen(4)]);

% plot the feedback on the second monitor
% set(gcf,'position',[1025,1,1024,768]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% this is the general BCI loop where realtime incoming data is handled
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
while true

    % determine number of samples available in buffer
    hdr = ft_read_header(cfg.headerfile, 'cache', true);

    % see whether new samples are available
    newsamples = (hdr.nSamples*hdr.nTrials-prevSample);

    if newsamples>=blocksize

        % determine the samples to process
        if strcmp(cfg.bufferdata, 'last') && count==0
            begsample  = hdr.nSamples*hdr.nTrials - blocksize + 1;
            endsample  = hdr.nSamples*hdr.nTrials;
        elseif strcmp(cfg.bufferdata, 'last')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        elseif strcmp(cfg.bufferdata, 'first')
          begsample  = prevSample+1;
          endsample  = prevSample+blocksize ;
        else
            error('unsupported value for cfg.bufferdata');
        end

        % this allows overlapping data segments
        if overlap && (begsample>overlap) %#ok`<BDLGI>`
          begsample = begsample - overlap;
          endsample = endsample - overlap;
        end

        % remember up to where the data was read
        prevSample  = endsample;
        count       = count + 1;
        % fprintf('processing segment %d from sample %d to %d\n', count, begsample, endsample);

        % 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 hilbert phase sinchronisation from the data %
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        % put the data in a fieldtrip-like raw structure
        data.trial{1} = dat;
        data.time{1}  = offset2time(begsample, hdr.Fs, endsample-begsample+1);
        data.label    = hdr.label(chanindx);
        data.hdr      = hdr;
        data.fsample  = hdr.Fs;

        % correction of EOG based on algoritm fro
        % Author: German Gomez-Herrero
        %         german.gomezherrero@ieee.org
        %         http://www.cs.tut.fi/~gomezher/index.htm
        %         Institute of Signal Processing
        %         Tampere University of Technology, 2007
        % Reference
        % [1] P. He et al., Med. Biol. Comput. 42 (2004), 407-412
        % [2] S. Haykin. Adaptive Filter Theory, (1996), Prentice Hall

        data.trial{2}(9,:)=data.trial{1}(9,:)-data.trial{1}(10,:);
        data.trial{2}(11,:)=data.trial{1}(11,:)-data.trial{1}(12,:);
        opt.refdata=[data.trial{2}(9,:);data.trial{2}(11,:)];
        [data.trial{3}] = crls_regression(data.trial{1}(5:8,:), opt);

        % Build a FIR filter for EMG correction
        N    = 150;      % Order
        gammaband = [35 45];
        emgband = [60 80];
        emgfnband = [60 499];
        flag = 'scale';  % Sampling Flag
        Beta = 0.9;      % Window Parameter
        win = kaiser(N+1, Beta);

        % Correction between EMG and EEG based on Sheer D.E. "Biofeedback training
        % of 40-Hz eeg and behavior", pp. 325-362, on Behavior and
        % brain electrical activity (1975), Plenum Press. New York

        gammafilter  = fir1(N, gammaband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltgamma1 = filtfilt(gammafilter,1,data.trial{1}(5,:));
        datfiltgamma2 = filtfilt(gammafilter,1,data.trial{1}(6,:));
        datfiltgamma3 = filtfilt(gammafilter,1,data.trial{1}(7,:));
        datfiltgamma4 = filtfilt(gammafilter,1,data.trial{1}(8,:));
        datfiltgamma = [datfiltgamma1; datfiltgamma2; datfiltgamma3; datfiltgamma4];

        emgfilter = fir1(N, emgband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemg1 = filtfilt(emgfilter,1,data.trial{1}(5,:));
        datfiltemg2 = filtfilt(emgfilter,1,data.trial{1}(6,:));
        datfiltemg3 = filtfilt(emgfilter,1,data.trial{1}(7,:));
        datfiltemg4 = filtfilt(emgfilter,1,data.trial{1}(8,:));
        datfiltemg = [datfiltemg1; datfiltemg2; datfiltemg3; datfiltemg4];

        datfiltemgsqr = datfiltemg.^2;
        datfiltgammasqr = datfiltgamma.^2;
        datfiltcrossqr = (datfiltemg.*datfiltgamma).^2;

        correction = datfiltgammasqr-(datfiltcrossqr./datfiltemgsqr);

        data.trial{3}(5:8,:) = datfiltgamma - correction;

        % Find the EMG on forehead and neck
        data.trial{2}(1,:)=data.trial{1}(1,:)-data.trial{1}(2,:);   % Frontal electrods
        data.trial{2}(3,:)=data.trial{1}(3,:)-data.trial{1}(4,:);   % Neck electrods
        emgfnfilter = fir1(N, emgfnband/(hdr.Fs/2), 'bandpass', win, flag);
        datfiltemgf = filter(emgfnfilter,1,data.trial{2}(1,:));
        datfiltemgn = filter(emgfnfilter,1,data.trial{2}(3,:));
        datfiltemgf = abs(datfiltemgf);
        datfiltemgn = abs(datfiltemgn);
        extrMaxValuef = datfiltemgf(find(diff(sign(diff(datfiltemgf)))==-2)+1);
        extrMaxValuen = datfiltemgn(find(diff(sign(diff(datfiltemgn)))==-2)+1);
        extrMaxIndexf =   find(diff(sign(diff(datfiltemgf)))==-2)+1;
        extrMaxIndexn =   find(diff(sign(diff(datfiltemgn)))==-2)+1;
        upf = extrMaxValuef;
        upn = extrMaxValuen;
        upf_t = data.time{1}(extrMaxIndexf);
        upn_t = data.time{1}(extrMaxIndexn);
        upf = interp1(upf_t,upf,data.time{1},'linear');
        upn = interp1(upn_t,upn,data.time{1},'linear');
        emgfmean = nanmean (upf'); %#ok`<UDIM>`
        emgnmean = nanmean (upn'); %#ok`<UDIM>`
        % plot(data.time{1},upf,'r')

        % Istantaneous (proto)phase difference found via Hilbert
        % Based on Pikovsky, A. R. (2001). Synchronization. A Universal
        % Concept In Nonlinear Sciences. Cambridge: Cambridge University
        % Press, pag. 368 A2.7

        % crate the data needed for phase coherence index
        chan1=data.trial{3}(5,:); % F3
        chan2=data.trial{3}(6,:); % F4
        chan3=data.trial{3}(7,:); % P3
        chan4=data.trial{3}(8,:); % P4
        chan1h = hilbert(chan1);
        chan2h = hilbert(chan2);
        chan3h = hilbert(chan3);
        chan4h = hilbert(chan4);
        chan1hi = imag(chan1h);
        chan2hi = imag(chan2h);
        chan3hi = imag(chan3h);
        chan4hi = imag(chan4h);

        % find the istantaneous left hemisphere (proto)phase difference
        phil = atan2(((chan1hi .* chan3)-(chan1 .* chan3hi)),...
                     ((chan1 .* chan3)+(chan1hi .* chan3hi)));

        % find the istantaneous right hemisphere (proto)phase difference
        phir = atan2(((chan2hi .* chan4)-(chan2 .* chan4hi)),...
                     ((chan2 .* chan4)+(chan2hi .* chan4hi)));

        % find the right hemisphere synchronization index
        sumsinr = sum(sin(phir))/blocksize;
        sumcosr = sum(cos(phir))/blocksize;
        rhor = sqrt(sumsinr.^2 + sumcosr.^2);

        % find the left hemisphere synchronization index
        sumsinl = sum(sin(phil))/blocksize;
        sumcosl = sum(cos(phil))/blocksize;
        rhol = sqrt(sumsinl.^2 + sumcosl.^2);

        % update the array for the mean
        vlmean(i)=rhol;
        vrmean(i)=rhor;
        frontmeans(i)=emgfmean;
        neckmeans(i)=emgnmean;

        i=i+1;

    end % if enough new samples

    % screen for the subject
    clf;
    set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[-50 50],'YColor',[0.8 0.8 0.8],...
        'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
        'DataAspectRatio',[0.2 0.1 3],'OuterPosition', [-0.0175 0.185 1 0.605],...
        'Color',[0.8 0.8 0.8],'CLim', [-50 50]);
    annotation(fig1,'textbox',[0.18 0.35 0.6292 0.1929],...
        'String',{'Rilassati mantenendo gli occhi aperti'},...
        'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
        'Color',[1 0 0]);
    drawnow ;

    if count > (nummaxarray-1)
        close all hidden;
        break;
    end

end % while true

% print the mean baseline
baselinel = mean (vlmean);
fprintf ('mean baseline left= %d \n', baselinel);
baseliner = mean (vrmean);
fprintf ('mean baseline right= %d \n', baseliner);
baselinef = mean (frontmeans);
fprintf ('mean baseline forehead= %d \n', baselinef);
baselinen = mean (neckmeans);
fprintf ('mean baseline neck= %d \n', baselinen);
stf = std(frontmeans);
fprintf ('Std baseline front= %d \n', stf);
stn = std(neckmeans);
fprintf ('Std baseline neck= %d \n', stn);

save means.mat;   % save to the working folder - (edit jonaweber 2010)