main.m

EOG Correction program is used to work out the amount of ocular potential recorded by the electroenc

    epoch1_count=0; epoch2_count=epoch1_count; epoch3_count=epoch1_count; epoch4_count=epoch1_count; epoch5_count=epoch1_count; 

    for rowidx=1:length(type)
        if type(rowidx)==stimulus_code(1)
            epoch1_count=epoch1_count+1;
        elseif type(rowidx)==stimulus_code(2)
            epoch2_count=epoch2_count+1;
        elseif type(rowidx)==stimulus_code(3)
            epoch3_count=epoch3_count+1;
        elseif type(rowidx)==stimulus_code(4)
            epoch4_count=epoch4_count+1;
        elseif type(rowidx)==stimulus_code(5)
            epoch5_count=epoch5_count+1;
        end        
    end

    %% II. Title: Data Preprocessing
    %% 1. Epoch and Baseline Removal
    % EM data: -200ms to 600ms (type 1 - 4)
    %
    % Blink data: -200ms to 800ms (type 5)
    %
    % Epoch data types only from 1 to 5

    tl=200e-3;      % lower time range
    th=600e-3;      % upper time range for EM data
    tbl=800e-3;     % upper time range for blink data

    nl=tl/per;      % no of datapoints for lower range
    nh=th/per;      % no of datapoints for upper range for EM data
    nbl=tbl/per;    % no of datapoints for upper range for blink data
    no=datapnts+1;  % mean or the centre of epoch

    range=(nh+nl);          % range of epoch for EM data
    range_blink=(nbl+nl);   % range of epoch for blink data

    % Initialisation for the epoch data
    epoch1=zeros((range+2)*epoch1_count,length(ch_label));
    epoch2=zeros((range+2)*epoch2_count,length(ch_label));
    epoch3=zeros((range+2)*epoch3_count,length(ch_label));
    epoch4=zeros((range+2)*epoch4_count,length(ch_label));
    epoch5=zeros((range_blink+2)*epoch5_count,length(ch_label));

    % Initialisation of yss data
    yss1=zeros(epoch1_count,length(ch_label));
    yss2=zeros(epoch2_count,length(ch_label));
    yss3=zeros(epoch3_count,length(ch_label));
    yss4=zeros(epoch4_count,length(ch_label));
    yss5=zeros(epoch5_count,length(ch_label));
    
    % Maybe these data types can be changed to single instead. 
    rowstart1=1; rowstart2=rowstart1; rowstart3=rowstart1; rowstart4=rowstart1; rowstart5=rowstart1;   % init counts
    ssidx1=0; ssidx2=ssidx1; ssidx3=ssidx1; ssidx4=ssidx1; ssidx5=ssidx1;   % init counts
    
    keep('hMainGui','data','ch_label','type','range',...
        'range_blink','rowstart1','rowstart2','rowstart3','rowstart4','rowstart5',...
        'no','nl','nh','nbl','ssidx1','ssidx2','ssidx3','ssidx4','ssidx5',...
        'yss1','yss2','yss3','yss4','yss5','epoch1','epoch2','epoch3','epoch4','epoch5','in','handles',...
        'hObject','eventdata','stimulus_code');
    
    try
        for data_rowidx=1:length(type)
            workbar(data_rowidx/length(type),'Step 1 of 4...','Program Status');
            if type(data_rowidx)==stimulus_code(1)
                    rowend1=range+rowstart1;
                    epoch1(rowstart1:rowend1,:)=data(no(data_rowidx)-nl:no(data_rowidx)+nh,:);
                    ssidx1=ssidx1+1;
                    yss1(ssidx1,:)=mean(epoch1(((range+1)*(ssidx1-1)+ssidx1:(range+1)*(ssidx1-1)+ssidx1+(nl-1)),:));     % 2. baseline removal
                    for colidx=1:length(ch_label)
                        epoch1((range+1)*(ssidx1-1)+ssidx1:(range+1)*(ssidx1-1)+ssidx1+range,colidx)=epoch1((range+1)*(ssidx1-1)+ssidx1:(range+1)*(ssidx1-1)+ssidx1+range,colidx)-yss1(ssidx1,colidx);
                    end
                    epoch1(rowend1+1,:)=nan;
                    rowstart1=rowend1+2;
            elseif type(data_rowidx)==stimulus_code(2)
                    rowend2=range+rowstart2;
                    epoch2(rowstart2:rowend2,:)=data(no(data_rowidx)-nl:no(data_rowidx)+nh,:);
                    ssidx2=ssidx2+1;
                    yss2(ssidx2,:)=mean(epoch2(((range+1)*(ssidx2-1)+ssidx2:(range+1)*(ssidx2-1)+ssidx2+(nl-1)),:));
                    for colidx=1:length(ch_label)
                        epoch2((range+1)*(ssidx2-1)+ssidx2:(range+1)*(ssidx2-1)+ssidx2+range,colidx)=epoch2((range+1)*(ssidx2-1)+ssidx2:(range+1)*(ssidx2-1)+ssidx2+range,colidx)-yss2(ssidx2,colidx);
                    end
                    epoch2(rowend2+1,:)=nan;
                    rowstart2=rowend2+2;
            elseif type(data_rowidx)==stimulus_code(3)
                    rowend3=range+rowstart3;
                    epoch3(rowstart3:rowend3,:)=data(no(data_rowidx)-nl:no(data_rowidx)+nh,:);
                    ssidx3=ssidx3+1;
                    yss3(ssidx3,:)=mean(epoch3(((range+1)*(ssidx3-1)+ssidx3:(range+1)*(ssidx3-1)+ssidx3+(nl-1)),:));
                    for colidx=1:length(ch_label)
                        epoch3((range+1)*(ssidx3-1)+ssidx3:(range+1)*(ssidx3-1)+ssidx3+range,colidx)=epoch3((range+1)*(ssidx3-1)+ssidx3:(range+1)*(ssidx3-1)+ssidx3+range,colidx)-yss3(ssidx3,colidx);
                    end
                    epoch3(rowend3+1,:)=nan;
                    rowstart3=rowend3+2;
            elseif type(data_rowidx)==stimulus_code(4)
                    rowend4=range+rowstart4;
                    epoch4(rowstart4:rowend4,:)=data(no(data_rowidx)-nl:no(data_rowidx)+nh,:);
                    ssidx4=ssidx4+1;
                    yss4(ssidx4,:)=mean(epoch4(((range+1)*(ssidx4-1)+ssidx4:(range+1)*(ssidx4-1)+ssidx4+(nl-1)),:));
                    for colidx=1:length(ch_label)
                        epoch4((range+1)*(ssidx4-1)+ssidx4:(range+1)*(ssidx4-1)+ssidx4+range,colidx)=epoch4((range+1)*(ssidx4-1)+ssidx4:(range+1)*(ssidx4-1)+ssidx4+range,colidx)-yss4(ssidx4,colidx);
                    end
                    epoch4(rowend4+1,:)=nan;
                    rowstart4=rowend4+2;
            elseif type(data_rowidx)==stimulus_code(5)
                    rowend5=range_blink+rowstart5;
                    epoch5(rowstart5:rowend5,:)=data(no(data_rowidx)-nl:no(data_rowidx)+nbl,:);
                    ssidx5=ssidx5+1;
                    yss5(ssidx5,:)=mean(epoch5(((range_blink+1)*(ssidx5-1)+ssidx5:(range_blink+1)*(ssidx5-1)+ssidx5+(nl-1)),:));
                    for colidx=1:length(ch_label)
                        epoch5((range_blink+1)*(ssidx5-1)+ssidx5:(range_blink+1)*(ssidx5-1)+ssidx5+range_blink,colidx)=epoch5((range_blink+1)*(ssidx5-1)+ssidx5:(range_blink+1)*(ssidx5-1)+ssidx5+range_blink,colidx)-yss5(ssidx5,colidx);
                    end
                    epoch5(rowend5+1,:)=nan;
                    rowstart5=rowend5+2;
            end
        end
    catch
    end
    
    % Exit the workbar if data_rowidx<length(type)
    
    if data_rowidx<length(type)
        data_rowidx=length(type);
        workbar(data_rowidx/length(type),'Step 1 of 4...','Program Status');
    end

    % 2. Averaging of epochs data

    % Initialisation of em data (epoch mean)
    em1=yss1; em2=yss2; em3=yss3; em4=yss4; em5=yss5;
    
    keep('hMainGui','ch_label','range','range_blink',...
        'ssidx1','ssidx2','ssidx3','ssidx4','ssidx5',...
        'epoch1','epoch2','epoch3','epoch4','epoch5','in','handles',...
        'em1','em2','em3','em4','em5','hObject','eventdata');
    
    % Initialisation of avg data
    avg1=zeros(range+1,length(ch_label));
    avg2=avg1;
    avg3=avg1;
    avg4=avg1;
    avg5=zeros(range_blink+1,length(ch_label));
    avg_count=0;

    for data_rowidx=1:range+1
        for avg_rowidx=1:ssidx1
            em1(avg_rowidx,:)=epoch1((avg_rowidx-1)*(range+1)+avg_rowidx+(data_rowidx-1),1:length(ch_label));
            if avg_rowidx==ssidx1
                avg1(data_rowidx,:)=mean(em1,1);
                avg_count(1)=1;
            end
        end
        for avg_rowidx=1:ssidx2
            em2(avg_rowidx,:)=epoch2((avg_rowidx-1)*(range+1)+avg_rowidx+(data_rowidx-1),1:length(ch_label));
            if avg_rowidx==ssidx2
                avg2(data_rowidx,:)=mean(em2,1);
                avg_count(2)=2;
            end
        end
        for avg_rowidx=1:ssidx3
            em3(avg_rowidx,:)=epoch3((avg_rowidx-1)*(range+1)+avg_rowidx+(data_rowidx-1),1:length(ch_label));
            if avg_rowidx==ssidx3
                avg3(data_rowidx,:)=mean(em3,1);
                avg_count(3)=3;
            end
        end
        for avg_rowidx=1:ssidx4
            em4(avg_rowidx,:)=epoch4((avg_rowidx-1)*(range+1)+avg_rowidx+(data_rowidx-1),1:length(ch_label));
            if avg_rowidx==ssidx4
                avg4(data_rowidx,:)=mean(em4,1);
                avg_count(4)=4;
            end
        end
    end
    for data_rowidx=1:range_blink+1
        for avg_rowidx=1:ssidx5
            em5(avg_rowidx,:)=epoch5((avg_rowidx-1)*(range_blink+1)+avg_rowidx+(data_rowidx-1),1:length(ch_label));
            if avg_rowidx==ssidx5
                avg5(data_rowidx,:)=mean(em5,1);
                avg_count(5)=5;
            end
        end
    end
      
    % 3. Appending average data

    keep('hMainGui','ch_label','avg1','avg2','avg3','avg4','avg5','avg_count','in','handles','hObject','eventdata');
    
    avg=0;
       
    for avg_rowidx=1:length(avg_count)
        switch(avg_count(avg_rowidx))
            case 1
                avg=avg1;
            case 2
                avg=[avg; avg2];
            case 3
                avg=[avg; avg3];
            case 4
                avg=[avg; avg4];    % just up to 4 for the EM data
            case 5
                blink=avg5;         % blink data
        end
    end

    % avg is the EM (Eye Movement data)

    %% III. Title: Correction Factors Computation
    %% 1. Creating channel VE, HE, and RE for EM data
    
    keep('hMainGui','ch_label','avg','blink','avg_count','in','handles','hObject','eventdata');
       
    data_lim=5*2^20;    % limit data to 5MB
    
    [avg_size,colsize]=size(avg);  % colsize = length(ch_label)
        
    data_byte=avg_size*colsize*8;     % 8 for double type of data
       
    if data_byte>data_lim
        num_blocks=ceil(data_byte/data_lim);
        row_size=ceil(avg_size/num_blocks);
    else
        num_blocks=1;
        row_size=avg_size;
    end

    row_start=1;
    for block=1:num_blocks
        row_end=(row_start-1)+row_size;
        if row_end>avg_size
            row_end=avg_size;
        end
        for i=1:length(ch_label)
            ch.(ch_label{i})=avg(row_start:row_end,i);
            y=[ch_label{i},'(row_start:row_end,1)=ch.',ch_label{i},';'];
            eval(y);
        end
        if block~=num_blocks
            row_start=row_end+1;
        end
    end
    
    VE=zeros(avg_size,1); HE=VE; RE=VE; % Pre-allocation with zeros
    
    VE=eval(in{1}); % Assigning "in" variable to VE, HE, and RE
    HE=eval(in{2});
    RE=eval(in{3});

    %% 2.  Calculating Bv, Bh from EM data
    
    EM=zeros(avg_size,colsize+1);
    EM=[avg RE];
    
    keep('hMainGui','ch_label','avg_size','colsize','EM','blink','avg_count','in','handles','VE','HE','data_lim','hObject','eventdata');
    
    v1=zeros(avg_size,3); coeff1=zeros(3,colsize+1); Bv=zeros(1,colsize+1); Bh=zeros(1,colsize+1);  % Pre-allocation with zeros
    
    v1=[ones(size(VE)) VE HE];  % creating vandermode matrix
    coeff1=v1\EM;
    Bv=coeff1(2,:);
    Bh=coeff1(3,:);

    %% 3. Calculating Br from blink data
    % Creating VE, HE, and RE for blink data.
    %
    % 1st half of correction
    %
    % $$corrdata=data-\beta_v*VE-\beta_h*HE$$

    keep('hMainGui','ch_label','blink','avg_count','in','handles','Bv','Bh','data_lim','hObject','eventdata');

    [blink_size,colsize]=size(blink);  % colsize = length(ch_label)

    data_byte=blink_size*colsize*8;     % 8 for double type of data

    if data_byte>data_lim
        num_blocks=ceil(data_byte/data_lim);
        row_size=ceil(blink_size/num_blocks);
    else
        num_blocks=1;
        row_size=blink_size;
    end

    row_start=1;
    for block=1:num_blocks
        row_end=(row_start-1)+row_size;
        if row_end>blink_size
            row_end=blink_size;
        end
        for i=1:length(ch_label)
            ch.(ch_label{i})=blink(row_start:row_end,i);
            y=[ch_label{i},'(row_start:row_end,1)=ch.',ch_label{i},';'];
            eval(y);
        end
        if block~=num_blocks
            row_start=row_end+1;
        end
    end

    VE=zeros(blink_size,1); HE=VE; RE=VE; % Pre-allocation with zeros

    VE=eval(in{1});
    HE=eval(in{2});
    RE=eval(in{3});

    keep('hMainGui','handles','avg_count','ch_label','blink','in','Bv','Bh','VE','HE','RE',...
        'num_blocks','row_size','data_lim','colsize','blink_size','hObject','eventdata');

    blink=[blink RE];

    corr1=zeros(blink_size,colsize);

    row_start=1;
    for block=1:num_blocks
        row_end=(row_start-1)+row_size;
        if row_end>blink_size
            row_end=blink_size;
        end
        for colidx=1:colsize+1
            if colidx~=colsize+1
                corr1(row_start:row_end,colidx)=blink(row_start:row_end,colidx)-VE(row_start:row_end)*Bv(colidx)-HE(row_start:row_end)*Bh(colidx); % 1st half of correction
            else
                RE(row_start:row_end,1)=blink(row_start:row_end,colidx)-VE(row_start:row_end)*Bv(colidx)-HE(row_start:row_end)*Bh(colidx);  % RE channel
            end
        end
        if block~=num_blocks
            row_start=row_end+1;
        end
    end

    keep('hMainGui','blink_size','colsize','handles','avg_count','ch_label','corr1','RE','Bv','Bh','in','hObject','eventdata');

    v2=zeros(blink_size,2); coeff2=zeros(2,colsize); Br=zeros(1,colsize);  % Pre-allocation with zeros

    v2=[ones(size(RE)) RE];
    coeff2=v2\corr1;
    Br=coeff2(2,:);

    %% IV. Title: Correction for To Be Corrected (TBC) data
    %% 1. Creating channel VE, HE, and RE for the tbc_data
    keep('hMainGui','handles','avg_count','ch_label','Bv','Bh','Br','in','hObject','eventdata');
    
    tbc_data=getappdata(hMainGui,'tbc_data');
    
    [tbc_size,colsize]=size(tbc_data);  % colsize = length(ch_label)