ft_realtime_hilbert.m

function ft_realtime_hilbert()

% FT_REALTIME_HILBERT is a neurofeedback application based on Hilbert phase estimation.
%
% Use as
%   ft_realtime_hilbert()
% 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;

% This is to save subject data
starttime= DATESTR(now, 30);
subjname = input('Insert the subject name.>>>>>', 's');
trialdata=strcat(starttime,subjname);

realtime_baseline(); %This it will launch the baseline script

load means;  % get means from file created by realtime_baseline

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 = ft_checkconfig(cfg, 'dataset2files', 'yes');
cfg = ft_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
vectordim = 360;         % vector dimention for the lasts 3 minutes of data
vl=zeros(vectordim,1);  % the visual vector of data for left synchrony
vr=zeros(vectordim,1);  % the visual vector of data for right synchrony
vl1=vl(1);              % the first element of the vector for the visual feedback
vr1=vr(1);              % the first element of the vector for the visual feedback

% This is the ratio above which accept the coherence
meanl = baselinel;           % It has to be calculated from the baseline
meanr = baseliner;           % It has to be calculated from the baseline

% This is the mean above which no EMG feedback has given
meanf = baselinef+(stf*2);             % It has to be calculated from the baseline
meann = baselinen+(stn*2);             % It has to be calculated from the baseline

% This is used to plot the feedback step
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]);

% fig2=figure;

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

    % Create a matrix that contains the rhos found
    M = [vl vr]; %#ok`<NASGU>`

    % Create a matrix for the visual feedback
    K = [vl1 vr1];

    % 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);

        % Give the visual feedback
        if (meanf>emgfmean) && (meann>emgnmean)
            clf;
            visual = bar3(K,0.3);
            view([-90 0]);
            grid off;
            shading interp;
            for i = 1:length(visual)
                zdata = get(visual(i),'Zdata');
                set(visual(i),'Cdata',zdata,'EdgeColor','none')
                colormap hot;
            end
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[-1 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'Color',[0.8 0.8 0.8],'CLim', [-1 1]);

            % Create colorbar
            colorbar([0.5 0.148 0.02 0.73],'ZColor',[0.8 0.8 0.8],'YTick',[],...
                'YColor',[0.8 0.8 0.8],'XColor',[0.8 0.8 0.8]);

            % Create textboxes
            annotation('textbox',[0.496 0.89 0.03 0.04],'String',{'+'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.50 0.11 0.02 0.04],'String',{'-'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.25 0.49 0.05 0.07],'String',{'Emisfero','sinistro'},...
                'HorizontalAlignment','center','FontSize',14,'FitBoxToText','off','EdgeColor','none');
            annotation('textbox',[0.735 0.49 0.05 0.07],'String',{'Emisfero','destro'},...
                'HorizontalAlignment','center','FontSize',14,'FitBoxToText','off','EdgeColor','none');

            %create a copy of the plot for the experimenter
            % h1=gcf;
            % h2=figure;
            % objects=allchild(h1);
            % fig2=copyobj(get(h1,'children'),h2);

            % force MATLAB to update the figure
            drawnow ;

            % add the step to the array for the feedback and upgrade the M
            % for the final performance plot

            vl = vl([end 1:end-1]);
            vl(1)=rhol-meanl;
            vl1=vl(1);

            vr = vr([end 1:end-1]);
            vr(1)=rhor-meanr;
            vr1=vr(1);

        elseif  emgfmean>meanf
            clf;
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[0 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'OuterPosition', [-0.0175 0.185 1 0.605],...
                'Color',[0.8 0.8 0.8],'CLim', [-50 50]);
            annotation(fig1,'textbox',[0.20 0.35 0.6292 0.1929],...
                'String',{'Rilassa i muscoli della fronte'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
                'Color',[1 0 0]);
            drawnow ;

        elseif  emgnmean>meann
            clf;
            set(gca,'ZColor',[0.8 0.8 0.8],'Zlim',[0 1],'YColor',[0.8 0.8 0.8],...
                'Ylim',[0 3],'XColor',[0.8 0.8 0.8],'Xlim',[0.85 1.15],...
                'OuterPosition', [-0.0175 0.185 1 0.605],...
                'Color',[0.8 0.8 0.8],'CLim', [-1 1]);
            annotation(fig1,'textbox',[0.20 0.35 0.6292 0.1929],...
                'String',{'Rilassa i muscoli del collo'},...
                'HorizontalAlignment','center','FontSize',20,'FitBoxToText','off','EdgeColor','none',...
                'Color',[1 0 0]);
            drawnow ;

        end % end feedbacks

    end % if enough new samples

    if count > vectordim
        cd 'C:\TEST\hilbert\risultati';
        save (trialdata);
        save hilbert.mat;
        close all hidden;
        cd 'C:\TEST';
        break;
    end

end % while true