image thumbnail

eogui – a software to analyze electro-oculogram (EOG) recordings

by

 

10 Aug 2011 (Updated )

detects blinks and saccadic eye movements in EOG recordings

analyzeeog.m
function [sakkaden, blinks, varargout]=analyzeeog(input,parameter,varargin)
% erkennt die Ereignise innerhalb des Signals
% Parameter:
%       input = n*3 Matrix, n = Anzahl der Samples, diskretes Signal
%       parameter = Struktur mit den Parameter
%
% Ausgabe:
%       [sakkade, blinks]
%           sakkade = Liste aller erkannten Sakkaden
%           blinks = Liste aller erkannten Blinks
%
% $Id: analyzeeog.m 15 2005-03-10 12:53:40Z maik $


%------------------------------------------------------------------
% Problem 1:
% Die Verrechnung mit dem Eingangsoffset stimmt noch nicht.
% realpos=pos+toffset
% Falls der Eingangsoffset jedoch in der richtigen Freq. ist, passt es
%
% pos in samplenr. <-!!!!!  muss noch in msec umgerechnet werden
% realpos in msec
% toffset in msec
%------------------------------------------------------------------

mode=1;
if(nargin==3)
    mode=varargin{1};
end


sakkaden=[];
blinks=[];
algodata=[];
grobdebug=[];
artefakte=[];


%
%  Samples
%  |
%  |  +-----------+
%  +->|Filterstufe|-+---------------------------------+
%     |    1      | |                                 |
%     +-----------+ |                                 |
%                   |  +------------+    +---------+  |
%                   +->| Filterstufe|--->|Grobsuche|  |
%                      |     2      |    |         |  |
%                      +------------+    +---------+  |
%                                              |      |
%                                              \/     \/
%    /-----------\      +------------+     +-----------+
%    |  Ist      |<-----| Parameter- |<----| Feinsuche |
%    |  Blink    |      | bestimmung |     |           |
%    |  moeglich |      +------------+     +-----------+
%    |  ?        |
%    \-----------/
%          |   |         NEIN
%       JA |   +-----------------------------+
%          |                                 |
%          \/                                \/
%    +---------------------------+    +---------------------------------+
%    | 1. Auf-,Absakk. bestimmen |    | 1. Bestimme restliche Attribute |
%    | 2. Rest. Attr. bestimmen  |    | 2. speichere Sakkaden           |
%    | 3. Blink speichern        |    |                                 |
%    +---------------------------+    +---------------------------------+



%prfe input
%-------------------------------------
if( length(input(1,:)) <3)
    disp('Fehler, es werden mindestens 3 Kanle bentigt (Zeit,X,Y)');
    return;
end




%prfe parameter
%-------------------------------------
%parameter.f1(n).para  % parameter fr den eingesetzten Filter
%parameter.f1(n).type  % Art des eingesetzten Filter (z.Z. mean,median,notch)

%parameter.f2(n).para  % parameter fr den eingesetzten Filter
%parameter.f2(n).type  % Art des eingesetzten Filter (z.Z. mean,median,notch)
%parameter.grobsuche.plateaufenster


% berechne gefilterte Signale und passe die verschiebung (Einschwingphasen) an
% sig1 = leicht gesmooth
% sig2 = stark geglttet zur Grobsuche
channel1.sig0=input(:,2)*parameter.wfakt.x;
[channel1.sig1 offset]=feval(parameter.filteralgo,input(:,2)*parameter.wfakt.x,parameter.filter1);


channel1.sig0=channel1.sig0(1+offset:end);

channel1.offset=offset;
[channel1.sig2 offset]=feval(parameter.filteralgo,channel1.sig1,parameter.filter2);
channel1.offset=channel1.offset+offset;
channel1.sig0=channel1.sig0(1+offset:end);
channel1.sig1=channel1.sig1(1+offset:end);
channel1.sig3=diff(channel1.sig1)*parameter.sfreq; % Speed in Grad /sec

channel1.sig0=channel1.sig0(2:end);
channel1.sig1=channel1.sig1(2:end); % schneide ersten sample ab, wegen diff
channel1.sig2=channel1.sig2(2:end); % schneide ersten sample ab, wegen diff
channel1.offset=channel1.offset+1;


channel2.sig0=input(:,3)*parameter.wfakt.y;
[channel2.sig1 offset]=feval(parameter.filteralgo,input(:,3)*parameter.wfakt.y,parameter.filter1);
channel2.sig0=channel2.sig0(1+offset:end);
channel2.offset=offset;
[channel2.sig2 offset]=feval(parameter.filteralgo,channel2.sig1,parameter.filter2);
channel2.offset=channel2.offset+offset;
channel2.sig0=channel2.sig0(1+offset:end);
channel2.sig1=channel2.sig1(1+offset:end);
channel2.sig3=diff(channel2.sig1)*parameter.sfreq; % Speed in Grad /sec
channel2.sig4=diff(channel2.sig2)*parameter.sfreq; % Speed in Grad /sec

channel2.sig0=channel2.sig0(2:end);
channel2.sig1=channel2.sig1(2:end); % schneide ersten sample ab, wegen diff
channel2.sig2=channel2.sig2(2:end); % schneide ersten sample ab, wegen diff
channel2.offset=channel2.offset+1;


if(mode==1) disp('Preprocessing finished...'); end
offset=channel2.offset; %zustzlicher offset fr Filterverschiebungen.
zeitachse=double(input(1+offset:end,1)); % fr die Darstellung




%------------------------------------------------------------
%suche Blinks
rohBlinks=[];
tblinks = findBlinks(channel2.sig4,parameter);
tcount=size(tblinks);
for i=1:tcount(1)

    wfakt=parameter.wfakt.y;
    blink.bstart1=zeitachse(tblinks(i,1));
    blink.bende1=zeitachse(tblinks(i,2));
    blink.bstartamp1 = channel2.sig2(tblinks(i,1));

%     if (tblinks(i,1)-(parameter.grobsuche.plateaufenster/2)<1)
%         blink.bstartamp1 = channel2.sig1(tblinks(i,1));
%     else
%         blink.bstartamp1 = channel2.sig1(tblinks(i,1)-round(parameter.grobsuche.plateaufenster/4));
%     end


    blink.bendeamp1 = channel2.sig2(tblinks(i,2));
    blink.bamp1 =(blink.bendeamp1-blink.bstartamp1);

    blink.bstart2=zeitachse(tblinks(i,3));
    blink.bende2=zeitachse(tblinks(i,4));
    blink.bstartamp2 =channel2.sig2(tblinks(i,3));

%     if(tblinks(i,4)+(parameter.grobsuche.plateaufenster/2) > length(channel2.sig1))
%         blink.bendeamp2 =channel2.sig1(tblinks(i,4));
%     else
%         blink.bendeamp2 =channel2.sig1(tblinks(i,4)+round(parameter.grobsuche.plateaufenster/4));
%     end
    blink.bendeamp2 =channel2.sig2(tblinks(i,4));
    blink.bamp2 =(blink.bendeamp2-blink.bstartamp2);


    blink.qmaxspeedpos1=zeitachse(tblinks(i,5));
    blink.qmaxspeedpos2=zeitachse(tblinks(i,6));
    blink.qmaxspeed1=channel2.sig3(tblinks(i,5));
    blink.qmaxspeed2=channel2.sig3(tblinks(i,6));


    % rechne speedpos in Prozent um ...
    blink.qmaxspeedpos1_proz=...
    (blink.qmaxspeedpos1-blink.bstart1)/...
        (blink.bende1-blink.bstart1)*100;

    blink.qmaxspeedpos2_proz=...
    (blink.qmaxspeedpos2-blink.bstart2)/...
        (blink.bende2-blink.bstart2)*100;

    blink.flag='true';

    % prfe Plausi
    if(blink.bamp1 <parameter.grobsuche.y.schwellwert)
        continue;
    end

    if(abs(blink.bamp2) <parameter.grobsuche.y.schwellwert)
        continue;
    end

    if(blink.bende1-blink.bstart1 <parameter.plausi.min_dauer)
        continue;
    end

    if(blink.bende2-blink.bstart2 <parameter.plausi.min_dauer)
        continue;
    end



    if(blink.bamp1>=parameter.blinkamp) &&...
        (abs(blink.bamp2)>=parameter.blinkamp)
        rohBlinks=[rohBlinks blink]; % speichere blink
    end


end
%------------------------------------------------------------



% Untersuche komplettes Signal, Hauptloop
samplecount=length(channel2.sig2);
position=1;
laststat=-inf;  %letzte Statusaktuallisierung -inf (wegen buttonbug)
refresh=samplecount*3/100;          %alle ca 5 % aktuallisieren

%%hhh = waitbar(0,'Fortschritt: 0%');
if(mode==1)
    hhh = waitbar(0,'Progress: 0%','CreateCancelBtn','set(gcf,''UserData'',''cancel'' );');
end

while(position<samplecount-parameter.grobsuche.plateaufenster) % durchlaufe alle Samples

    % Statusausgabe alle n Samples !

    if(mode==1 && laststat+refresh<=position)
        name=get(hhh,'UserData');
        if(~ishandle(hhh) || strcmp(name,'cancel'))
            warning('Preprocessing cancelled');
            if(ishandle(hhh)) delete(hhh), end;
            if(nargout==3)
                varargout{1}=[];
            end
            sakkaden=[];
            blinks=[];
            return;
        end
        waitbar(position/samplecount,hhh,sprintf('Progress: %d%%',round(100*position/samplecount)));
        %[position samplecount]
        laststat=position;
    end

    blink_found=false; % ist im aktuellen Segment ein Blink ?

    % wurde in einem Kanal eine Schwellwertberschreitung gefunden ?
     lastblinkpos=1;
     if ((abs(channel1.sig2(position) - channel1.sig2(position+parameter.grobsuche.plateaufenster)) >parameter.grobsuche.x.schwellwert) || ...
             ((abs(channel2.sig2(position) - channel2.sig2(position+parameter.grobsuche.plateaufenster)) >parameter.grobsuche.y.schwellwert)))

                %-------------------------------------------------------------
                % Grobsuche, bestimmt lnge und mglichen Typ
                %-------------------------------------------------------------
                [len xflg yflg artefakt blinkflg] = grobsuche(channel1.sig2,channel2.sig2,channel1.sig3,channel2.sig3,...
                    position,parameter.grobsuche);

                % wurde das Dateiende erkannt, sofort abbrechen
                if(artefakt{1}==4)
                    disp(artefakt{2});
                    break;
                end

                %prfe ob Blink im Bereich ist
                %-------------------------------------
                for i=lastblinkpos:length(rohBlinks)
                   blink=rohBlinks(i);

                   if(blink.bstart1>zeitachse(position+len)) break; end; % Blink kommt erst  nach dem Bereich -> ende

                   % links berlappt
                   if(blink.bstart1>zeitachse(position) && blink.bstart1<zeitachse(position+len))
                       yflg=false;
                       blink_found=true;
                       if(artefakt{1}~=0)
                           rohBlinks(i).flag='invalide';
                       end
                       lastblinkpos=i;
                   end

                   % rechts berlappt
                   if(blink.bende2>zeitachse(position) && blink.bende2<zeitachse(position+len))
                       yflg=false;
                       blink_found=true;
                       if(artefakt{1}~=0)
                           rohBlinks(i).flag='invalide';
                       end
                       lastblinkpos=i;
                   end

                   % innerhalb
                   if(blink.bstart1<=zeitachse(position) && blink.bende2>=zeitachse(position+len))
                       yflg=false;
                       blink_found=true;
                       if(artefakt{1}~=0)
                           rohBlinks(i).flag='invalide';
                       end
                       lastblinkpos=i;
                   end


                end


                % wurde ein Artefakt erkannt, keine Auswertung
                if(artefakt{1}~=0)
                    artefakte=[artefakte; [zeitachse(position) zeitachse(position+len)]];
                    position=position+len;
                    continue;
                end

                %Auswertung
                sakk.start.amp.x.var =channel1.sig2(position);
                sakk.start.amp.y.var =channel2.sig2(position);

                sakk.ende.amp.x.var = channel1.sig2(position+len);   %merke die Endamplitude X Kanal des
		        sakk.ende.amp.y.var = channel2.sig2(position+len);   %merke die Endamplitude X Kanal des

                subList=[];
                if( (xflg==1) && (yflg==1))
                    %---------------------------------------------------------------
                    %  wenn in X und Y Kanal eine Sakkade gefunden wurde,
                    %  wird die Sakkade mit der groessten Amplitude gewaehlt.
                    %---------------------------------------------------------------
                    ampx =abs(sakk.start.amp.x.var - sakk.ende.amp.x.var);
                    ampy =abs(sakk.start.amp.y.var - sakk.ende.amp.y.var);
                    if ampx > ampy
                        yflg=0;
                    else
                        xflg=0;
                    end
                end  % xflg AND yflg


                %-------------------------------------------------------------------
                %  werte X Kanal aus
                %-------------------------------------------------------------------

                if(xflg ==1)
                    subList=feinauswertung(channel1.sig1,channel1.sig3,position,len,parameter.feinsuche.x.schwellwert);
                    if(~isempty(subList))
                        subList=checkplausi(subList,channel1.sig1,channel1.sig2,parameter.plausi,parameter.grobsuche.x.schwellwert);
                    end

                end

                %-------------------------------------------------------------------
                %  werte Y Kanal aus
                %-------------------------------------------------------------------
                if(yflg ==1)
                    subList=feinauswertung(channel2.sig1,channel2.sig3,position,len,parameter.feinsuche.y.schwellwert);
                    if(~isempty(subList))
                        subList=checkplausi(subList,channel2.sig1,channel2.sig2,parameter.plausi,parameter.grobsuche.x.schwellwert);
                    end
                end


                %-------------------------------------------------------------------
                %  speichere Sakkaden
                %-------------------------------------------------------------------
               if (~isempty(subList))


                        %----------------------------------------
                        %   berechen Amplituden
                        %----------------------------------------
%                         subList(1).startampx=sakk.start.amp.x.var;
%                         subList(1).startampy=sakk.start.amp.y.var;
%                         subList(length(subList)).endampx=sakk.ende.amp.x.var;
%                         subList(length(subList)).endampy=sakk.ende.amp.y.var;
%
%                         for i=1:length(subList)-1
%                             subList(i).endampx=channel1.sig1(subList(i).ende);
%                             subList(i).endampy=channel2.sig1(subList(i).ende);
%                             subList(i+1).startampx=channel1.sig1(subList(i+1).start);
%                             subList(i+1).startampy=channel2.sig1(subList(i+1).start);
%                         end

                        for i=1:length(subList)
                            subList(i).endampx=channel1.sig1(subList(i).ende);
                            subList(i).endampy=channel2.sig1(subList(i).ende);
                            subList(i).startampx=channel1.sig1(subList(i).start);
                            subList(i).startampy=channel2.sig1(subList(i).start);
                        end

                        %----------------------------------------
                        %   speicher Sakkaden
                        %----------------------------------------
                        for i=1:length(subList)
                            sakkade.channel=(xflg~=1);

                            sakkade.start= zeitachse(subList(i).start);
                            sakkade.ende= zeitachse(subList(i).ende);
                            sakkade.startampx =subList(i).startampx;
                            sakkade.endeampx =subList(i).endampx;
                            sakkade.ampx=(sakkade.endeampx-sakkade.startampx);

                            if(blink_found)
                                sakkade.startampy =0;
                                sakkade.endeampy =0;
                                sakkade.ampy=0;
                            else
                                sakkade.startampy =subList(i).startampy;
                                sakkade.endeampy =subList(i).endampy;
                                sakkade.ampy=(sakkade.endeampy-sakkade.startampy);
                            end


                            sakkade.qmaxspeedpos = zeitachse(subList(i).maxspeedpos);

                            % rechne speedpos Prozent um ...
                            sakkade.qmaxspeedpos_proz=...
                            (sakkade.qmaxspeedpos-sakkade.start)/...
                                (sakkade.ende-sakkade.start)*100;

                            % Geschwindigkeit umrechnen
                            sakkade.qmaxspeed = subList(i).maxspeed;


                            % Sakkade mu mindestens 2 Grad berbrcken
                            %if(sakkade.ampx>2 ||sakkade.ampy>2)
                                sakkaden=[sakkaden sakkade]; % speichere sakkaden
                            %end

                        end
               end


         % speichert die grobbereiche fr die mgliche Debugausgabe
         grobdebug=[grobdebug;[zeitachse(position) zeitachse(position+len)]];
         position=position+len;
     else % Schwellwert wurde nicht berschritten, ein sample weiter
         position=position+1;
     end
end % Hauptschleife


% Filtere ungltige Blinks aus
if ~isempty(rohBlinks)
    blinks=rohBlinks(find(~strcmp({rohBlinks(:).flag},'invalide')));
end

% for i=1:length(rohBlinks)
%    if(~strcmp(rohBlinks(i).flag,'invalide'))
%        blinks=[blinks rohBlinks(i)];
%    end
% end

if mode==1
    delete(hhh);
    disp('Final processing finished...');
end

if(nargout==3)
    algodata.channel1=channel1;
    algodata.channel2=channel2;
    algodata.zeitachse=zeitachse;
    algodata.artefakte=artefakte;
    algodata.grobdebug=grobdebug;
    varargout{1}=algodata;
end

Contact us