Realtime neurofeedback application based on Hilbert phase estimation
This page describes the neurofeedback application that was developed by Marco Rotonda for his MSc thesis “Controllo volontario della sincronizzazione di fase intraemisferica nella banda gamma eeg mediante neurofeedback” (Voluntary control of intrahemispheric phase sychronization in the EEG gamma band using neurofeedback).
The thesis (in Italian) which describes all details and the results is available here. The signal analysis done are available here
Neurofeedback, as the word suggest, is feeding back to the subject it’s neuronal activity analyzed by it’s EEG activity in (almost) real time. In this training we give back to the subjest a visual stimulation (vertical bars) when it’s gamma phase synchronization between F3-P3 and F4-P4 increased over the baseline (bars up) or decreased (bars down). Here is the feedback
If you scroll the script you can see that this image is a particular point of view of a 3D visualization. I you change the POV you could see all the data from the last 3 minutes.
Procedure.
Get the FieldtripBufferDemo from the workshop_bci2000 folder on the fieldtrip ftp server and . Start BCI2000 and modifiy the config. The script works originally with 40 channels (see the variable 'lab' in the realtime_hilbert script), the samplingrate used is 1000Hz. Use a blocksize of 1000 samples. Ensure BCI2000 runs for more than 2 minutes in the application tab. Set config and start BCI2000. Start Matlab and in the shell type realtime_hilbert.
Now it will start the function realtime_baseline that for 2 minutes will record the subject baseline and store the results for realtime_hilbert. Basically there are 12 channels: 4 for EEG (I used F3-P3-F4-P4), 4 for the EOG, 4 for the EMG (front and neck). In this script, those are selected from the forty originally recorded channels you can find in the variable 'lab'. The algorithm will care to take away both EOG and EMG artifacts.
Source code for realtime_hilbert
function realtime_hilbert() % realtime_hilbert is a neurofeedback application based on hilbert phase estimation. % % Use as % realtime_hilbert() % with the following configuration options % 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 from: % 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 % References: % [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
Source code for realtime_baseline
function realtime_baseline() % realtime_hilbert is a neurofeedback application based on hilbert phase estimation. % % Use as % realtime_hilbert() % with the following configuration options % 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 from: % 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 % References: % [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)
127.0.0.1
