Code covered by the BSD License  

Highlights from
Brain-Machine Interface (BMI) based on Electroencephalography (EEG)

image thumbnail

Brain-Machine Interface (BMI) based on Electroencephalography (EEG)

by

 

Real-Time Discrete Wavelet Transform and ANFIS classifier for Brain-Machine Interface based on EEG

bmi_three_channels.m
function varargout = bmi_three_channels(varargin)
% BMI_THREE_CHANNELS Application M-file for bmi_three_channels.fig

% The explanation of the project methodology and results is presented on:
% http://www.youtube.com/watch?v=4IodfA_fHUM

% Moreover, the signal processing algorithm and pattern recognition system
% are presented in the IEEE conference publication: 
% Eduardo Lpez-Arce Vivas, Alejandro Garca-Gonzlez, Ivn Figueroa, 
% and Rita Fuentes. Discrete Wavelet Transform and ANFIS Classifier for 
% Brain-Machine Interface based on EEG. International Conference on Human 
% System Interaction, 2013.
% http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6577814&tag=1

% BMI_THREE_CHANNELS is an aplication developed for three different
% purposes: EEG and ECG three channels monitoring, and for a Brain-Machine
% Interface (BMI) based on one EEG bipolar connection: O1-P3.
% The biosignals were acquired using self-made intrumentation and
% converting from analog data to digital data with the DAQ: USB-6009.

% For the BMI application:
% The Brain-Machine Interface consists of a robotic hand controlled by 
% two servomotors. The robotic hand closes and opens 
% whenever a subject closes or opens its eyes, respectively. The EEG signal
% was acquired using the NI USB-6009 Hardware and the Matlab Data 
% Acquisition Toolbox. The signal was processed by a Discrete-Wavelet 
% Transform algorithm to analyze the data on-line with help of the Wavelet 
% Toolbox. An ANFIS was used to classify the data previously processed. 
% The ANFIS was created and trained using the Fuzzy Logic Toolbox. A 
% Graphical User Interface was created to plot and save the EEG signal, 
% analyze the data offline, and compare the output of the ANFIS with the 
% actual state of the subject's eyes.

% AUTHOR: EDUARDO LPEZ ARCE VIVAS (eduardo.lopezarce@gmail.com)

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',          mfilename, ...
                   'gui_Singleton',     gui_Singleton, ...
                   'gui_OpeningFcn',    @bmi_three_channels_OpeningFcn, ...
                   'gui_OutputFcn',     @bmi_three_channels_OutputFcn, ...
                   'gui_LayoutFcn',     [], ...
                   'gui_Callback',      []);
if nargin && ischar(varargin{1})
   gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% clc;
% End initialization code - DO NOT EDIT

% --- Executes just before bmi_three_channels is made visible.
function bmi_three_channels_OpeningFcn(hObject, ~, handles, varargin)
% Choose default command line output for bmi_three_channels
handles.output = hObject;
% Declare global variables to be used
global vid hImage register

% Set textbox background: white
set(handles.text11,'String','','BackgroundColor','white');
% Variables initialization
register=1; % Register = 1 to divide webcam in two images (eyes and hand)
eeg_p=0;

% Set main webcam properties and initialize video
guidata(hObject, handles);
axes(handles.webcam)
vid = videoinput('winvideo');
vidRes = get(vid, 'VideoResolution');
nBands = get(vid, 'NumberOfBands');
hImage = image( zeros(vidRes(2), vidRes(1), nBands) );
preview(vid,hImage); 
axis off;

% Second video is done with the same webcam
% Get snapshot for static image in the second video box
axes(handles.webcam2);
image(getsnapshot(vid));
set(gca,'xlim',[300,900],'ylim',[600 1000]);
axis off;
pause(0.1);
clc;

function plot_button_Callback(hObject, ~, handles, varargin)
% Declare and initialize global variables.
global caso register fs analyze data estado1 estado2 vid temp_adc p1 p2 p3 output total_pt proc_t
estado1=0;
estado2=0;
analyze=0;
temp_adc=0;
contador=1;

% Get value when start button has been pressed
while (get (hObject, 'Value')&&contador) 
    
% Get user input from GUI
fs = str2double(get(handles.f1_input,'String'));
temp = str2double(get(handles.f2_input,'String'));
duration = str2double(get(handles.t_input,'String'));
caso=get(handles.menu_registros, 'Value');
np=fs*temp;

% NI DAQ properties for analog input channels
ni=daq.createSession('ni');         % Create ni session
channels_i=0:2;                     % Declare channels (0 to 2)
ni.addAnalogInputChannel('Dev1',channels_i,'Voltage');  % Create input obj
set(ni.Channels,'InputType','Differential');
set(ni.Channels,'Range',[-5 5]);    % Voltage range

% Set the sample rate and samples per trigger
ni.Rate = fs;                       % Sampling frequency 
ni.NumberOfScans = fs*duration+1;   % Total of samples to acquire
data=[];

% Add listener to get data when it is available in the buffer
lh = ni.addlistener('DataAvailable',...
        @(src,event) adc_interrupt(event));
% Set number of data to start listener event
ni.NotifyWhenDataAvailableExceeds = 10;
ni.startBackground();

% NI DAQ properties for digital output channels
ni_o=daq.createSession('ni');
channels_o=0:1;                     % Select output channels (0 to 1)
ni_o.addAnalogOutputChannel('Dev1',channels_o,'Voltage'); % Create output obj
figure(1);
close 1;
% Variables initialization
max_cD4=0;          % Maximum value of wavelet coefficient at level 4
cD4=0;              % Wavelet coefficient at level 4
ojos=0;             % Classifier output (ojos = 0, eyes closed; 1, open)
casos=1;            % Temporal variable
proc_t=0;           % Processing time temporal variable     
times=0;            % Temporal variable
total_pt=0;         % Processing time
output=[];          % Data vector
estado2=3.3;        % Output voltage
temp_adc=0;         % Adc value
window_a=0.075;     % Wavelet window value in seconds   
data_plot=[];       % Data temporal variable

% Read ANFIS system, after training it off-line
fuzzy=readfis('ANFIS_EEG_OR');

switch(caso)
    case 1 %EEG
%       Hide one channel from GUI and adjust the other two (EEG and EOG)
        cla(handles.chan2,'reset');
        axes(handles.chan2);
        axis off;
        set(handles.chan1,'Position',[0.088 0.550 0.671 0.357]);
        set(handles.chan3,'Position',[0.088 0.077 0.671 0.357]);
%       Get data until stop time or until user presses the stop button
        while (ni.IsRunning&&get (hObject, 'Value'))
%           Save data acquired with the listener-event function:
%           adc_interrupt
            data_temp=data;
            data_plot=data;
%           Start stopwatch timer
            tic;
%           Start calibration
            if (~temp_adc)&&(length(data_temp(:,1))>=(fs+1))
%               Analyze last second of the signal
                ult_seg=data_temp((length(data_temp(:,1))-fs):length(data_temp(:,1)),...
                    1);
%               Compute wavelet coefficients and perform normalization
                [coefs,L]=wtAnalysis(ult_seg,5,'db4',350,0);
                cD4=coefs(3,1:L(3));
                abs_cD4=abs(cD4);
                if max(abs_cD4)>max_cD4
                    max_cD4=max(abs_cD4);
                end
%           Start acquisition
            elseif temp_adc
               estado2=0;
%              Analyze last second of the signal
               ult_seg=data_temp((length(data_temp(:,1))-fs):length(data_temp(:,1)),...
                    1);
%              Compute wavelet coefficients
               [coefs,L]=wtAnalysis(ult_seg,5,'db4',350,0);
               cD4=coefs(3,1:L(3));
%              Feature extraction for classifier input
               round_c4=round(length(cD4)*window_a);
               c4_temp=[];
               for j=0:(round_c4-1)
                   c4_temp=[cD4(length(cD4)-j),c4_temp];
               end
               cD4=c4_temp;
               abs_cD4=abs(cD4);
               mean_cD4=mean(abs_cD4);
               norm_cD4=max(abs_cD4)/max_cD4;
               if norm_cD4>1
                   norm_cD4=1;
               end
%              Evaluate ANFIS classifier               
               estado1=evalfis([norm_cD4,mean_cD4],fuzzy);
               if(estado1>2.6)
                   ojos=1;
                   estado1=3.3;
               else
                   ojos=0;
                   estado1=0;
               end
%              Print results obtained
               fprintf(1,'max: %.3f current: %.3f norm: %.2f ',...
                   max_cD4,max(abs(cD4)),norm_cD4);
               fprintf(1,'mean: %.2f eyes: %u time: %s\r',...
                   mean_cD4,ojos,num2str(length(data_temp(:,1))/fs,3));
               output=[output;max_cD4,max(cD4),norm_cD4,ojos,...
                   length(data_temp(:,1))/fs,mean_cD4];
            end
%           Send information to microcontroller
            ni_o.outputSingleScan([estado1 estado2]);
%           Determine processing and classification time
            proc_temp=toc;
            total_pt=total_pt+proc_temp;
            times=times+1;
            if proc_temp>proc_t
                proc_t=proc_temp;
            end
%           Determine time vector for first plot
            if(length(data_temp(:,1))>np)
                temp_plot1=data_temp(length(data_temp)-np:length(data_temp),1);
            else
                temp_plot1=data_temp(:,1);
            end
%           Plot acquired data of the EEG channel
            plot(handles.chan1,temp_plot1,'r');
            title(handles.chan1,'EEG Signal: O1-P3',...
                'fontsize',12,'FontWeight','bold');
            xlim(handles.chan1,[1,np]);
            xlabel(handles.chan1,'Samples');
            ylabel(handles.chan1,'Voltage (mV)');
%           Determine time vector for second plot
            if(length(data_temp(:,3))>np) 
                temp_plot3=data_temp(length(data_temp)-np:length(data_temp),3);
            else
                temp_plot3=data_temp(:,3);
            end
%           Plot acquired data of the EOG channel
            plot(handles.chan3,temp_plot3,'b');
            xlim(handles.chan3,[1,np]);
            title(handles.chan3,'EOG Signal: Right Eye',...
                'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan3,'Samples');
            ylabel(handles.chan3,'Voltage (mV)');
%           Update plots
            drawnow; 
%           Get snapshot for the video in the second window
            if register
                axes(handles.webcam2);
                image(getsnapshot(vid));
                set(gca,'xlim',[300,900],'ylim',[600 1000]);
                axis off;
            end
        end
%       Set button to start; user can press it again.
        set(hObject,'Value',0)
%       Close sessions
        ni_o.outputSingleScan([0 0]);
        ni_o.stop;
        ni_o.delete;
        lh.delete;
        ni.stop;    
        ni.delete;
%       Create data vectors and calculate processing time
        data=data_plot;
        y1=data_plot(:,1);
        y3=data_plot(:,3);
        tx=(0:1/fs:(length(y1)-1)/fs)';
        total_pt=total_pt/times;
%       Print processing time information
        fprintf(1,'Average Processing Time: %.4f seconds.\r',...
           total_pt);
        fprintf(1,'Maximum Processing Time: %.4f seconds.\r\n',...
           proc_t);
        clear time;
%       Plot all acquired data from EEG channel
        p1=plot(handles.chan1,tx,y1,'r');
        title(handles.chan1,'EEG Signal: O1-P3','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan1,[0,max(tx)]);
        ylim(handles.chan1,[min(y1),max(y1)]);
        xlabel(handles.chan1,'Time (sec)');
        ylabel(handles.chan1,'Voltage (mV)');
%       Plot all acquired data from EOG channel
        p3=plot(handles.chan3,tx,y3,'b');
        title(handles.chan3,'EOG Signal: Right Eye','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan3,[0,max(tx)]);
        ylim(handles.chan3,[min(y3),max(y3)]);
        xlabel(handles.chan3,'Time (sec)');
        ylabel(handles.chan3,'Voltage (mV)');
        
        
%   The comments of the EEG and ECG monitoring code (case 2 and 3) are
%   equivalent to the comments of the BMI application.
    case 2 
        set(handles.chan1,'Position',[0.088 0.721 0.650 0.223]);
        set(handles.chan2,'Position',[0.088 0.398 0.650 0.223]);
        set(handles.chan3,'Position',[0.088 0.075 0.650 0.221]);                         %BMI
        while (ni.IsRunning&&get (hObject, 'Value'))
            ni_o.outputSingleScan([estado1 estado2]);
            data_temp=data;
            if(length(data_temp(:,1))>np) 
                temp_plot1=data_temp(length(data_temp)-np:length(data_temp),1);
            else
                temp_plot1=data_temp(:,1);
            end
            plot(handles.chan1,temp_plot1,'r');
            xlim(handles.chan1,[1,np]);        
            title(handles.chan1,'EEG Signal: O1-P3',...
            'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan1,'Samples');
            ylabel(handles.chan1,'Voltage (mV)');
            
             if(length(data_temp(:,2))>np) 
                temp_plot2=data_temp(length(data_temp)-np:length(data_temp),2);
            else
                temp_plot2=data_temp(:,2);
            end
            plot(handles.chan2,temp_plot2,'b');
            xlim(handles.chan2,[1,np]);
            title(handles.chan2,'EOG Signal: Left Eye',...
                'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan2,'Samples');
            ylabel(handles.chan2,'Voltage (mV)');
            
            if(length(data_temp(:,3))>np) 
                temp_plot3=data_temp(length(data_temp)-np:length(data_temp),3);
            else
                temp_plot3=data_temp(:,3);
            end
            plot(handles.chan3,temp_plot3,'k');
            xlim(handles.chan3,[1,np]);
            title(handles.chan3,'EOG Signal: Right Eye',...
                'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan3,'Samples');
            ylabel(handles.chan3,'Voltage (mV)');
            
            drawnow; % forces MATLAB to update the figure
            axes(handles.webcam2);
            image(getsnapshot(vid));
            set(gca,'xlim',[300,900],'ylim',[600 1000]);
            axis off;
        end
        set(hObject,'Value',0)
        ni_o.outputSingleScan([0 0]);
        ni_o.stop;
        ni_o.delete;
        lh.delete;
        % disp(ni.ScansAcquired);
        ni.stop;
        ni.delete;
        y1=data(:,1);
        y2=data(:,2);
        y3=data(:,3);
        tx=(0:1/fs:(length(y1)-1)/fs)';
        
        p1=plot(handles.chan1,tx,y1,'r');
        title(handles.chan1,'EEG Signal: O1-P3','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan1,[0,max(tx)]);
        ylim(handles.chan1,[min(y1),max(y1)]);
        xlabel(handles.chan1,'Time (sec)');
        ylabel(handles.chan1,'Voltage (mV)');

        p2=plot(handles.chan2,tx,y2,'b');
        title(handles.chan2,'EOG Signal: Left Eye','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan2,[0,max(tx)]);
        ylim(handles.chan2,[min(y2),max(y2)]);
        xlabel(handles.chan2,'Time (sec)');
        ylabel(handles.chan2,'Voltage (mV)');

        p3=plot(handles.chan3,tx,y3,'k');
        title(handles.chan3,'EOG Signal: Right Eye','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan3,[0,max(tx)]);
        ylim(handles.chan3,[min(y3),max(y3)]);
        xlabel(handles.chan3,'Time (sec)');
        ylabel(handles.chan3,'Voltage (mV)');
    case 3 %ECG (Bipolares)
        set(handles.chan1,'Position',[0.088 0.721 0.650 0.223]);
        set(handles.chan2,'Position',[0.088 0.398 0.650 0.223]);
        set(handles.chan3,'Position',[0.088 0.075 0.650 0.221]);
        while (ni.IsRunning&&get (hObject, 'Value'))
            ni_o.outputSingleScan([estado1 estado2]);
            data_temp=data;
            if(length(data_temp(:,1))>np) 
                temp_plot1=data_temp(length(data_temp)-np:length(data_temp),1);
            else
                temp_plot1=data_temp(:,1);
            end
            plot(handles.chan1,temp_plot1,'r');
            xlim(handles.chan1,[1,np]);
            title(handles.chan1,'ECG Signal: Lead I',...
            'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan1,'Samples');
            ylabel(handles.chan1,'Voltage (mV)');
            
            if(length(data_temp(:,2))>np) 
                temp_plot2=data_temp(length(data_temp)-np:length(data_temp),2);
            else
                temp_plot2=data_temp(:,2);
            end
            plot(handles.chan2,temp_plot2,'b');
            xlim(handles.chan2,[1,np]);
            title(handles.chan2,'ECG Signal: Lead II',...
                'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan2,'Samples');
            ylabel(handles.chan2,'Voltage (mV)');
            
            if(length(data_temp(:,3))>np) 
                temp_plot3=data_temp(length(data_temp)-np:length(data_temp),3);
            else
                temp_plot3=data_temp(:,3);
            end
            plot(handles.chan3,temp_plot3,'k');
            xlim(handles.chan3,[1,np]);
            title(handles.chan3,'ECG Signal: Lead III',...
                'fontsize',12,'FontWeight','bold');
            xlabel(handles.chan3,'Samples');
            ylabel(handles.chan3,'Voltage (mV)');
            
            drawnow; % forces MATLAB to update the figure
        end
        set(hObject,'Value',0)
        ni_o.outputSingleScan([0 0]);
        ni_o.stop;
        ni_o.delete;
        lh.delete;
        % disp(ni.ScansAcquired);
        ni.stop;
        ni.delete;
        y1=data(:,1);
        y2=data(:,2);
        y3=data(:,3);
        tx=(0:1/fs:(length(y1)-1)/fs)';
        
        p1=plot(handles.chan1,tx,y1,'r');
        title(handles.chan1,'ECG Signal: Lead I','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan1,[0,max(tx)]);
        ylim(handles.chan1,[min(y1),max(y1)]);
        xlabel(handles.chan1,'Time (sec)');
        ylabel(handles.chan1,'Voltage (mV)');

        p2=plot(handles.chan2,tx,y2,'b');
        title(handles.chan2,'ECG Signal: Lead II','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan2,[0,max(tx)]);
        ylim(handles.chan2,[min(y2),max(y2)]);
        xlabel(handles.chan2,'Time (sec)');
        ylabel(handles.chan2,'Voltage (mV)');

        p3=plot(handles.chan3,tx,y3,'k');
        title(handles.chan3,'ECG Signal: Lead III','fontsize',12,...
            'FontWeight','bold');
        xlim(handles.chan3,[0,max(tx)]);
        ylim(handles.chan3,[min(y3),max(y3)]);
        xlabel(handles.chan3,'Time (sec)');
        ylabel(handles.chan3,'Voltage (mV)');
end
set(handles.plot_button,'String','Start')
pause(0.5);
end


% --- Outputs from this function are returned to the command line.
function varargout = bmi_three_channels_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


function f1_input_Callback(hObject, eventdata, handles)
% Validate that the text in the f1 field converts to a real number
f1 = str2double(get(hObject,'String'));
if isnan(f1) || ~isreal(f1)  
    % isdouble returns NaN for non-numbers and f1 cannot be complex
    % Disable the Plot button and change its string to say why
    set(handles.plot_button,'String','Cannot plot f1')
    set(handles.plot_button,'Enable','off')
    % Give the edit text box focus so user can correct the error
    uicontrol(hObject)
else 
    % Enable the Plot button with its original name
    set(handles.plot_button,'String','Start/Stop')
    set(handles.plot_button,'Enable','on')
end


function f2_input_Callback(hObject, eventdata, handles)
% Validate that the text in the f2 field converts to a real number
f2 = str2double(get(hObject,'String'));
if isnan(f2) ...          % isdouble returns NaN for non-numbers
        || ~isreal(f2)    % f1 should not be complex
    % Disable the Plot button and change its string to say why
    set(handles.plot_button,'String','Cannot plot f2')
    set(handles.plot_button,'Enable','off')
    % Give the edit text box focus so user can correct the error
    uicontrol(hObject)
else 
    % Enable the Plot button with its original name
    set(handles.plot_button,'String','Start/Stop')
    set(handles.plot_button,'Enable','on')
end


function t_input_Callback(hObject, eventdata, handles)
% Disable the Plot button ... until proven innocent
set(handles.plot_button,'Enable','off')
t = str2double(get(hObject,'String'));
if isnan(t) ...          % isdouble returns NaN for non-numbers
        || ~isreal(t)    % f1 should not be complex
    % Disable the Plot button and change its string to say why
    set(handles.plot_button,'String','Cannot plot f2')
    set(handles.plot_button,'Enable','off')
    % Give the edit text box focus so user can correct the error
    uicontrol(hObject)
else 
    % Enable the Plot button with its original name
    set(handles.plot_button,'String','Start/Stop')
    set(handles.plot_button,'Enable','on')
end


% --- Executes on selection change in menu_registros.
function menu_registros_Callback(hObject, eventdata, handles)


% --- Executes during object creation, after setting all properties.
function menu_registros_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject, 'String', {'EEG', 'BMI (EEG-EOG)', 'ECG (Bipolar)'});


% --- Executes on button press in analizar.
function analizar_Callback(hObject, eventdata, handles)
% THIS FUNCTION ANALIZE THEACQUIRED DATA USING SHORT-TIME FOURIER 
% TRANSFORM (STFT)

% Declare global variables
global caso

% Do it only for EEG signals not ECG
if(caso==1)
    eeg_analyzer
elseif(caso==2)
    eeg_analyzer
end

% --- Executes on button press in pushbutton5.
function pushbutton5_Callback(hObject, eventdata, handles)
% THIS FUNCTION WILL SAVE THE ACQUIRED DATA AS .MAT FILE
global analyze
analyze=0;
save_data


% --------------------------------------------------------------------
function uipushtool3_ClickedCallback(hObject, eventdata, handles)
% THIS FUNCTION WILL SAVE THE GUI IMAGE AS .EPS IN COLOR AND BLACK & WHITE
global p2 caso
% For the BMI hide the channel 2 (it is not used)
if caso~=1
   set(p2, 'LineWidth', 1);
end
% Save the GUI as image
saveas(handles.chan1,'GUI_BW.eps');
axes(handles.chan1);    
shohid = get(0,'ShowHiddenHandles');
set(0,'ShowHiddenHandles','on')
saveas(gcf,'GUI_COLOR.eps','epsc')
set(0,'ShowHiddenHandles',shohid);

Contact us