%%% AVIV CHANGE BACK ALL OCCASIONS OF LFP_017 ETC. TO RAW_017


% updates from 100115
% Given the immensity of the number of nights we still have to take care
% of, we decided to move to a different strategy - using a high to low
% frequency ratio, combined with EMG RMS to decide which stage is it.
% So, the following changes were done. The normal staging is done for some
% 30 nights to ensure we eye-balled the data enough, and a parallel
% automatic staging is done for the same segments. Then, it will be
% compared and we will see if the automatic method is good enough.
% Practically:
% 1. Segments are now 10 seconds to precent pollution.
% 2. For every new segment, the calculation of the EMG RMS and high to low
%    ratio, as well as the number of EOG "events" (number of events
%    that exceed n times the std, are recorded.
% 3. We will scatter all the dots and decide the boundaries, and then do it
% automatically.
% 4. Another possibility is to so this by hand for the first part of the
% nigt and then automatically to the later part.

%TODO THE ISSUE OF SAMPLE BORDERS!
% when I have my system:
% 1. line: eval(['clear ',name_EEG_FR(1:end-8),' ',name_EEG_FL(1:end-8),'
% 2. line: load(polysom_analog(1,ii).name,'CLFP*') - according to the actual vars.
% 3. same for line eval(['polysom_analog(1,ii).EEG_FR=',name_EEG_FR,'-',name_Ground_1,';']); and same lines following.
% 4. same for line: name_EEG_FR='CLFP_017.Samples';.. and following.
% 5. line fs=CLFP_017.KHz*1000; % sampling frequency in Hz. Don't forget to fill the right one.

%% Important Remarks:
% 0. This code does the following: 
%    a. Receives all the raw .mat files recorded in a daily session.
%    b. extracts 16 analog traces (details below) from the raw data.
%    c. Substracts the ground from the EEG/EOG data and EMG_4 from all the 
%       rest to receive 13 referenced traces.
%    d. For every such referenced trace, concatenates traces from all files
%       recorded that day to receive 13 final, long, referenced traces.
%    e. Presents to the user, through GUI, the traces in 30 seconds
%       periods, and recives a staging result.
%    f. Creates an excel list of the 30-sec periods, with specific bagin-
%       and end-sample for every period. The begin- and end-samples are the
%       samples in the specific file which was used, so the table contains
%       also the relevant number of file.
%    g. Feeds the relevant data into the abridged mat file, to be used for
%       analysis, created by arrange_digital_old_data, along with stager 
%       name and date.
% % 1. This code is meant to be used just after the daily recording. It will
% retrieve the relevant data from the raw .mat files so the different
% traces (e.g. EEG_FR, EEG_FL etc.) used for staging are not saved anywhere as such.
% 2. Currently, the code is adapted to the version (Shiran's set up) where
% there are no gaps between one file and the next (i.e. if file 12 ends in
% sample 23000 then file 13 starts in 23001. This is the case in Shiran's
% dataset, and for this setting this code is adapted. In Eitan's dataset,
% "leaps" are possible between files (i.e 12 ends in 23000 and 13 starts in
% 23010 or even in 22998). The lines where this could be changed are _____
% and if I see that I'm working with Eitan, I sould be correcting it.

% The traces are recorded in the following order:
% 1.  EEG_FR (Frontal Right)
% 2.  EEG_FL (Frontal Left)
% 3.  EEG_CR (Central Right)
% 4.  EEG_CL (Central Left) 
% 5.  EEG_OR (Occipital Right) 
% 6.  EEG_OL (Occipital Left)
% 7.  EOG_UR (Upper Right (eye))
% 8.  EOG_UL (Upper Left (eye))
% 9.  EOG_LR (Lower Right (eye))
% 10. EOG_LL (Lower Left (eye))
% 11. EMG_1
% 12. EMG_2
% 13. EMG_3
% 14. EMG_4 (different EMG electrodes)
% (15.-16. GROUND - to substract from relevant files.

% The files used for the recording are all files recorded that day,
% regardless of whether I had units. That is for me to see and report the
% general architecture of the monkey's sleep.

% Because I wrote this program before I had any data, I used LFP data from
% SNr recordings by Mati. I changed the vars names to match the EEG_FR etc.
% system. The code starts in the starting point where the following
% variables are in the workspace (all from "polysom_ddmmyy"): 
% - BeginFile: A string with the file name of the first file.
% - EndFile  : A string with the file name of the last file.
% - RecBegin : The number of the first file.
% - RecEnd   : The number of the last file.
% - polysom_analog: A 1X5 structure, every cell having 14 analog variables:
%   EEG_FR, EEG_FL, EEG_CR, EEG_CL, EEG_OR, EEG_OL,
%   EOG_UR, EEG_UL, EEG_LR, EEG_LL,
%   EMG_1, EMG_2, EMG_3, EMG_4.

% Another prblem I had to address in Eitan's dataset, which I don't know if
% it will recurr wherever I'll be recording, is that there is usually a
% gap, of usually -50 to 15000 samples (one third of a second at most, that
% is) between the end of one file to the beginning of the next.

% Lastly, The algorithm for generating the initial data used by this code:
% 1. For every night, find the earliest file in which the first unit was
%    recorded.
% 2. Find the latest file in which the last unit was recorded.
% 3. Stage the entire length, between the beginning of the 5th file before 
%    the beginning of the first and the end of the 5th file after the 
%    end of the last file. for example, if one unit was recorded in the 
%    files 9-18, another in 14-38, and another during 19-39, then we start
%    staging for the night from (9-5) until (39+5)= 4 until 44.
% 4. The creating of the structure of files relevant for the
%    polysomnography data is to be done by the "Save Compact File" in the
%    Polysomnographic Data Retrieval GUI. This code will generate the
%    aforementioned "polysom_ddmmyy" file containing:
%    - FilesPaths    : A vector (1Xn) of strings with the paths of every file recorded to be staged.
%    - BeginFile     : A double, the number of the first file.
%    - EndFile       : A double, the number of the last file.
%    - RecBegin      : a 1Xn vector, where cell (1,i) is the sample in which the file i started.
%    - RecEnd        : a 1Xn vector, where cell (1,i) is the sample in which the file i ended.
%      note: samples thould be a continuum across files. i.e. if file 1
%      starts in sample 2000 and ends in sample 9000 (for exmple), then
%      the next file, file 2, should start in sample 9001. If not, and the
%      recording syste, lost samples in moving betwenn files, or didn't
%      save them, then these samples will be left blank in the staging
%      screen, i.e. if 30 sec = 7000 samles, then if it turns out that the 
%      next 30-sec period should include the gap between files, then it
%      will be samples 5000 until 12000 even if some samples in the border
%      betwenn the files are indeed abscent.
%    - polysom_analog: A 1Xn structure (n being the number of files), every
%      cell having 14 analog variables:
%      EEG_FR, EEG_FL, EEG_CR, EEG_CL, EEG_OR, EEG_OL,
%      EOG_UR, EEG_UL, EEG_LR, EEG_LL,
%      EMG_1, EMG_2, EMG_3, EMG_4.
%    - path_to_save  : A path, containing the initial polysom_ddmmyy file,
%      in which the excel and ma files with the actual staging, will be
%      saved.
%    - date: the date of the recording night, in the form ddmmyy, like in
%      the name of the "polysom_ddmmyy" file.
%    - monkey_name: a string.
%    Two other variables are to be saved post-staging, in this very file, 
%    but should be initialized in the creation of polysom_ddmmyy, so this 
%    code should receive polysom_ddmmyy with them in it:
%    - segments_already_staged - The number of 30-sec segments, out of the
%      entire length, that have been staged.
%    - Daily_Data_Cell (A cell of (1+1440,6)); the first dimension is the 
%      number of 30-second segments I will record per day, because I want 
%      the excel sheet to be for one day. 12 hours (max) of recording times
%      60 minutes times 2 times 30 sec in a minute is 1440. +1 is for the 
%      headline row, the first one. Initialized as Nan.
%    - in the future, I want to have an option where I could insert the
%      actual spikes as dots of various colors (for separate units) to the
%      staging screen where I see the analog traces, to enable me to see
%      the connection, if there is, between spikes and k complexes,
%      spindles etc.

function varargout = Polysomnography_GUI(varargin)
% POLYSOMNOGRAPHY_GUI MATLAB code for Polysomnography_GUI.fig
%      POLYSOMNOGRAPHY_GUI, by itself, creates a new POLYSOMNOGRAPHY_GUI or raises the existing
%      singleton*.
%
%      H = POLYSOMNOGRAPHY_GUI returns the handle to a new POLYSOMNOGRAPHY_GUI or the handle to
%      the existing singleton*.
%
%      POLYSOMNOGRAPHY_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in POLYSOMNOGRAPHY_GUI.M with the given input arguments.
%
%      POLYSOMNOGRAPHY_GUI('Property','Value',...) creates a new POLYSOMNOGRAPHY_GUI or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before Polysomnography_GUI_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to Polysomnography_GUI_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help Polysomnography_GUI

% Last Modified by GUIDE v2.5 15-Oct-2016 18:05:38

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @Polysomnography_GUI_OpeningFcn, ...
                   'gui_OutputFcn',  @Polysomnography_GUI_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
% End initialization code - DO NOT EDIT

function Polysomnography_GUI_OpeningFcn(hObject, eventdata, handles, varargin)
handles.output = hObject;
folder_for_staging=varargin{1,1};
path_to_save=varargin{1,2};
matfiles=varargin{1,3};
if folder_for_staging(1)~=0 && isdir(folder_for_staging) && ~isempty(dir([folder_for_staging,'\*.mat']))
    if strcmp(path_to_save,'Original Data Directory')
        path_to_save=folder_for_staging;
    end
end
handles.path_to_save=path_to_save;
handles.folder_for_staging=folder_for_staging;
handles.matfiles=matfiles;
guidata(hObject, handles);

function varargout = Polysomnography_GUI_OutputFcn(hObject, eventdata, handles) 
varargout{1} = handles.output;

% Introduction of GUI-generated variables into the GUI
Overwrite_Previous_bool=0;
Create_New_bool=0;
Staging_Protocol_bool=0;
Summary_Graph_bool=0;  

handles.Overwrite_Previous_bool=Overwrite_Previous_bool;
handles.Create_New_bool=Create_New_bool;
handles.Staging_Protocol_bool=Staging_Protocol_bool;
handles.Summary_Graph_bool=Summary_Graph_bool;

guidata(hObject,handles);

function Stage_And_Move_Callback(hObject, eventdata, handles)
% This callback sould:
% 1. Create a new Excel file OR change an existing one
% 2. In the Excel file, in only one book, create the following table:
%   date_of_recording | starting file | ending file | No. of sample on starting file | No. of sample on ending file | Stage
%   ...                 ...             ...           ...                       ...                     ...  
% 3. The excel file should be called Staging_For_Date_dd/mm/yyyy, in a pre-
%    specified path (path_to_save).
% 4. The callback will receive the user's staging and insert it in the
%    adequate place in the excel sheet as well as in a specific matlab file, 
%    and the other required data as well.
% 5. After updating and saving, the function will change the view to the
%    next 30 second screen
    
BeginFile=handles.BeginFile;
EndFile=handles.EndFile;
FilesPaths=handles.FilesPaths;
RecBegin=handles.RecBegin;
RecEnd=handles.RecEnd;
recording_date=handles.recording_date;
monkey_name=handles.monkey_name;
global polysomnography
global bit_resolution
global measure_box
global num_stds_EOG
global SEF
global SEF_1 %%^^%%
bit_resolution=1.9073; % uV per value
global SA_1 SA_2 SA_3 SA_4 SA_5 SA_6 SA_7
polysom_analog=handles.polysom_analog;
Daily_Data_Cell=handles.Daily_Data_Cell;
path_to_save=handles.path_to_save;
segments_already_staged=handles.segments_already_staged;
segments_tracings=handles.segments_tracings;
fs=handles.fs;
sec_2_stage=handles.sec_2_stage;
samples_per_segment=handles.samples_per_segment;
total_segment_num=handles.total_segment_num;
extra=handles.extra;
global insert_midstage
Daily_Data_Cell{segments_already_staged+2,1}=recording_date; % Only to check
Daily_Data_Cell{segments_already_staged+2,2}=BeginFile+segments_tracings(3,1+segments_already_staged-extra)-1;
Daily_Data_Cell{segments_already_staged+2,3}=BeginFile+segments_tracings(4,1+segments_already_staged-extra)-1;
Daily_Data_Cell{segments_already_staged+2,4}=segments_tracings(5,1+segments_already_staged-extra);
Daily_Data_Cell{segments_already_staged+2,5}=segments_tracings(6,1+segments_already_staged-extra);
Daily_Data_Cell{segments_already_staged+2,6}=handles.Choose_Stage_string;

midStage=0;
if  (insert_midstage && isfield(handles, 'Start_Sec_mid') &&  isfield(handles, 'Choose_Stage_mid'))
    Start_Sec_mid=handles.Start_Sec_mid; 
    Choose_Stage_mid=handles.Choose_Stage_mid;
    segments_already_staged=segments_already_staged+1;
    rel_sample=floor((segments_tracings(2,segments_already_staged-extra)-segments_tracings(1,segments_already_staged-extra)+1)/30*Start_Sec_mid)+Daily_Data_Cell{segments_already_staged+2-1,4}-1;
    Daily_Data_Cell{segments_already_staged+2,1}=Daily_Data_Cell{segments_already_staged+2-1,1};
    Daily_Data_Cell{segments_already_staged+2,2}=find(RecBegin(1,Daily_Data_Cell{segments_already_staged+2-1,2}-Daily_Data_Cell{2,2}+1)-1+rel_sample<=RecEnd,1,'first')+Daily_Data_Cell{2,2}-1;
    Daily_Data_Cell{segments_already_staged+2,3}=Daily_Data_Cell{segments_already_staged+2-1,3};
    Daily_Data_Cell{segments_already_staged+2,4}=rel_sample;
    Daily_Data_Cell{segments_already_staged+2,5}=Daily_Data_Cell{segments_already_staged+2-1,5};
    Daily_Data_Cell{segments_already_staged+2,6}=Choose_Stage_mid;
    
    Daily_Data_Cell{segments_already_staged+2-1,3}=find(RecBegin(1,Daily_Data_Cell{segments_already_staged+2-1,2}-Daily_Data_Cell{2,2}+1)-1+(rel_sample-1)<=RecEnd,1,'first')+Daily_Data_Cell{2,2}-1;
    Daily_Data_Cell{segments_already_staged+2-1,5}=rel_sample-1;
    extra=extra+1;
    midStage=1;
    total_segment_num=total_segment_num+1;
end

if mod(segments_already_staged,50)==0 || segments_already_staged==total_segment_num-1 || (extra==1 && segments_already_staged==1)
    if (segments_already_staged==0 || (extra==1 && segments_already_staged==1))
        mkdir(handles.path_to_save,['Staging Folder for Date ',recording_date]); % Create a new Folder to which the data is saved - both excel file and var file.
    end
    excel_file_name=['Staging for date_',recording_date]; 
    xlswrite([handles.path_to_save,'\Staging Folder for Date ',recording_date,'\',excel_file_name,'.xls'],Daily_Data_Cell); % save to excel
    save([handles.path_to_save,'\Staging Folder for Date ',recording_date,'\',excel_file_name], 'Daily_Data_Cell'); % save to cell
    save([handles.path_to_save,'\polysom_',recording_date,'.mat'],'BeginFile','EndFile','FilesPaths','RecBegin','RecEnd','polysom_analog','recording_date','monkey_name','path_to_save','Daily_Data_Cell','segments_already_staged');
end

% display data
switch midStage
    case 0
%         disp(['begin file is:
%         ',num2str(Daily_Data_Cell{segments_already_staged+2,2})]) %         commented 101716
%         disp(['end file is: ',num2str(Daily_Data_Cell{segments_already_staged+2,3})])
          disp(['   '])        % added 101716
          disp(['number of staged segments is: ',num2str(segments_already_staged+1)])
%         disp(['start sample is: ',num2str(Daily_Data_Cell{segments_already_staged+2,4})])
%         disp(['end sample is: ',num2str(Daily_Data_Cell{segments_already_staged+2,5})])
        disp(['Chosen stage is: ',handles.Choose_Stage_string])
    case 1
        disp('In this Section two segments were staged.')
        disp('1st Stage:')
        disp(['begin file is: ',num2str(Daily_Data_Cell{segments_already_staged+1,2})])
        disp(['end file is: ',num2str(Daily_Data_Cell{segments_already_staged+1,3})])
        disp(['number of staged segments is: ',num2str(segments_already_staged)])
        disp(['start sample is: ',num2str(Daily_Data_Cell{segments_already_staged+1,4})])
        disp(['end sample is: ',num2str(Daily_Data_Cell{segments_already_staged+1,5})])
        disp(['Chosen stage is: ',handles.Choose_Stage_string])

        disp('2nd Stage:')
        disp(['begin file is: ',num2str(Daily_Data_Cell{segments_already_staged+2,2})])
        disp(['end file is: ',num2str(Daily_Data_Cell{segments_already_staged+2,3})])
        disp(['number of staged segments is: ',num2str(segments_already_staged+1)])
        disp(['start sample is: ',num2str(Daily_Data_Cell{segments_already_staged+2,4})])
        disp(['end sample is: ',num2str(Daily_Data_Cell{segments_already_staged+2,5})])
        disp(['Chosen stage is: ',Choose_Stage_mid])        
end

handles.Daily_Data_Cell=Daily_Data_Cell;

if segments_already_staged==total_segment_num-1
    errordlg('All segments are Staged','Error'); 
     if handles.update_ddmmyy==1
         name=[handles.path_to_ddmmyy_string,'\Data_',recording_date,'.mat'];
         load(name);
         num_of_stages=6; stage_matrix={'W','N1','N2','N3','R','U'}; % changed 102715 (added U)
         for stage=1:num_of_stages
             relevant_indices=find(strcmp(stage_matrix(1,stage),Daily_Data_Cell(2:end,6)));
             ind=find(diff(relevant_indices)==1)+1;
             openers=relevant_indices;
             openers(ind)=[];
             ind=find(diff(relevant_indices)==1);
             closers=relevant_indices;
             closers(ind)=[];
             stage_openers=nan(1,length(openers));
             for pp=1:length(openers)
                 stage_openers(1,pp)=datafiles(Daily_Data_Cell{openers(pp)+1,2}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{openers(pp)+1,4}/fs*(generalInfo.SAMPLING*1000));
             end
             stage_closers=nan(1,length(closers));
             for ss=1:length(closers)
                 stage_closers(1,ss)=datafiles(Daily_Data_Cell{closers(ss)+1,3}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{closers(ss)+1,5}/fs*(generalInfo.SAMPLING*1000));
             end
             eval(['polysomnography.stages.',stage_matrix{1,stage},'.begin_sample=stage_openers;']);
             eval(['polysomnography.stages.',stage_matrix{1,stage},'.end_sample=stage_closers;']);
             disp(stage)
         end
     segments_already_staged=segments_already_staged+1;
     handles.segments_already_staged=segments_already_staged;

     polysomnography.Staging_date=handles.date_of_staging;
     polysomnography.Stager=handles.stager_name_string;
     
     % create *all* segments openers and closers
     all_begin=Daily_Data_Cell{2:end,4}; all_end=Daily_Data_Cell{2:end,5};
     for aa=1:size(Daily_Data_Cell,1)-1
         all_begin(aa,1)=datafiles(Daily_Data_Cell{aa+1,2}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{aa+1,4}/fs*(generalInfo.SAMPLING*1000));
         all_end(aa,1)=datafiles(Daily_Data_Cell{aa+1,3}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{aa+1,5}/fs*(generalInfo.SAMPLING*1000));
     end
     %correction that is due to the "dead samples" that are created because
     %of /fs*(generalInfo.SAMPLING*1000)) - because
     %fs<generalInfo.SAMPLING*1000, a group of undefined "dead samples"
     %occur. added 231014 MIGHT DO PROBLEMS IN neuron profile
     all_begin_new=all_begin;
     all_end_new=all_end;
     all_begin_new(1,1)=datafiles(Daily_Data_Cell{1+1,2}-BeginFile+1,1).start_sample;
     all_end_new(end,1)=datafiles(Daily_Data_Cell{end,2}-BeginFile+1,1).end_sample;
     
     for runner=1:length(all_begin_new)-1
         midpoint=round((all_end_new(runner,1)+all_begin_new(runner+1,1))/2);
         all_end_new(runner,1)=midpoint;
         all_begin_new(runner+1,1)=midpoint+1;
     end
     
     A=polysomnography.stages.W.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.W.begin_sample=A; clear A;
     A=polysomnography.stages.W.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.W.begin_sample=A; clear A;
     A=polysomnography.stages.W.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.W.end_sample=A; clear A;
     A=polysomnography.stages.W.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.W.end_sample=A; clear A;
     
     A=polysomnography.stages.N1.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N1.begin_sample=A; clear A;
     A=polysomnography.stages.N1.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N1.begin_sample=A; clear A;
     A=polysomnography.stages.N1.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N1.end_sample=A; clear A;
     A=polysomnography.stages.N1.end_sample; A=changem(A,all_end_new,all_end);  polysomnography.stages.N1.end_sample=A; clear A;
     
     A=polysomnography.stages.N2.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N2.begin_sample=A; clear A;
     A=polysomnography.stages.N2.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N2.begin_sample=A; clear A;
     A=polysomnography.stages.N2.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N2.end_sample=A; clear A;
     A=polysomnography.stages.N2.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N2.end_sample=A; clear A;
     
     A=polysomnography.stages.N3.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N3.begin_sample=A; clear A;
     A=polysomnography.stages.N3.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N3.begin_sample=A; clear A;
     A=polysomnography.stages.N3.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N3.end_sample=A; clear A;
     A=polysomnography.stages.N3.end_sample; A=changem(A,all_end_new,all_end);polysomnography.stages.N3.end_sample=A; clear A;
     
     A=polysomnography.stages.R.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.R.begin_sample=A; clear A;
     A=polysomnography.stages.R.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.R.begin_sample=A; clear A;
     A=polysomnography.stages.R.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.R.end_sample=A; clear A;
     A=polysomnography.stages.R.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.R.end_sample=A; clear A;

     A=polysomnography.stages.U.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.U.begin_sample=A; clear A;
     A=polysomnography.stages.U.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.U.begin_sample=A; clear A;
     A=polysomnography.stages.U.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.U.end_sample=A; clear A;
     A=polysomnography.stages.U.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.U.end_sample=A; clear A;     
     
     all_begin=all_begin_new; clear all_begin_new
     all_end=all_end_new; clear all_end_new
     polysomnography.all_beginnings=all_begin;
     polysomnography.all_ends=all_end;

    % compute general statistics for the session
    % Total Sleep Time
    polysomnography.Stats.Total_Sleep_Time=(sum(polysomnography.stages.N1.end_sample-polysomnography.stages.N1.begin_sample)+sum(polysomnography.stages.N2.end_sample-polysomnography.stages.N2.begin_sample)+sum(polysomnography.stages.N3.end_sample-polysomnography.stages.N3.begin_sample)+sum(polysomnography.stages.R.end_sample-polysomnography.stages.R.begin_sample))/44000;
    % Sleep Efficacy
    polysomnography.Stats.Sleep_Efficacy=polysomnography.Stats.Total_Sleep_Time/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000);
    % Sleep Fragmentation
    % I'm still not sure as to how to compute this
    % REM Latency
    polysomnography.Stats.REM_Latency=polysomnography.stages.R.begin_sample(1,1)-min([polysomnography.stages.W.begin_sample(1,1),polysomnography.stages.N1.begin_sample(1,1),polysomnography.stages.N2.begin_sample(1,1),polysomnography.stages.N3.begin_sample(1,1),polysomnography.stages.R.begin_sample(1,1)])/44000;
    % Total REM (precentage)
    polysomnography.Stats.Total_REM=(sum(polysomnography.stages.R.end_sample-polysomnography.stages.R.begin_sample)/44000)/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000);
    % Total SWS
    polysomnography.Stats.Total_SWS=(sum(polysomnography.Events.SWS.end-polysomnography.Events.SWS.start)/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000));
    % And see some more in my Email
     
     save (name,'electrode','generalInfo','datafiles','polysomnography','LickBlinkData');
     end
end

    % compute measures and classification and save them
    % EEG power ratio
    rms_thresh_EEG=10; % threshold rms for "non-real" EEG signals.
    num_of_EEG_segs=0;
    window_length=2^12;
    ratio_1=0;
    ratio_2=0;
    ratio_3=0; % added Aviv 101316 to create a better ration. This will have the ratio of the SEF to the lows
    SEF_50=0;
    SEF_result=zeros(10,1); %%^^%%
    name=[handles.path_to_ddmmyy_string,'\Data_',recording_date,'.mat'];
    load(name,'LickBlinkData');
    if rms(SA_1)>=rms_thresh_EEG
        [pxx1,f1] = pwelch(SA_1,hann(window_length),window_length/2,window_length,fs);
        % power ratio 1
        lows1_1=nanmean(pxx1(f1>=0.1 & f1<=5)); % GPe
        highs1_1=nanmean(pxx1(f1>=20 & f1<=30)); % GPe
%         lows1_1=nanmean(pxx1((f1>=0.1 & f1<=5) | (f1>=10 & f1<=17))); % striatum
%         highs1_1=nanmean(pxx1(f1>=17 & f1<=25));  % striatum     
        % power ratio 2
        lows1_2=nanmean(pxx1(f1>=0.05 & f1<=15));
        highs1_2=nanmean(pxx1(f1>=20 & f1<=35));        
        ratio_1=ratio_1+highs1_1/lows1_1;        
        ratio_2=ratio_2+highs1_2/lows1_2;
        % SEF_50
        pxx1=pxx1/sum(pxx1); % normalize power spectra
        edge_freq=SEF*sum(pxx1);
        SEF_50=SEF_50+f1(find(cumsum(pxx1)>=edge_freq,1,'first'));
        ratio_3=ratio_3+f1(find(cumsum(pxx1)>=edge_freq,1,'first'))/lows1_1; % added Aviv 101316
        e1=SEF_1(1,1)*sum(pxx1); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f1(find(cumsum(pxx1)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx1); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f1(find(cumsum(pxx1)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx1); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f1(find(cumsum(pxx1)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx1); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f1(find(cumsum(pxx1)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx1); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f1(find(cumsum(pxx1)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx1); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f1(find(cumsum(pxx1)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx1); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f1(find(cumsum(pxx1)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx1); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f1(find(cumsum(pxx1)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx1); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f1(find(cumsum(pxx1)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx1); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f1(find(cumsum(pxx1)>=e10,1,'first')); %%^^%%        
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_2)>=rms_thresh_EEG
        [pxx2,f2] = pwelch(SA_2,hann(window_length),window_length/2,window_length,fs);
        lows2_1=nanmean(pxx2(f2>=0.1 & f2<=5)); % GPe
        highs2_1=nanmean(pxx2(f2>=20 & f2<=30)); % GPe
%         lows2_1=nanmean(pxx2((f2>=0.1 & f2<=5) | (f2>=10 & f2<=17))); % striatum
%         highs2_1=nanmean(pxx2(f2>=17 & f2<=25));  % striatum      
        lows2_2=nanmean(pxx2(f2>=0.05 & f2<=15));
        highs2_2=nanmean(pxx2(f2>=20 & f2<=35));        
        ratio_1=ratio_1+highs2_1/lows2_1;
        ratio_2=ratio_2+highs2_2/lows2_2;
        % SEF_50
        pxx2=pxx2/sum(pxx2); % normalize power spectra
        edge_freq=SEF*sum(pxx2);
        SEF_50=SEF_50+f2(find(cumsum(pxx2)>=edge_freq,1,'first')); 
        ratio_3=ratio_3+f2(find(cumsum(pxx2)>=edge_freq,1,'first'))/lows2_1; % added Aviv 101316
        e1=SEF_1(1,1)*sum(pxx2); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f2(find(cumsum(pxx2)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx2); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f2(find(cumsum(pxx2)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx2); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f2(find(cumsum(pxx2)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx2); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f2(find(cumsum(pxx2)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx2); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f2(find(cumsum(pxx2)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx2); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f2(find(cumsum(pxx2)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx2); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f2(find(cumsum(pxx2)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx2); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f2(find(cumsum(pxx2)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx2); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f2(find(cumsum(pxx2)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx2); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f2(find(cumsum(pxx2)>=e10,1,'first')); %%^^%%           
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_3)>=rms_thresh_EEG
        [pxx3,f3] = pwelch(SA_3,hann(window_length),window_length/2,window_length,fs);
        lows3_1=nanmean(pxx3(f3>=0.1 & f3<=5)); % GPe
        highs3_1=nanmean(pxx3(f3>=20 & f3<=30)); % GPe
%         lows3_1=nanmean(pxx3((f3>=0.1 & f3<=5) | (f3>=10 & f3<=17))); % striatum
%         highs3_1=nanmean(pxx3(f3>=17 & f3<=25)); % striatum       
        lows3_2=nanmean(pxx3(f3>=0.05 & f3<=15));
        highs3_2=nanmean(pxx3(f3>=20 & f3<=35));
        ratio_1=ratio_1+highs3_1/lows3_1;        
        ratio_2=ratio_2+highs3_2/lows3_2;
        % SEF_50
        pxx3=pxx3/sum(pxx3); % normalize power spectra
        edge_freq=SEF*sum(pxx3);
        SEF_50=SEF_50+f3(find(cumsum(pxx3)>=edge_freq,1,'first'));   
        ratio_3=ratio_3+f3(find(cumsum(pxx3)>=edge_freq,1,'first'))/lows3_1; % added Aviv 101316
        e1=SEF_1(1,1)*sum(pxx3); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f3(find(cumsum(pxx3)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx3); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f3(find(cumsum(pxx3)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx3); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f3(find(cumsum(pxx3)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx3); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f3(find(cumsum(pxx3)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx3); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f3(find(cumsum(pxx3)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx3); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f3(find(cumsum(pxx3)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx3); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f3(find(cumsum(pxx3)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx3); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f3(find(cumsum(pxx3)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx3); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f3(find(cumsum(pxx3)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx3); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f3(find(cumsum(pxx3)>=e10,1,'first')); %%^^%%           
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_4)>=rms_thresh_EEG
        [pxx4,f4] = pwelch(SA_4,hann(window_length),window_length/2,window_length,fs);
        lows4_1=nanmean(pxx4(f4>=0.1 & f4<=5)); % GPe
        highs4_1=nanmean(pxx4(f4>=20 & f4<=30)); % GPe
%         lows4_1=nanmean(pxx4((f4>=0.1 & f4<=5) | (f4>=10 & f4<=17))); % striatum
%         highs4_1=nanmean(pxx4(f4>=17 & f4<=25)); % striatum       
        lows4_2=nanmean(pxx4(f4>=0.05 & f4<=15));
        highs4_2=nanmean(pxx4(f4>=20 & f4<=35));
        ratio_1=ratio_1+highs4_1/lows4_1;        
        ratio_2=ratio_2+highs4_2/lows4_2;        
        % SEF_50
        pxx4=pxx4/sum(pxx4); % normalize power spectra
        edge_freq=SEF*sum(pxx4);
        SEF_50=SEF_50+f4(find(cumsum(pxx4)>=edge_freq,1,'first')); 
        ratio_3=ratio_3+f4(find(cumsum(pxx4)>=edge_freq,1,'first'))/lows4_1; % added Aviv 101316
        e1=SEF_1(1,1)*sum(pxx4); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f4(find(cumsum(pxx4)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx4); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f4(find(cumsum(pxx4)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx4); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f4(find(cumsum(pxx4)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx4); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f4(find(cumsum(pxx4)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx4); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f4(find(cumsum(pxx4)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx4); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f4(find(cumsum(pxx4)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx4); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f4(find(cumsum(pxx4)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx4); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f4(find(cumsum(pxx4)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx4); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f4(find(cumsum(pxx4)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx4); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f4(find(cumsum(pxx4)>=e10,1,'first')); %%^^%%           
        num_of_EEG_segs=num_of_EEG_segs+1;
    end    
    ratio_1=ratio_1/num_of_EEG_segs;
    ratio_2=ratio_2/num_of_EEG_segs;
    ratio_3=ratio_3/num_of_EEG_segs;
    SEF_50=SEF_50/num_of_EEG_segs;
    SEF_result=SEF_result/num_of_EEG_segs; %%^^%%
    if segments_already_staged==0
        measure_box(1,1)=log10(ratio_1);
    else
        measure_box(size(measure_box,1)+1,1)=log10(ratio_1);
    end
    
    measure_box(end,2)=log10(ratio_2);
    
    measure_box(end,3)=SEF_50;
    measure_box(end,8)=SEF_result(1,1); %%^^%%
    measure_box(end,9)=SEF_result(2,1); %%^^%%
    measure_box(end,10)=SEF_result(3,1); %%^^%%
    measure_box(end,11)=SEF_result(4,1); %%^^%%
    measure_box(end,12)=SEF_result(5,1); %%^^%%
    measure_box(end,13)=SEF_result(6,1); %%^^%%
    measure_box(end,14)=SEF_result(7,1); %%^^%%
    measure_box(end,15)=SEF_result(8,1); %%^^%%
    measure_box(end,16)=SEF_result(9,1); %%^^%%
    measure_box(end,17)=SEF_result(10,1); %%^^%%
    measure_box(end,18)=log10(ratio_3); %%^^%% Added Aviv 101316
        
    % EOG events number
    num_of_EOG_segs=0;
    rms_thresh_EOG=10;

    crossnum=0;
    if rms(SA_5)>=rms_thresh_EOG
        all_above_points1=find(SA_5>=nanmean(SA_5)+std(SA_5)*num_stds_EOG);
        all_below_points1=find(SA_5<=nanmean(SA_5)-std(SA_5)*num_stds_EOG);
        crossnum=crossnum+length(find(diff(all_above_points1)>1))+1; % +1 for the ever first one that is missed by this method
        crossnum=crossnum+length(find(diff(all_below_points1)>1))+1;
        num_of_EOG_segs=num_of_EOG_segs+1;
    end
    if rms(SA_6)>=rms_thresh_EOG
        all_above_points2=find(SA_6>=nanmean(SA_6)+std(SA_6)*num_stds_EOG);
        all_below_points2=find(SA_6<=nanmean(SA_6)-std(SA_6)*num_stds_EOG);
        crossnum=crossnum+length(find(diff(all_above_points2)>1))+1; % +1 for the ever first one that is missed by this method
        crossnum=crossnum+length(find(diff(all_below_points2)>1))+1;
        num_of_EOG_segs=num_of_EOG_segs+1;
    end    
    
    measure_box(end,4)=crossnum/num_of_EOG_segs;
    
    % EMG root mean square
    measure_box(end,5)=rms(SA_7);
    
    % fraction of eye_open (=1)
    load(name,'datafiles');
    file_num=Daily_Data_Cell{segments_already_staged+2,2}-Daily_Data_Cell{2,2}+1;
    segment_start=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,4}-1)*(44000/fs)+1;
    segment_end=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,5})*(44000/fs);
    if segment_end<segment_start
        file_num=file_num+1;
        segment_end=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,5})*(44000/fs);
    end        
    relevant_blink=LickBlinkData.EyeState(find(LickBlinkData.SNR_TimeStamp>=segment_start & LickBlinkData.SNR_TimeStamp<=segment_end));
    if ~isnan((length(relevant_blink)-sum(relevant_blink))/length(relevant_blink))
        measure_box(end,6)=(length(relevant_blink)-sum(relevant_blink))/length(relevant_blink);
    else
        measure_box(end,6)=1.5;
    end
    
    if strcmp(handles.Choose_Stage_string,'W')
        measure_box(end,7)=1; % wake
    elseif strcmp(handles.Choose_Stage_string,'N2')
        measure_box(end,7)=2; % NREM
    elseif strcmp(handles.Choose_Stage_string,'R')
        measure_box(end,7)=3; % REM
    elseif strcmp(handles.Choose_Stage_string,'Unclassified')
        measure_box(end,7)=0; % Unclassified        
    end
    
%     disp(['RATIO 1 is: ',num2str(measure_box(end,1))]); % commented 101716
%     disp(['RATIO 2 is: ',num2str(measure_box(end,2))]); % commented 101716
%     disp(['RATIO 3 is: ',num2str(measure_box(end,18))]); % Added Aviv 101316 commented 101816
    disp(['SEF-70 is: ',num2str(measure_box(end,13))]);
    disp(['EMG rms is: ',num2str(measure_box(end,5))]);
    disp(['Fraction of eye open is: ',num2str(measure_box(end,6))]);

if segments_already_staged~=total_segment_num
    segments_already_staged=segments_already_staged+1; % every time this button is pushed, a segment was staged.
    handles.segments_already_staged=segments_already_staged;

    EEG_trace_1=handles.EEG_trace_1;
    EEG_trace_2=handles.EEG_trace_2;
    EEG_trace_3=handles.EEG_trace_3;
    EEG_trace_4=handles.EEG_trace_4;

    EOG_trace_1=handles.EOG_trace_1;
    EOG_trace_2=handles.EOG_trace_2;

    EMG_trace_1=handles.EMG_trace_1;

    for all_plotter=1:7
        if all_plotter<5 % EEG data traces
            if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
                eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EEG_trace_',num2str(all_plotter),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
            end
            trace_to_show=double(trace_to_show);
            eval(['plot(handles.EEG_',num2str(all_plotter),',1:(segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1),BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP''),''k'');']);
            eval(['axes(handles.EEG_',num2str(all_plotter),');']);
            ylabel('Voltage \muV')
            ylim([-250 250])
            % save trace for Spectral Analysis 
            trace_to_show=double(trace_to_show);
            eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);
        elseif (all_plotter<7 && all_plotter>=5)
            if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
                eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EOG_trace_',num2str(all_plotter-4),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
            end
%             if all_plotter==6
%                 trace_to_show=-trace_to_show;
%             end
            trace_to_show=double(trace_to_show);
            eval(['plot(handles.EOG_',num2str(all_plotter-4),',1:(segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1),BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP''),''b'');']);
            eval(['axes(handles.EOG_',num2str(all_plotter-4),');']);
            if all_plotter==5 %from here changed 100115
                hold on;
                plot(1:length(SA_5),(nanmean(SA_5)+num_stds_EOG*std(SA_5))*ones(1,length(SA_5)),'r','LineStyle','--');
                hold on;
                plot(1:length(SA_5),(nanmean(SA_5)-num_stds_EOG*std(SA_5))*ones(1,length(SA_5)),'r','LineStyle','--');
            elseif all_plotter==6
                hold on;
                plot(1:length(SA_6),(nanmean(SA_6)+num_stds_EOG*std(SA_6))*ones(1,length(SA_6)),'r','LineStyle','--');
                hold on;
                plot(1:length(SA_6),(nanmean(SA_6)-num_stds_EOG*std(SA_6))*ones(1,length(SA_6)),'r','LineStyle','--');                
            end %until here changed 100115
            hold off;            
            ylabel('Voltage \muV')
            ylim([-400 400])
            % save trace for Spectral Analysis 
            trace_to_show=double(trace_to_show);
            eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);            
        elseif all_plotter==7
            if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
            elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
               eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EMG_trace_',num2str(all_plotter-6),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']); 
            end
            trace_to_show=double(trace_to_show);
            eval(['plot(handles.EMG_',num2str(all_plotter-6),',1:(segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1),BWfilter_designer(trace_to_show, fs, [10,100], [5,120], ''BP''),''r'');']);
            axes(handles.EMG_1);
            ylabel('Voltage \muV')
            ylim([-80 80])
            % save trace for Spectral Analysis 
            trace_to_show=double(trace_to_show);
            eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [10,100], [5,120], ''BP'');']);
        end
        
    numer=floor((segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1)/10); % changed 100115
%     set(gca,'Xtick',[1*numer 2*numer 3*numer 4*numer 5*numer 6*numer 7*numer 8*numer 9*numer 10*numer 11*numer 12*numer 13*numer 14*numer 15*numer 16*numer 17*numer 18*numer 19*numer 20*numer 21*numer 22*numer 23*numer 24*numer 25*numer 26*numer 27*numer 28*numer 29*numer 30*numer]) % changed 100115
set(gca,'Xtick',[1*numer 2*numer 3*numer 4*numer 5*numer 6*numer 7*numer 8*numer 9*numer 10*numer]) % changed 100115
%     set(gca,'XTickLabel',char('1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20','21','22','23','24','25','26','27','28','29','30')); % changed 100115
    set(gca,'XTickLabel',char('1','2','3','4','5','6','7','8','9','10')); % changed 100115
    xlim([0 segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)])
%     axis 'auto y'
    grid on
    end
end

% SAVING
% segments_already_staged and Daily_Data_Cell should be updated into the
% original polysom_ddmmyy file (that holds the initial data! note you don't
% erase the initial data and save well!!)
handles.extra=extra;
handles.total_segment_num=total_segment_num;
guidata(hObject,handles);

function Choose_EEG_1_Callback(hObject, eventdata, handles)
Choose_EEG_1=get(hObject,'Value');
switch Choose_EEG_1-1
       case 1
              EEG_trace_1='EEG_FR';
       case 2
              EEG_trace_1='EEG_FL';
       case 3
              EEG_trace_1='EEG_CR';
       case 4
              EEG_trace_1='EEG_CL';
       case 5
              EEG_trace_1='EEG_OR';
       case 6
              EEG_trace_1='EEG_OL';
       case 0 % this means the 'Choose EEG_1 Channel' option is chosen
              errordlg('Proper channel should be provided in EEG_1','Error');
end
handles.EEG_trace_1=EEG_trace_1;
guidata(hObject,handles);

function Choose_EEG_1_CreateFcn(hObject, eventdata, handles)

if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EEG_1 Channel';'Frontal(R)';'Frontal(L)';'Central(R)';'Central(L)';'Occipital(R)';'Occipital(L)';});

function Choose_EEG_2_Callback(hObject, eventdata, handles)
Choose_EEG_2=get(hObject,'Value');
switch Choose_EEG_2-1
       case 1
              EEG_trace_2='EEG_FR';
       case 2
              EEG_trace_2='EEG_FL';
       case 3
              EEG_trace_2='EEG_CR';
       case 4
              EEG_trace_2='EEG_CL';
       case 5
              EEG_trace_2='EEG_OR';
       case 6
              EEG_trace_2='EEG_OL';
       case 0 
              errordlg('Proper channel should be provided in EEG_2','Error');
end
handles.EEG_trace_2=EEG_trace_2;
guidata(hObject,handles);

function Choose_EEG_2_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EEG_2 Channel';'Frontal(R)';'Frontal(L)';'Central(R)';'Central(L)';'Occipital(R)';'Occipital(L)';});

function Choose_EEG_3_Callback(hObject, eventdata, handles)
Choose_EEG_3=get(hObject,'Value');
switch Choose_EEG_3-1
       case 1
              EEG_trace_3='EEG_FR';
       case 2
              EEG_trace_3='EEG_FL';
       case 3
              EEG_trace_3='EEG_CR';
       case 4
              EEG_trace_3='EEG_CL';
       case 5
              EEG_trace_3='EEG_OR';
       case 6
              EEG_trace_3='EEG_OL';
       case 0 
              errordlg('Proper channel should be provided in EEG_3','Error');
end
handles.EEG_trace_3=EEG_trace_3;
guidata(hObject,handles);

% --- Executes during object creation, after setting all properties.
function Choose_EEG_3_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EEG_3 Channel';'Frontal(R)';'Frontal(L)';'Central(R)';'Central(L)';'Occipital(R)';'Occipital(L)';});

function Choose_EEG_4_Callback(hObject, eventdata, handles)
Choose_EEG_4=get(hObject,'Value');
switch Choose_EEG_4-1
       case 1
              EEG_trace_4='EEG_FR';
       case 2
              EEG_trace_4='EEG_FL';
       case 3
              EEG_trace_4='EEG_CR';
       case 4
              EEG_trace_4='EEG_CL';
       case 5
              EEG_trace_4='EEG_OR';
       case 6
              EEG_trace_4='EEG_OL';
       case 0 
              errordlg('Proper channel should be provided in EEG_4','Error');
end
handles.EEG_trace_4=EEG_trace_4;
guidata(hObject,handles);

function Choose_EEG_4_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EEG_4 Channel';'Frontal(R)';'Frontal(L)';'Central(R)';'Central(L)';'Occipital(R)';'Occipital(L)';});

function Choose_EOG_1_Callback(hObject, eventdata, handles)
Choose_EOG_1=get(hObject,'Value');
switch Choose_EOG_1-1
       case 1
              EOG_trace_1='EOG_UR';
       case 2
              EOG_trace_1='EOG_UL';
       case 3
              EOG_trace_1='EOG_LR';
       case 4
              EOG_trace_1='EOG_LL';
       case 0 
              errordlg('Proper channel should be provided in EOG_1','Error');
end
handles.EOG_trace_1=EOG_trace_1;
guidata(hObject,handles);

function Choose_EOG_1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EOG_1 Channel';'Upper Right';'Upper Left';'Lower Right';'Lower Left';});

function Choose_EOG_2_Callback(hObject, eventdata, handles)
Choose_EOG_2=get(hObject,'Value');
switch Choose_EOG_2-1
       case 1
              EOG_trace_2='EOG_UR';
       case 2
              EOG_trace_2='EOG_UL';
       case 3
              EOG_trace_2='EOG_LR';
       case 4
              EOG_trace_2='EOG_LL';
       case 0 
              errordlg('Proper channel should be provided in EOG_2','Error');
end
handles.EOG_trace_2=EOG_trace_2;
guidata(hObject,handles);

function Choose_EOG_2_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EOG_2 Channel';'Upper Right';'Upper Left';'Lower Right';'Lower Left';});

function Choose_EMG_1_Callback(hObject, eventdata, handles)
Choose_EMG_1=get(hObject,'Value');
switch Choose_EMG_1-1
       case 1
              EMG_trace_1='EMG_1';
       case 2
              EMG_trace_1='EMG_2';
       case 3
              EMG_trace_1='EMG_3';
       case 4
              EMG_trace_1='EMG_4';
       case 0 
              errordlg('Proper channel should be provided in EMG_1','Error');
end
handles.EMG_trace_1=EMG_trace_1;
guidata(hObject,handles);

function Choose_EMG_1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose EMG_1 Channel';'EMG_1';'EMG_2';'EMG_3';'EMG_4';});

function Choose_Stage_Callback(hObject, eventdata, handles)
list=get(hObject,'String');
val=get(hObject,'Value');
Choose_Stage_string=list{val};
handles.Choose_Stage_string=Choose_Stage_string;
guidata(hObject,handles);

function Choose_Stage_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose Sleep Stage';'W';'N1';'N2';'N3';'R';'Unclassified'});

function Staging_Path_Callback(hObject, eventdata, handles)
Staging_Path_string=get(hObject,'String');
handles.Staging_Date_string=Staging_Path_string;

if strcmp(Staging_Path_string,'Polysom Path for Staging')
    errordlg('Path for Staging must be specified','Error');
end

guidata(hObject,handles);

function Staging_Path_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Choose_Spec_File_Callback(hObject, eventdata, handles)
Choose_Spec_File_string=get(hObject,'String');
handles.Choose_Spec_File_string=Choose_Spec_File_string;
guidata(hObject,handles);

function Choose_Spec_File_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function Overwrite_Previous_Callback(hObject, eventdata, handles)
Overwrite_Previous_bool=get(hObject,'Value');
handles.Overwrite_Previous_bool=Overwrite_Previous_bool;
guidata(hObject,handles);

function Create_New_Callback(hObject, eventdata, handles)
Create_New_bool=get(hObject,'Value');
handles.Create_New_bool=Create_New_bool;
guidata(hObject,handles);

function Staging_Protocol_Callback(hObject, eventdata, handles)
Staging_Protocol_bool=get(hObject,'Value');
handles.Staging_Protocol_bool=Staging_Protocol_bool;
guidata(hObject,handles);

function Summary_Graph_Callback(hObject, eventdata, handles)
Summary_Graph_bool=get(hObject,'Value');
handles.Summary_Graph_bool=Summary_Graph_bool;
guidata(hObject,handles);

function Load_Traces_Callback(hObject, eventdata, handles)
path_to_save=handles.path_to_save;
folder_for_staging=handles.folder_for_staging;
matfiles=handles.matfiles;
global SWS;
global bit_resolution
bit_resolution=1.9073;
global REMs
global K_complexes;
global Arousal;
global sleep_spindles;
global REM_without_atonia;
global Repetitive_Muscle;
SWS=0;
REMs=0;
K_complexes=0;
Arousal=0;
sleep_spindles=0;
REM_without_atonia=0;
Repetitive_Muscle=0;
extra=0;
insert_midstage=0;
global SA_1 SA_2 SA_3 SA_4 SA_5 SA_6 SA_7
global measure_box % changed 100115
global num_stds_EOG
global SEF
global SEF_1 %%^^%%
num_stds_EOG=2.25;
SEF=0.5;
SEF_1=[0.4,0.45,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9];
% Here, the names of the relevant file in the RAW data are filled. This
% will be constant for all the recording time, so fill only once. % CHANGE BACK 110315% CHANGE BACK 110315% CHANGE BACK 110315
% name_EEG_FL='CRAW_023.Samples'; % Nymeria % change back to name_EEG_FL='CRAW_023.Samples'; % Nymeria
name_EEG_FL='CRAW_020.Samples'; % Daenerys
name_EEG_FR='CRAW_017.Samples'; name_EEG_CR='CRAW_018.Samples'; name_EEG_CL='CRAW_021.Samples'; name_EEG_OR='CRAW_019.Samples'; name_EEG_OL='CRAW_022.Samples';
name_EOG_UR='CRAW_034.Samples'; name_EOG_UL='CRAW_034.Samples'; name_EOG_LR='CRAW_035.Samples'; name_EOG_LL='CRAW_035.Samples'; % CHANGE BACK FIRST TWO TO 34!!
name_EMG_1='CRAW_036.Samples'; name_EMG_2='CRAW_036.Samples'; name_EMG_3='CRAW_036.Samples'; name_EMG_4='CRAW_036.Samples';
name_Ground_1='CRAW_022.Samples'; name_Ground_2='CRAW_022.Samples'; % TODO NEW: NO NEED TO SUBTRACT GROUND. IT IS ALREADY SUBTRACTED.

% name_EEG_FL='CRAW_023'; % Nymeria
% % name_EEG_FL='CRAW_020'; % Daenerys
% name_EEG_FR='CRAW_017'; name_EEG_CR='CRAW_018'; name_EEG_CL='CRAW_021'; name_EEG_OR='CRAW_019'; name_EEG_OL='CRAW_022';
% name_EOG_UR='CRAW_034'; name_EOG_UL='CRAW_034'; name_EOG_LR='CRAW_035'; name_EOG_LL='CRAW_035';
% name_EMG_1='CRAW_036.'; name_EMG_2='CRAW_036'; name_EMG_3='CRAW_036'; name_EMG_4='CRAW_036';
% name_Ground_1='CRAW_022'; name_Ground_2='CRAW_022'; % TODO NEW: NO NEED TO SUBTRACT GROUND. IT IS ALREADY SUBTRACTED.

% name_EEG_FL='CLFP_023.Samples'; % Nymeria
% name_EEG_FL='CLFP_020.Samples'; % Daenerys
% name_EEG_FR='CLFP_017.Samples'; name_EEG_CR='CLFP_018.Samples'; name_EEG_CL='CLFP_021.Samples'; name_EEG_OR='CLFP_019.Samples'; name_EEG_OL='CLFP_022.Samples';  %% % CHANGE BACK 110315
% name_EOG_UR='CLFP_034.Samples'; name_EOG_UL='CLFP_034.Samples'; name_EOG_LR='CLFP_035.Samples'; name_EOG_LL='CLFP_035.Samples';
% name_EMG_1='CLFP_036.Samples'; name_EMG_2='CLFP_036.Samples'; name_EMG_3='CLFP_036.Samples'; name_EMG_4='CLFP_036.Samples';
% name_Ground_1='CLFP_022.Samples'; name_Ground_2='CLFP_022.Samples'; % TODO NEW: NO NEED TO SUBTRACT GROUND. IT IS ALREADY SUBTRACTED.
measure_box=nanmean(1,18); %%^^%% % changed 100115. This will hold both the EEG, EOG and EMG measures about the specific semgnet, and its subsequent staging (by human). Changed 101316 Aviv to allow ratio_3

FilesPaths=cell(1,length(matfiles));

    polysom_analog=struct('name',[]);
    for ii=1:length(matfiles)
        eval(['clear ',name_EEG_FR(1:end-8),' ',name_EEG_FL(1:end-8),' ',name_EEG_CR(1:end-8),' ',name_EEG_CL(1:end-8),' ',name_EEG_OR(1:end-8),' ',name_EEG_OL(1:end-8),' ',name_EOG_UR(1:end-8),' ',name_EOG_UL(1:end-8),' ',name_EOG_LR(1:end-8),' ',name_EOG_LL(1:end-8),' ',name_EMG_1(1:end-8),' ',name_EMG_2(1:end-8),' ',name_EMG_3(1:end-8),' ',name_EMG_4(1:end-8),' ',name_Ground_1(1:end-8),';']);
        filename=matfiles(ii,1).name;
        if ~isempty(filename)
            disp(ii)
            polysom_analog(1,ii).name=[folder_for_staging,'\',filename(1,:)];
            FilesPaths{1,ii}=polysom_analog(1,ii).name;
%             load(polysom_analog(1,ii).name,'CLFP*')
            load(polysom_analog(1,ii).name,'CRAW*') % CHANGE BACK 110315
            eval(['polysom_analog(1,ii).EEG_FR=',double(name_EEG_FR),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EEG_FL=',double(name_EEG_FL),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EEG_CR=',double(name_EEG_CR),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EEG_CL=',double(name_EEG_CL),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EEG_OR=',double(name_EEG_OR),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EEG_OL=',double(name_EEG_OL),'*bit_resolution;']);             
%             eval(['polysom_analog(1,ii).EOG_UR=',name_EOG_UR,'*bit_resolution+',name_EOG_LR,'*bit_resolution;']);  % to make the EOG signal cleaner.
%             eval(['polysom_analog(1,ii).EOG_UL=',name_EOG_UL,'*bit_resolution+',name_EOG_LL,'*bit_resolution;']); 
%             eval(['polysom_analog(1,ii).EOG_LR=',name_EOG_LR,'*bit_resolution+',name_EOG_UR,'*bit_resolution;']); 
%             eval(['polysom_analog(1,ii).EOG_LL=',name_EOG_LL,'*bit_resolution+',name_EOG_UL,'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EOG_UR=',double(name_EOG_UR),'*bit_resolution;']);  
            eval(['polysom_analog(1,ii).EOG_UL=',double(name_EOG_UL),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EOG_LR=',double(name_EOG_LR),'*bit_resolution;']); 
            eval(['polysom_analog(1,ii).EOG_LL=',double(name_EOG_LL),'*bit_resolution;']); 

            eval(['polysom_analog(1,ii).EMG_1=',double(name_EMG_1),'*bit_resolution;']);
            eval(['polysom_analog(1,ii).EMG_2=',double(name_EMG_2),'*bit_resolution;']);
            eval(['polysom_analog(1,ii).EMG_3=',double(name_EMG_3),'*bit_resolution;']);
        end;
    end;

    segments_already_staged=0;
    BeginFile=str2double(matfiles(1,1).name(end-10:end-9)); EndFile=str2double(matfiles(end,1).name(end-10:end-9));
    RecBegin=nan(1,length(matfiles)); RecEnd=nan(1,length(matfiles));
    RecBegin(1,1)=1; RecEnd(1,1)=length(polysom_analog(1,1).EEG_FR);
    for jj=2:length(matfiles)
        RecBegin(1,jj)=RecBegin(1,jj-1)+length(polysom_analog(1,jj-1).EEG_FR);
        RecEnd(1,jj)=RecEnd(1,jj-1)+length(polysom_analog(1,jj).EEG_FR);
    end
% - BeginFile: A string with the file name of the first file.
% - EndFile  : A string with the file name of the last file.
% - RecBegin : The number of the first file.
% - RecEnd   : The number of the last file.
% - polysom_analog: A 1X5 structure, every cell having 14 analog variables:
%   EEG_FR, EEG_FL, EEG_CR, EEG_CL, EEG_OR, EEG_OL,
%   EOG_UR, EEG_UL, EEG_LR, EEG_LL,
%   EMG_1, EMG_2, EMG_3, EMG_4.    
    
% % When you check if this works, you ought to make sure that when the user
% % wishes to change the electrode they inspect, the changing is updated here
% % and the user is further presented, in the future screens, the same
% % electrode they requested.

% create structure containing all the relevant files

fs=CRAW_017.KHz*1000; % CHANGE BACK TO 1000 !!!!!!! 220415!!! sampling frequency in Hz. Don't forget to fill the right one. 
% fs=CRAW_017_KHz*1000; % CHANGE BACK TO 1000 !!!!!!! 220415!!! sampling frequency in Hz. Don't forget to fill the right one. 
% fs=CLFP_017.KHz*1000;  % CHANGE BACK
sec_2_stage=10; % length of segment in seconds % changed 100115
samples_per_segment=floor(sec_2_stage*fs);

% create segments tracings (in sample number), based on the number of
% files, fs, and number of seconds per staged segment.
total_segment_num=ceil((RecEnd(1,end)-RecBegin(1,1))/samples_per_segment);

segments_tracings=nan(6,total_segment_num); % 1st row - the first sample of segmment i, 2nd row - the last sample, 3rd and 4th rows - the original data files in which each of the beginning and end samples of the segment (respectively) are. 
for kk=1:1:total_segment_num
    segments_tracings(1,kk)=RecBegin(1,1)+floor((kk-1)*samples_per_segment);
    segments_tracings(2,kk)=min(RecEnd(1,end),RecBegin(1,1)+floor(kk*samples_per_segment-1));
    segments_tracings(3,kk)=find(RecEnd>=segments_tracings(1,kk),1,'first');
    segments_tracings(4,kk)=find(RecEnd>=segments_tracings(2,kk),1,'first');
    segments_tracings(5,kk)=segments_tracings(1,kk)-RecBegin(1,segments_tracings(3,kk))+1;
    segments_tracings(6,kk)=segments_tracings(2,kk)-RecBegin(1,segments_tracings(4,kk))+1;
end

EEG_trace_1=handles.EEG_trace_1;
EEG_trace_2=handles.EEG_trace_2;
EEG_trace_3=handles.EEG_trace_3;
EEG_trace_4=handles.EEG_trace_4;

EOG_trace_1=handles.EOG_trace_1;
EOG_trace_2=handles.EOG_trace_2;

EMG_trace_1=handles.EMG_trace_1;

% Correction of the lengths of polysom_analog. Sometimes, because of saving
% issues, some of the vars in polysom_analog have more samples than the
% "official" sample count for the file. Here I trim these vars. 081516
file_lengths=RecEnd-RecBegin+1;
for zzz=1:length(polysom_analog)
    if length(polysom_analog(1,zzz).EEG_FR)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_FR=polysom_analog(1,zzz).EEG_FR(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_FR)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_FR=[polysom_analog(1,zzz).EEG_FR,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_FR))];        
    end
    if length(polysom_analog(1,zzz).EEG_FL)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_FL=polysom_analog(1,zzz).EEG_FL(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_FL)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_FL=[polysom_analog(1,zzz).EEG_FL,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_FL))];                
    end 
    if length(polysom_analog(1,zzz).EEG_CR)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_CR=polysom_analog(1,zzz).EEG_CR(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_CR)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_CR=[polysom_analog(1,zzz).EEG_CR,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_CR))];                
    end
    if length(polysom_analog(1,zzz).EEG_CL)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_CL=polysom_analog(1,zzz).EEG_CL(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_CL)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_CL=[polysom_analog(1,zzz).EEG_CL,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_CL))];                
    end    
    if length(polysom_analog(1,zzz).EEG_OR)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_OR=polysom_analog(1,zzz).EEG_OR(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_OR)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_OR=[polysom_analog(1,zzz).EEG_OR,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_OR))];                        
    end    
    if length(polysom_analog(1,zzz).EEG_OL)>file_lengths(zzz)
        polysom_analog(1,zzz).EEG_OL=polysom_analog(1,zzz).EEG_OL(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EEG_OL)<file_lengths(zzz)
        polysom_analog(1,zzz).EEG_OL=[polysom_analog(1,zzz).EEG_OL,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EEG_OL))];                        
    end
    if length(polysom_analog(1,zzz).EOG_UR)>file_lengths(zzz)
        polysom_analog(1,zzz).EOG_UR=polysom_analog(1,zzz).EOG_UR(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EOG_UR)<file_lengths(zzz)
        polysom_analog(1,zzz).EOG_UR=[polysom_analog(1,zzz).EOG_UR,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EOG_UR))];                        
    end  
    if length(polysom_analog(1,zzz).EOG_UL)>file_lengths(zzz)
        polysom_analog(1,zzz).EOG_UL=polysom_analog(1,zzz).EOG_UL(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EOG_UL)<file_lengths(zzz)
        polysom_analog(1,zzz).EOG_UL=[polysom_analog(1,zzz).EOG_UL,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EOG_UL))];                        
    end
    if length(polysom_analog(1,zzz).EOG_LR)>file_lengths(zzz)
        polysom_analog(1,zzz).EOG_LR=polysom_analog(1,zzz).EOG_LR(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EOG_LR)<file_lengths(zzz)
        polysom_analog(1,zzz).EOG_LR=[polysom_analog(1,zzz).EOG_LR,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EOG_LR))];                        
    end  
    if length(polysom_analog(1,zzz).EOG_LL)>file_lengths(zzz)
        polysom_analog(1,zzz).EOG_LL=polysom_analog(1,zzz).EOG_LL(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EOG_LL)<file_lengths(zzz)
        polysom_analog(1,zzz).EOG_LL=[polysom_analog(1,zzz).EOG_LL,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EOG_LL))];                        
    end
    if length(polysom_analog(1,zzz).EMG_1)>file_lengths(zzz)
        polysom_analog(1,zzz).EMG_1=polysom_analog(1,zzz).EMG_1(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EMG_1)<file_lengths(zzz)
        polysom_analog(1,zzz).EMG_1=[polysom_analog(1,zzz).EMG_1,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EMG_1))];                        
    end    
    if length(polysom_analog(1,zzz).EMG_2)>file_lengths(zzz)
        polysom_analog(1,zzz).EMG_2=polysom_analog(1,zzz).EMG_2(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EMG_2)<file_lengths(zzz)
        polysom_analog(1,zzz).EMG_2=[polysom_analog(1,zzz).EMG_2,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EMG_2))];                                
    end    
    if length(polysom_analog(1,zzz).EMG_3)>file_lengths(zzz)
        polysom_analog(1,zzz).EMG_3=polysom_analog(1,zzz).EMG_3(1:file_lengths(zzz));
    elseif length(polysom_analog(1,zzz).EMG_3)<file_lengths(zzz)
        polysom_analog(1,zzz).EMG_3=[polysom_analog(1,zzz).EMG_3,nan(1,file_lengths(zzz)-length(polysom_analog(1,zzz).EMG_3))];                                
    end 
end
% end 081516
    
for all_plotter=1:7
    if all_plotter<5 % EEG data traces
        if segments_tracings(3,1+segments_already_staged)==segments_tracings(4,1+segments_already_staged)
            eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged):segments_tracings(6,1+segments_already_staged));']);
        else
           eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged)).(EEG_trace_',num2str(all_plotter),')(1,1:segments_tracings(6,1+segments_already_staged))];']);
        end
        trace_to_show=double(trace_to_show);
        eval(['plot(handles.EEG_',num2str(all_plotter),',1:(segments_tracings(2,1+segments_already_staged)-segments_tracings(1,1+segments_already_staged)+1),BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP''),''k'');']);
        eval(['axes(handles.EEG_',num2str(all_plotter),');']);
        ylim([-200 200])  
        ylabel('Voltage \muV')
        % save trace for Spectral Analysis 
        trace_to_show=double(trace_to_show);
        eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);        
    elseif (all_plotter<7 && all_plotter>=5)
        if segments_tracings(3,1+segments_already_staged)==segments_tracings(4,1+segments_already_staged)
           eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged):segments_tracings(6,1+segments_already_staged));']);
        else
           eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged)).(EOG_trace_',num2str(all_plotter-4),')(1,1:segments_tracings(6,1+segments_already_staged))];']);
        end
%         if all_plotter==6
%             trace_to_show=-trace_to_show;
%         end
        trace_to_show=double(trace_to_show);
        eval(['plot(handles.EOG_',num2str(all_plotter-4),',1:(segments_tracings(2,1+segments_already_staged)-segments_tracings(1,1+segments_already_staged)+1),BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP''),''b'');']);
        eval(['axes(handles.EOG_',num2str(all_plotter-4),');']);
            if all_plotter==5 %from here changed 100115
                hold on;
                plot(1:length(SA_5),(nanmean(SA_5)+num_stds_EOG*std(SA_5))*ones(1,length(SA_5)),'r','LineStyle','--');
                hold on;
                plot(1:length(SA_5),(nanmean(SA_5)-num_stds_EOG*std(SA_5))*ones(1,length(SA_5)),'r','LineStyle','--');
            elseif all_plotter==6
                hold on;
                plot(1:length(SA_6),(nanmean(SA_6)+num_stds_EOG*std(SA_6))*ones(1,length(SA_6)),'r','LineStyle','--');
                hold on;
                plot(1:length(SA_6),(nanmean(SA_6)-num_stds_EOG*std(SA_6))*ones(1,length(SA_6)),'r','LineStyle','--');                
            end %until here changed 100115 
            hold off;        
        ylabel('Voltage \muV')
        ylim([-400, 400])
        % save trace for Spectral Analysis 
        trace_to_show=double(trace_to_show);
        eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);                    
    elseif all_plotter==7
        if segments_tracings(3,1+segments_already_staged)==segments_tracings(4,1+segments_already_staged)
           eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged):segments_tracings(6,1+segments_already_staged));']);
        else
           eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged)).(EMG_trace_',num2str(all_plotter-6),')(1,1:segments_tracings(6,1+segments_already_staged))];']);
        end
        trace_to_show=double(trace_to_show);
        eval(['plot(handles.EMG_',num2str(all_plotter-6),',1:(segments_tracings(2,1+segments_already_staged)-segments_tracings(1,1+segments_already_staged)+1),BWfilter_designer(trace_to_show, fs, [10,100], [5,120], ''BP''),''r'');']);
        axes(handles.EMG_1);
        ylabel('Voltage \muV')
        ylim([-80 80])
        % save trace for Spectral Analysis 
        trace_to_show=double(trace_to_show);
        eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [10,100], [5,120], ''BP'');']);        
    end
    
numer=floor((segments_tracings(2,1+segments_already_staged)-segments_tracings(1,1+segments_already_staged)+1)/10);  % changed 100115
% set(gca,'Xtick',[1*numer 2*numer 3*numer 4*numer 5*numer 6*numer 7*numer 8*numer 9*numer 10*numer 11*numer 12*numer 13*numer 14*numer 15*numer 16*numer 17*numer 18*numer 19*numer 20*numer 21*numer 22*numer 23*numer 24*numer 25*numer 26*numer 27*numer 28*numer 29*numer 30*numer]) % changed 100115
set(gca,'Xtick',[1*numer 2*numer 3*numer 4*numer 5*numer 6*numer 7*numer 8*numer 9*numer 10*numer]) % changed 100115
% set(gca,'XTickLabel',char('1','2','3','4','5','6','7','8','9','10','11','12','13','14','15','16','17','18','19','20','21','22','23','24','25','26','27','28','29','30')); % changed 100115
set(gca,'XTickLabel',char('1','2','3','4','5','6','7','8','9','10')); % changed 100115
xlim([0 segments_tracings(2,1+segments_already_staged)-segments_tracings(1,1+segments_already_staged)])
grid on
end

% Introduction of polysom_ddmmyy variables into the GUI
handles.BeginFile=BeginFile;
handles.EndFile=EndFile;
handles.FilesPaths=FilesPaths; % inserted into polysom
handles.RecBegin=RecBegin;
handles.RecEnd=RecEnd;
handles.polysom_analog=polysom_analog;
handles.segments_already_staged=segments_already_staged; % This is the running var for the different segments being staged.
        Daily_Data_Cell{1,1}='Recording Date';
        Daily_Data_Cell{1,2}='Starting File';
        Daily_Data_Cell{1,3}='Ending File';
        Daily_Data_Cell{1,4}='# of Sample on Starting file';
        Daily_Data_Cell{1,5}='# of Sample on Ending file';
        Daily_Data_Cell{1,6}='Stage';
handles.Daily_Data_Cell=Daily_Data_Cell;
global polysomnography
polysomnography=struct;
polysomnography.Events.SWS.start=[]; polysomnography.Events.SWS.end=[];
polysomnography.Events.REMs.start=[]; polysomnography.Events.REMs.end=[];
polysomnography.Events.sleep_spindles.start=[]; polysomnography.Events.sleep_spindles.end=[];
polysomnography.Events.K_complexes.start=[]; polysomnography.Events.K_complexes.end=[];
polysomnography.Events.Repetitive_Muscle.start=[]; polysomnography.Events.Repetitive_Muscle.end=[];
polysomnography.Events.REM_without_atonia.start=[]; polysomnography.Events.REM_without_atonia.end=[];
polysomnography.Events.Arousal.start=[]; polysomnography.Events.Arousal.end=[];

handles.segments_tracings=segments_tracings;
handles.fs=fs;
handles.extra=extra;
handles.total_segment_num=total_segment_num;
handles.sec_2_stage=sec_2_stage;
handles.samples_per_segment=samples_per_segment;
handles.path_to_save=path_to_save;
guidata(hObject,handles);

function recording_date_Callback(hObject, eventdata, handles)
recording_date=get(hObject,'String');
handles.recording_date=recording_date;

guidata(hObject,handles);

function recording_date_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function stager_name_Callback(hObject, eventdata, handles)
stager_name_string=get(hObject,'String');
handles.stager_name_string=stager_name_string;

guidata(hObject,handles);

function stager_name_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function update_ddmmyy_Callback(hObject, eventdata, handles)
update_ddmmyy=get(hObject,'Value');
handles.update_ddmmyy=update_ddmmyy;

guidata(hObject,handles);

function monkey_name_Callback(hObject, eventdata, handles)
monkey_name_string=get(hObject,'String');
handles.monkey_name_string=monkey_name_string;

guidata(hObject,handles);

function monkey_name_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function path_to_ddmmyy_Callback(hObject, eventdata, handles)
path_to_ddmmyy_string=get(hObject,'String');
handles.path_to_ddmmyy_string=path_to_ddmmyy_string;

guidata(hObject,handles);

function path_to_ddmmyy_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function date_of_staging_Callback(hObject, eventdata, handles)
date_of_staging=get(hObject,'String');
handles.date_of_staging=date_of_staging;

guidata(hObject,handles);

function date_of_staging_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function figure1_ResizeFcn(hObject, eventdata, handles)

function Choose_Stage_mid_Callback(hObject, eventdata, handles)
list_mid=get(hObject,'String');
val_mid=get(hObject,'Value');
Choose_Stage_mid=list_mid{val_mid};
handles.Choose_Stage_mid=Choose_Stage_mid;
guidata(hObject,handles);

function Choose_Stage_mid_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end
set(hObject,'String',{'Choose Sleep Stage';'W';'N1';'N2';'N3';'R';'Unclassified'});

function Start_Sec_mid_Callback(hObject, eventdata, handles)
Start_Sec_mid=str2double(get(hObject,'String'));
handles.Start_Sec_mid=Start_Sec_mid;
guidata(hObject,handles);

function Start_Sec_mid_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function insert_midstage_Callback(hObject, eventdata, handles)
global insert_midstage;
insert_midstage=get(hObject,'Value');
guidata(hObject,handles);

function Update_Events_Callback(hObject, eventdata, handles)
global startevent;
global endevent;
global bit_resolution
global polysomnography;
fs=handles.fs;
raw_fs=44000;
segments_already_staged=handles.segments_already_staged;
segments_tracings=handles.segments_tracings;
extra=handles.extra;
global SWS;
global REMs
global K_complexes;
global Arousal;
global sleep_spindles;
global REM_without_atonia;
global Repetitive_Muscle;
recording_date=handles.recording_date;
name=[handles.path_to_ddmmyy_string,'\Data_',recording_date,'.mat'];
load(name,'datafiles');

Event_Time_Start=datafiles(segments_tracings(3,1+segments_already_staged-extra),1).start_sample-1+floor(segments_tracings(5,1+segments_already_staged-extra)+startevent/30*(segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1)-1)/fs*raw_fs;
Event_Time_End=datafiles(segments_tracings(3,1+segments_already_staged-extra),1).start_sample-1+floor(segments_tracings(5,1+segments_already_staged-extra)+endevent/30*(segments_tracings(2,1+segments_already_staged-extra)-segments_tracings(1,1+segments_already_staged-extra)+1)-1)/fs*raw_fs;

if SWS
    polysomnography.Events.SWS.start=[polysomnography.Events.SWS.start,Event_Time_Start];
    polysomnography.Events.SWS.end=[polysomnography.Events.SWS.end,Event_Time_End];
elseif REMs
    polysomnography.Events.REMs.start=[polysomnography.Events.REMs.start,Event_Time_Start];
    polysomnography.Events.REMs.end=[polysomnography.Events.REMs.end,Event_Time_End];
elseif sleep_spindles
    polysomnography.Events.sleep_spindles.start=[polysomnography.Events.sleep_spindles.start,Event_Time_Start];
    polysomnography.Events.sleep_spindles.end=[polysomnography.Events.sleep_spindles.end,Event_Time_End];
elseif K_complexes
    polysomnography.Events.K_complexes.start=[polysomnography.Events.K_complexes.start,Event_Time_Start];
    polysomnography.Events.K_complexes.end=[polysomnography.Events.K_complexes.end,Event_Time_End];
elseif Repetitive_Muscle
    polysomnography.Events.Repetitive_Muscle.start=[polysomnography.Events.Repetitive_Muscle.start,Event_Time_Start];
    polysomnography.Events.Repetitive_Muscle.end=[polysomnography.Events.Repetitive_Muscle.end,Event_Time_End];
elseif REM_without_atonia
    polysomnography.Events.REM_without_atonia.start=[polysomnography.Events.REM_without_atonia.start,Event_Time_Start];
    polysomnography.Events.REM_without_atonia.end=[polysomnography.Events.REM_without_atonia.end,Event_Time_End];
elseif Arousal
    polysomnography.Events.Arousal.start=[polysomnography.Events.Arousal.start,Event_Time_Start];
    polysomnography.Events.Arousal.end=[polysomnography.Events.Arousal.end,Event_Time_End];
end

% --- Executes on button press in Repetitive_Muscle.
function Repetitive_Muscle_Callback(hObject, eventdata, handles)
global Repetitive_Muscle
Repetitive_Muscle=get(hObject,'Value');

% hObject    handle to Repetitive_Muscle (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of Repetitive_Muscle


% --- Executes on button press in REM_without_atonia.
function REM_without_atonia_Callback(hObject, eventdata, handles)
global REM_without_atonia
REM_without_atonia=get(hObject,'Value');

% hObject    handle to REM_without_atonia (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of REM_without_atonia


% --- Executes on button press in arousal.
function arousal_Callback(hObject, eventdata, handles)
global Arousal
Arousal=get(hObject,'Value');

% hObject    handle to arousal (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of arousal


% --- Executes on button press in sleep_spindles.
function sleep_spindles_Callback(hObject, eventdata, handles)
global sleep_spindles
sleep_spindles=get(hObject,'Value');

% hObject    handle to sleep_spindles (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of sleep_spindles


% --- Executes on button press in K_complexes.
function K_complexes_Callback(hObject, eventdata, handles)
global K_complexes
K_complexes=get(hObject,'Value');
% hObject    handle to K_complexes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of K_complexes


% --- Executes on button press in SWS.
function SWS_Callback(hObject, eventdata, handles)
global SWS
SWS=get(hObject,'Value');
% hObject    handle to SWS (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of SWS


% --- Executes on button press in REMs.
function REMs_Callback(hObject, eventdata, handles)
global REMs
REMs=get(hObject,'Value');

% hObject    handle to REMs (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hint: get(hObject,'Value') returns toggle state of REMs


function startevent_Callback(hObject, eventdata, handles)
global startevent
startevent=str2double(get(hObject,'String'));
% hObject    handle to startevent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of startevent as text
%        str2double(get(hObject,'String')) returns contents of startevent as a double


% --- Executes during object creation, after setting all properties.
function startevent_CreateFcn(hObject, eventdata, handles)
% hObject    handle to startevent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function endevent_Callback(hObject, eventdata, handles)
global endevent
endevent=str2double(get(hObject,'String'));

% hObject    handle to endevent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Hints: get(hObject,'String') returns contents of endevent as text
%        str2double(get(hObject,'String')) returns contents of endevent as a double


% --- Executes during object creation, after setting all properties.
function endevent_CreateFcn(hObject, eventdata, handles)
% hObject    handle to endevent (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called

% Hint: edit controls usually have a white background on Windows.
%       See ISPC and COMPUTER.
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

function show_movie_Callback(hObject, eventdata, handles)
recording_date=handles.recording_date;
segments_tracings=handles.segments_tracings;
segments_already_staged=handles.segments_already_staged;
extra=handles.extra;
BeginFile=handles.BeginFile;
fs=handles.fs;
matfiles=handles.matfiles;
fps=50;
polysom_analog=handles.polysom_analog;
filenum=BeginFile+segments_tracings(3,segments_already_staged+1-extra)-1;
% load(polysom_analog(1,segments_tracings(3,segments_already_staged+1-extra)).name,'SF_NORPIX_FRAME','CLFP_017'); % I use the *last* file that's been updated because my code only updates Daily_Data_Cell after the actual staging, and the movie should be seen before the staging is done. % CHANGE BACK TO RAW
load(polysom_analog(1,segments_tracings(3,segments_already_staged+1-extra)).name,'SF_NORPIX_FRAME','CRAW_017');
% % I use the *last* file that's been updated because my code only updates
% Daily_Data_Cell after the actual staging, and the movie should be seen before the staging is done. 
% relevant_sample=(segments_tracings(5,segments_already_staged+1-extra))/fs*44000+CLFP_017.TimeBegin*44000-1; % CHANGE BACK TO RAW
relevant_sample=(segments_tracings(5,segments_already_staged+1-extra))/fs*44000+CRAW_017.TimeBegin*44000-1;
if  size(SF_NORPIX_FRAME,2)>10000
    [~,ia,~]=unique(SF_NORPIX_FRAME(2,:)); ia=ia(find(ia==1):end); SF_NORPIX_FRAME=SF_NORPIX_FRAME(:,ia); SF_NORPIX_FRAME(2,:)=SF_NORPIX_FRAME(2,:)-(SF_NORPIX_FRAME(2,1)-1); % ERASE!
end
% easy_access_folder=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',recording_date(1:2),'-',recording_date(3:4),'-',recording_date(5:6),'_file#',num2str(filenum,'%02u')];
% if ~exist(easy_access_folder)
%     easy_access_folder=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',recording_date(1:2),'-',num2str(str2num(recording_date(3:4))+1),'-',recording_date(5:6),'_file#',num2str(filenum,'%02u')]; % change of day
%     if ~exist(easy_access_folder)
%         easy_access_folder=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',num2str(str2num(recording_date(1:2))+1,'%02u'),'-01-',recording_date(5:6),'_file#',num2str(filenum,'%02u')]; % change of month
%         if ~exist(easy_access_folder)
%             easy_access_folder=['K:\Daenerys\',recording_date,'\VIDEO\31-145\01-01-',num2str(str2num(recording_date(5:6))+1,'%02u'),'_file#',num2str(filenum,'%02u')]; % change of year
%         end
%     end
% end
% next_easy_access_folder=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',recording_date(1:2),'-',recording_date(3:4),'-',recording_date(5:6),'_file#',num2str(filenum+1,'%02u')];
% next_easy_access_folder_alternative1=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',recording_date(1:2),'-',num2str(str2num(recording_date(3:4))+1,'%02u'),'-',recording_date(5:6),'_file#',num2str(filenum+1,'%02u')];  % change of day
% next_easy_access_folder_alternative2=['K:\Daenerys\',recording_date,'\VIDEO\31-145\',num2str(str2num(recording_date(1:2))+1,'%02u'),'-01-',recording_date(5:6),'_file#',num2str(filenum+1,'%02u')]; % change of month
% next_easy_access_folder_alternative3=['K:\Daenerys\',recording_date,'\VIDEO\31-145\01-01-',num2str(str2num(recording_date(5:6))+1,'%02u'),'_file#',num2str(filenum+1,'%02u')]; % change of year
% right one use always:
recording_date_1='122715';
recording_date_2='122815';
easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',recording_date_2(3:4),'-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % replace every 071415 with "recording_date" in this next bulk.
if ~exist(easy_access_folder)
    easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',num2str(str2num(recording_date_1(3:4))+1),'-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % change of day
    if ~exist(easy_access_folder)
        easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',num2str(str2num(recording_date_1(1:2))+1,'%02u'),'-01-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % change of month
        if ~exist(easy_access_folder)
            easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\01-01-',num2str(str2num(recording_date_1(5:6))+1,'%02u'),'_file#',num2str(filenum,'%02u')]; % change of year
        end
    end
end
next_easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\',recording_date_2,'\',recording_date_1(1:2),'-',recording_date_2(3:4),'-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')];
next_easy_access_folder_alternative1=['I:\Daenerys\',recording_date_1,'\VIDEO\',recording_date_2,'\',recording_date_1(1:2),'-',num2str(str2num(recording_date_1(3:4))+1,'%02u'),'-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')];  % change of day
next_easy_access_folder_alternative2=['I:\Daenerys\',recording_date_1,'\VIDEO\',recording_date_2,'\',num2str(str2num(recording_date_1(1:2))+1,'%02u'),'-01-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')]; % change of month
next_easy_access_folder_alternative3=['I:\Daenerys\',recording_date_1,'\VIDEO\',recording_date_2,'\01-01-',num2str(str2num(recording_date_1(5:6))+1,'%02u'),'_file#',num2str(filenum+1,'%02u')]; % change of year

%only 062915
% recording_date_1='062115';
% recording_date_2='062215';
% easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',recording_date_2(3:4),'-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % replace every 071415 with "recording_date" in this next bulk.
% if ~exist(easy_access_folder)
%     easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',num2str(str2num(recording_date_1(3:4))+1),'-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % change of day
%     if ~exist(easy_access_folder)
%         easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',num2str(str2num(recording_date_1(1:2))+1,'%02u'),'-01-',recording_date_1(5:6),'_file#',num2str(filenum,'%02u')]; % change of month
%         if ~exist(easy_access_folder)
%             easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\01-01-',num2str(str2num(recording_date_1(5:6))+1,'%02u'),'_file#',num2str(filenum,'%02u')]; % change of year
%         end
%     end
% end
% next_easy_access_folder=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',recording_date_2(3:4),'-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')];
% next_easy_access_folder_alternative1=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',recording_date_1(1:2),'-',num2str(str2num(recording_date_1(3:4))+1,'%02u'),'-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')];  % change of day
% next_easy_access_folder_alternative2=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\',num2str(str2num(recording_date_1(1:2))+1,'%02u'),'-01-',recording_date_1(5:6),'_file#',num2str(filenum+1,'%02u')]; % change of month
% next_easy_access_folder_alternative3=['I:\Daenerys\',recording_date_1,'\VIDEO\1\',recording_date_2,'\01-01-',num2str(str2num(recording_date_1(5:6))+1,'%02u'),'_file#',num2str(filenum+1,'%02u')]; % change of year

[~,ind_of_frame]=min(abs(SF_NORPIX_FRAME(1,:)-relevant_sample));
list = ls(easy_access_folder);
list=list(3:end,:);
if strcmp(list(end,:),'FileList.mat               ')
    list=list(1:end-1,:);
end
total_frame_num=fps*10;
F=[];
if length(list)-ind_of_frame>=total_frame_num
    for ii=1:total_frame_num
        F(:,:,ii) = imread([easy_access_folder,'\',list(ind_of_frame+ii-1,:)],'jpg');
    end
%     disp(['file num is: ',num2str(filenum)]);
%     disp(['index of frame is: ',num2str(ind_of_frame)]);
%     disp(['relevant sample: ',num2str(segments_tracings(5,segments_already_staged+1-extra))]);    
    %implay(F, fps); % 101316 Alex's Advice
    h = implay(F,fps); % 101316 Alex's Advice
    set(findall (0, 'tag', 'spcui_scope_framework'), 'position', [150 150 700 550]); % 101316 Alex's Advice
    h.Visual.ColorMap.UserRange = 1; % 101316 Alex's Advice
    play(h.DataSource.Controls);  % 101316 Alex's Advice
else
    for ii=1:length(list)-ind_of_frame+1
        F(:,:,ii) = imread([easy_access_folder,'\',list(ind_of_frame+ii-1,:)],'jpg');
    end
    if isdir(next_easy_access_folder)
        list = ls(next_easy_access_folder); list=list(3:end,:);
        for jj=1:total_frame_num-ii
            F(:,:,ii+jj)=imread([next_easy_access_folder,'\',list(jj,:)],'jpg');
        end
    elseif isdir(next_easy_access_folder_alternative1)
        list = ls(next_easy_access_folder_alternative1);  list=list(3:end,:);           
        for jj=1:total_frame_num-ii
            F(:,:,ii+jj)=imread([next_easy_access_folder_alternative1,'\',list(jj,:)],'jpg');
        end      
    elseif isdir(next_easy_access_folder_alternative2)
        list = ls(next_easy_access_folder_alternative2); list=list(3:end,:);           
        for jj=1:total_frame_num-ii
            F(:,:,ii+jj)=imread([next_easy_access_folder_alternative2,'\',list(jj,:)],'jpg');
        end
    elseif isdir(next_easy_access_folder_alternative3)
        list = ls(next_easy_access_folder_alternative3); list=list(3:end,:);           
        for jj=1:total_frame_num-ii
            F(:,:,ii+jj)=imread([next_easy_access_folder_alternative3,'\',list(jj,:)],'jpg');
        end
    end
%     disp(['file num is: ',num2str(filenum)]); % all commented on 101716
%     disp(['index of frame is: ',num2str(ind_of_frame)]);
%     disp(['beginning is: ',num2str(ii)]);
%     disp(['remainder is: ',num2str(jj)]);
%     disp(['relevant sample: ',num2str(segments_tracings(5,segments_already_staged+1-extra))]);
    %implay(F, fps); % 101316 Alex's Advice
    h = implay(F,fps); % 101316 Alex's Advice
    set(findall (0, 'tag', 'spcui_scope_framework'), 'position', [150 150 700 550]); % 101316 Alex's Advice
    h.Visual.ColorMap.UserRange = 1; % 101316 Alex's Advice
    play(h.DataSource.Controls);  % 101316 Alex's Advice
    
end

function pushbutton7_Callback(hObject, eventdata, handles)
global SA_1 SA_2 SA_3 SA_4 SA_5 SA_6 SA_7
sec_2_stage=30;
fs=handles.fs;
spectral_edge_frequency=0.50;
%perform Spectral Analysis of the now presented traces

figure('units','normalized','outerposition',[0 0 1 1]);
subplot(2,16,1:4);
freq_vec=[0.5:0.5:40];
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_1,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
ylabel('Frequency (Hz)')
title('EEG 1','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)

l=legend(['sef-',num2str(spectral_edge_frequency)],'Location','NorthEast');
y=get(l,'Position');
y(1)=0.075; y(2)=0.87;
set(l,'Position',y)
subplot(2,16,5:8);
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_2,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
title('EEG 2','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
ylabel('')
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)

subplot(2,16,9:12);
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_3,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
title('EEG 3','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
ylabel('')
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)


subplot(2,16,13:16);
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_4,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
title('EEG 4','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
ylabel('')
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)


subplot(2,16,18:21);
freq_vec=[0.5:0.5:40];
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_5,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
ylabel('Frequency (Hz)')
title('EOG 1','FontWeight','Bold','Color','b','FontSize',14);
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)


subplot(2,16,22:25);
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_6,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
title('EOG 1','FontWeight','Bold','Color','b','FontSize',14);
ylabel('')
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)


subplot(2,16,28:31);
freq_vec=[1:2:200];
% compute the Spectral Edge Frequency-x
[s,w,t]=spectrogram(SA_7,1024,512,freq_vec,fs);
s=abs(s);
sef=nan(1,size(s,2));
for ii=1:size(s,2)
    s(:,ii)=s(:,ii)/sum(s(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
    edge_freq=spectral_edge_frequency*sum(s(:,ii));
    sef(1,ii)=freq_vec(find(cumsum(s(:,ii))>edge_freq*sum(s(:,ii)),1,'first'));
end
imagesc(t,w,10*log10(s)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
plot(t,sef,'k','LineWidth',2);
xlabel('Time (sec)')
title('EMG','FontWeight','Bold','Color','r','FontSize',14);
ylabel('Frequency (Hz)')
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/4;
x(1)=x(1)+x(3)/4;
x1(3)=x1(3)+5*x(3);
set(c,'Position',x)
set(gca,'position',x1)

suptitle('Spectral Analysis for Polysomnographic Measures')

function Auto_Staging_Callback(hObject, eventdata, handles)
global num_clusters
num_clusters=3;
global polygon_mat
polygon_mat=cell(3,3);
global temp_xy temp_xz temp_yz
global SA_1 SA_2 SA_3 SA_4 SA_5 SA_6 SA_7
BeginFile=handles.BeginFile;
EndFile=handles.EndFile;
FilesPaths=handles.FilesPaths;
RecBegin=handles.RecBegin;
RecEnd=handles.RecEnd;
recording_date=handles.recording_date;
monkey_name=handles.monkey_name;
global polysomnography
global bit_resolution
global measure_box
global num_stds_EOG
global SEF
global SEF_1
global iii
bit_resolution=1.9073; % uV per value
global color_matrix
color_matrix=[0,0,1;0,1,0;1,0,0];
global human_clustering
global post_hoc_mm
post_hoc_mm=0;
global wake_nrem_spect % creates a spectrogram for Wake and NREM, so that a better ratio measure could be made, based on the actual wake and NREM in the data
wake_nrem_spect=0;
wake_spect=[];
nrem_spect=[];
polysom_analog=handles.polysom_analog;
Daily_Data_Cell=handles.Daily_Data_Cell;
path_to_save=handles.path_to_save;
segments_already_staged=handles.segments_already_staged;
segments_tracings=handles.segments_tracings;
fs=handles.fs;
sec_2_stage=handles.sec_2_stage;
samples_per_segment=handles.samples_per_segment;
total_segment_num=handles.total_segment_num;
extra=handles.extra;
human_clustering=measure_box(:,7);
% for the remaining segments, i.e. those that are not classified by the
% user, continue filling measure box, only with no user classification.

% compute measures and classification and save them
% EEG power ratio
while segments_already_staged<total_segment_num
        
Daily_Data_Cell{segments_already_staged+2,1}=recording_date; % Only to check
Daily_Data_Cell{segments_already_staged+2,2}=BeginFile+segments_tracings(3,1+segments_already_staged-extra)-1;
Daily_Data_Cell{segments_already_staged+2,3}=BeginFile+segments_tracings(4,1+segments_already_staged-extra)-1;
Daily_Data_Cell{segments_already_staged+2,4}=segments_tracings(5,1+segments_already_staged-extra);
Daily_Data_Cell{segments_already_staged+2,5}=segments_tracings(6,1+segments_already_staged-extra);
    
    rms_thresh_EEG=10; % threshold rms for "non-real" EEG signals.
    num_of_EEG_segs=0;
    window_length=min(2^12,length(SA_1));
    ratio_1=0;
    ratio_2=0;
    ratio_3=0;
    SEF_50=0;
    SEF_result=zeros(10,1); %%^^%%    
    name=[handles.path_to_ddmmyy_string,'\Data_',recording_date,'.mat'];
    load(name,'LickBlinkData');
    if rms(SA_1)>=rms_thresh_EEG
        [pxx1,f1] = pwelch(SA_1,hann(window_length),window_length/2,window_length,fs);
        % power ratio 1
        lows1_1=nanmean(pxx1(f1>=0.1 & f1<=5));
        highs1_1=nanmean(pxx1(f1>=20 & f1<=30)); % GPe
%         lows1_1=nanmean(pxx1(f1>=0 & f1<=15)); % GPe
%         lows1_1=nanmean(pxx1((f1>=0.1 & f1<=5) | (f1>=10 & f1<=17))); % striatum
%         highs1_1=nanmean(pxx1(f1>=17 & f1<=25)); % striatum
        % power ratio 2
        lows1_2=nanmean(pxx1(f1>=0.05 & f1<=15));
        highs1_2=nanmean(pxx1(f1>=20 & f1<=35));        
        ratio_1=ratio_1+highs1_1/lows1_1;        
        ratio_2=ratio_2+highs1_2/lows1_2;        
        % SEF_50
        pxx1=pxx1/sum(pxx1); % normalize power spectra
        edge_freq=SEF*sum(pxx1);
        SEF_50=SEF_50+f1(find(cumsum(pxx1)>=edge_freq,1,'first'));
        ratio_3=ratio_3+f1(find(cumsum(pxx1)>=edge_freq,1,'first'))/lows1_1; % added Aviv 101316	
        e1=SEF_1(1,1)*sum(pxx1); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f1(find(cumsum(pxx1)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx1); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f1(find(cumsum(pxx1)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx1); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f1(find(cumsum(pxx1)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx1); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f1(find(cumsum(pxx1)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx1); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f1(find(cumsum(pxx1)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx1); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f1(find(cumsum(pxx1)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx1); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f1(find(cumsum(pxx1)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx1); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f1(find(cumsum(pxx1)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx1); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f1(find(cumsum(pxx1)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx1); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f1(find(cumsum(pxx1)>=e10,1,'first')); %%^^%%                
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_2)>=rms_thresh_EEG
        [pxx2,f2] = pwelch(SA_2,hann(window_length),window_length/2,window_length,fs);
        lows2_1=nanmean(pxx2(f2>=0.1 & f2<=5)); % GPe
        highs2_1=nanmean(pxx2(f2>=20 & f2<=30)); % GPe
%         lows2_1=nanmean(pxx2((f2>=0.1 & f2<=5) | (f2>=10 & f2<=17))); % striatum
%         highs2_1=nanmean(pxx2(f2>=17 & f2<=25)); % striatum
        lows2_2=nanmean(pxx2(f2>=0.05 & f2<=15));
        highs2_2=nanmean(pxx2(f2>=20 & f2<=35));        
        ratio_1=ratio_1+highs2_1/lows2_1;
        ratio_2=ratio_2+highs2_2/lows2_2;
        % SEF_50
        pxx2=pxx2/sum(pxx2); % normalize power spectra
        edge_freq=SEF*sum(pxx2);
        SEF_50=SEF_50+f2(find(cumsum(pxx2)>=edge_freq,1,'first'));
        ratio_3=ratio_3+f2(find(cumsum(pxx2)>=edge_freq,1,'first'))/lows2_1; % added Aviv 101316	        
        e1=SEF_1(1,1)*sum(pxx2); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f2(find(cumsum(pxx2)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx2); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f2(find(cumsum(pxx2)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx2); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f2(find(cumsum(pxx2)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx2); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f2(find(cumsum(pxx2)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx2); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f2(find(cumsum(pxx2)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx2); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f2(find(cumsum(pxx2)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx2); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f2(find(cumsum(pxx2)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx2); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f2(find(cumsum(pxx2)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx2); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f2(find(cumsum(pxx2)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx2); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f2(find(cumsum(pxx2)>=e10,1,'first')); %%^^%%                   
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_3)>=rms_thresh_EEG
        [pxx3,f3] = pwelch(SA_3,hann(window_length),window_length/2,window_length,fs);
        lows3_1=nanmean(pxx3(f3>=0.1 & f3<=5)); % GPe
        highs3_1=nanmean(pxx3(f3>=20 & f3<=30)); % GPe
%         lows3_1=nanmean(pxx3((f3>=0.1 & f3<=5) | (f3>=10 & f3<=17))); % striatum
%         highs3_1=nanmean(pxx3(f3>=17 & f3<=25));   % striatum     
        lows3_2=nanmean(pxx3(f3>=0.05 & f3<=15));
        highs3_2=nanmean(pxx3(f3>=20 & f3<=35));
        ratio_1=ratio_1+highs3_1/lows3_1;        
        ratio_2=ratio_2+highs3_2/lows3_2;
        % SEF_50
        pxx3=pxx3/sum(pxx3); % normalize power spectra
        edge_freq=SEF*sum(pxx3);
        SEF_50=SEF_50+f3(find(cumsum(pxx3)>=edge_freq,1,'first'));
        ratio_3=ratio_3+f3(find(cumsum(pxx3)>=edge_freq,1,'first'))/lows3_1; % added Aviv 101316	        
        e1=SEF_1(1,1)*sum(pxx3); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f3(find(cumsum(pxx3)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx3); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f3(find(cumsum(pxx3)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx3); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f3(find(cumsum(pxx3)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx3); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f3(find(cumsum(pxx3)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx3); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f3(find(cumsum(pxx3)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx3); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f3(find(cumsum(pxx3)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx3); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f3(find(cumsum(pxx3)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx3); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f3(find(cumsum(pxx3)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx3); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f3(find(cumsum(pxx3)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx3); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f3(find(cumsum(pxx3)>=e10,1,'first')); %%^^%%        
        num_of_EEG_segs=num_of_EEG_segs+1;
    end
    if rms(SA_4)>=rms_thresh_EEG
       [pxx4,f4] = pwelch(SA_4,hann(window_length),window_length/2,window_length,fs);
        lows4_1=nanmean(pxx4(f4>=0.1 & f4<=5)); % GPe
        highs4_1=nanmean(pxx4(f4>=20 & f4<=30)); % GPe
%         lows4_1=nanmean(pxx4((f4>=0.1 & f4<=5) | (f4>=10 & f4<=17))); % striatum
%         highs4_1=nanmean(pxx4(f4>=17 & f4<=25)); % striatum        
        lows4_2=nanmean(pxx4(f4>=0.05 & f4<=15));
        highs4_2=nanmean(pxx4(f4>=20 & f4<=35));
        ratio_1=ratio_1+highs4_1/lows4_1;        
        ratio_2=ratio_2+highs4_2/lows4_2;        
        % SEF_50
        pxx4=pxx4/sum(pxx4); % normalize power spectra
        edge_freq=SEF*sum(pxx4);
        SEF_50=SEF_50+f4(find(cumsum(pxx4)>=edge_freq,1,'first'));
        ratio_3=ratio_3+f4(find(cumsum(pxx4)>=edge_freq,1,'first'))/lows4_1; % added Aviv 101316	
        e1=SEF_1(1,1)*sum(pxx4); %%^^%%
        SEF_result(1,1)=SEF_result(1,1)+f4(find(cumsum(pxx4)>=e1,1,'first')); %%^^%%
        e2=SEF_1(1,2)*sum(pxx4); %%^^%%
        SEF_result(2,1)=SEF_result(2,1)+f4(find(cumsum(pxx4)>=e2,1,'first')); %%^^%%
        e3=SEF_1(1,3)*sum(pxx4); %%^^%%
        SEF_result(3,1)=SEF_result(3,1)+f4(find(cumsum(pxx4)>=e3,1,'first')); %%^^%%
        e4=SEF_1(1,4)*sum(pxx4); %%^^%%
        SEF_result(4,1)=SEF_result(4,1)+f4(find(cumsum(pxx4)>=e4,1,'first')); %%^^%%
        e5=SEF_1(1,5)*sum(pxx4); %%^^%%
        SEF_result(5,1)=SEF_result(5,1)+f4(find(cumsum(pxx4)>=e5,1,'first')); %%^^%%
        e6=SEF_1(1,6)*sum(pxx4); %%^^%%
        SEF_result(6,1)=SEF_result(6,1)+f4(find(cumsum(pxx4)>=e6,1,'first')); %%^^%%
        e7=SEF_1(1,7)*sum(pxx4); %%^^%%
        SEF_result(7,1)=SEF_result(7,1)+f4(find(cumsum(pxx4)>=e7,1,'first')); %%^^%%
        e8=SEF_1(1,8)*sum(pxx4); %%^^%%
        SEF_result(8,1)=SEF_result(8,1)+f4(find(cumsum(pxx4)>=e8,1,'first')); %%^^%%
        e9=SEF_1(1,9)*sum(pxx4); %%^^%%
        SEF_result(9,1)=SEF_result(9,1)+f4(find(cumsum(pxx4)>=e9,1,'first')); %%^^%%        
        e10=SEF_1(1,10)*sum(pxx4); %%^^%%
        SEF_result(10,1)=SEF_result(10,1)+f4(find(cumsum(pxx4)>=e10,1,'first')); %%^^%%        
        num_of_EEG_segs=num_of_EEG_segs+1;
    end    
    ratio_1=ratio_1/num_of_EEG_segs;
    ratio_2=ratio_2/num_of_EEG_segs;
    ratio_3=ratio_3/num_of_EEG_segs;
    SEF_50=SEF_50/num_of_EEG_segs;
    SEF_result=SEF_result/num_of_EEG_segs; %%^^%%    
    if segments_already_staged==0
        measure_box(1,1)=log10(ratio_1);
    else
        measure_box(size(measure_box,1)+1,1)=log10(ratio_1);
    end

    measure_box(end,2)=log10(ratio_2);
    measure_box(end,3)=SEF_50;
    measure_box(end,8)=SEF_result(1,1); %%^^%%
    measure_box(end,9)=SEF_result(2,1); %%^^%%
    measure_box(end,10)=SEF_result(3,1); %%^^%%
    measure_box(end,11)=SEF_result(4,1); %%^^%%
    measure_box(end,12)=SEF_result(5,1); %%^^%%
    measure_box(end,13)=SEF_result(6,1); %%^^%%
    measure_box(end,14)=SEF_result(7,1); %%^^%%
    measure_box(end,15)=SEF_result(8,1); %%^^%%
    measure_box(end,16)=SEF_result(9,1); %%^^%%
    measure_box(end,17)=SEF_result(10,1); %%^^%%    
    measure_box(end,18)=log10(ratio_3); % Added Aviv 101316

    % EOG events number
    num_of_EOG_segs=0;
    rms_thresh_EOG=10;
    crossnum=0;
    if rms(SA_5)>=rms_thresh_EOG
        all_above_points1=find(SA_5>=nanmean(SA_5)+std(SA_5)*num_stds_EOG);
        all_below_points1=find(SA_5<=nanmean(SA_5)-std(SA_5)*num_stds_EOG);
        crossnum=crossnum+length(find(diff(all_above_points1)>1))+1; % +1 for the ever first one that is missed by this method
        crossnum=crossnum+length(find(diff(all_below_points1)>1))+1;
        num_of_EOG_segs=num_of_EOG_segs+1;
    end
    if rms(SA_6)>=rms_thresh_EOG
        all_above_points2=find(SA_6>=nanmean(SA_6)+std(SA_6)*num_stds_EOG);
        all_below_points2=find(SA_6<=nanmean(SA_6)-std(SA_6)*num_stds_EOG);
        crossnum=crossnum+length(find(diff(all_above_points2)>1))+1; % +1 for the ever first one that is missed by this method
        crossnum=crossnum+length(find(diff(all_below_points2)>1))+1;
        num_of_EOG_segs=num_of_EOG_segs+1;
    end    

    measure_box(end,4)=crossnum/num_of_EOG_segs;

    % EMG root mean square
    measure_box(end,5)=rms(SA_7);

    % fraction of eye_open (=1)
    load(name,'datafiles');
    file_num=Daily_Data_Cell{segments_already_staged+2,2}-Daily_Data_Cell{2,2}+1;
    segment_start=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,4}-1)*(44000/fs)+1;
    segment_end=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,5})*(44000/fs);
    if segment_end<segment_start
        file_num=file_num+1;
        segment_end=datafiles(file_num,1).start_sample-1+(Daily_Data_Cell{segments_already_staged+2,5})*(44000/fs);
    end        
    relevant_blink=LickBlinkData.EyeState(find(LickBlinkData.SNR_TimeStamp>=segment_start & LickBlinkData.SNR_TimeStamp<=segment_end));
    if ~isnan((length(relevant_blink)-sum(relevant_blink))/length(relevant_blink))
        measure_box(end,6)=(length(relevant_blink)-sum(relevant_blink))/length(relevant_blink);
    else
        measure_box(end,6)=1.5;
    end

    if wake_nrem_spect
        min_rms_thresh=100;
        min_eye_open_thresh=0.2;
        nrem_segments=[308:311,325:336,398:400,449:453,495:502,554:555];
        freq_vec=[0.5:0.5:40];
        if ismember(size(measure_box,1),nrem_segments)
            s=[];
            if rms(SA_1)>=rms_thresh_EEG
                [s1,w,t]=spectrogram(SA_1,1024,512,freq_vec,fs);
                s1=abs(s1);
                for ii=1:size(s1,2)
                    s1(:,ii)=s1(:,ii)/sum(s1(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end
                s=cat(3,s,s1);
            end
             if rms(SA_2)>=rms_thresh_EEG
                [s2,w,t]=spectrogram(SA_2,1024,512,freq_vec,fs);
                s2=abs(s2);
                for ii=1:size(s2,2)
                    s2(:,ii)=s2(:,ii)/sum(s2(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end      
                s=cat(3,s,s2);
             end
             if rms(SA_3)>=rms_thresh_EEG
                [s3,w,t]=spectrogram(SA_3,1024,512,freq_vec,fs);
                s3=abs(s3);
                for ii=1:size(s3,2)
                    s3(:,ii)=s3(:,ii)/sum(s3(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end
                s=cat(3,s,s3);
             end
             if rms(SA_4)>=rms_thresh_EEG
                [s4,w,t]=spectrogram(SA_4,1024,512,freq_vec,fs);
                s4=abs(s4);
                for ii=1:size(s4,2)
                    s4(:,ii)=s4(:,ii)/sum(s4(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end    
                s=cat(3,s,s4);
             end
             s=nanmean(s,3);
             nrem_spect=cat(3,nrem_spect,s);
        elseif (measure_box(end,5)>=min_rms_thresh & measure_box(end,6)>=min_eye_open_thresh)
            s=[];
            if rms(SA_1)>=rms_thresh_EEG
                [s1,w,t]=spectrogram(SA_1,1024,512,freq_vec,fs);
                s1=abs(s1);
                for ii=1:size(s1,2)
                    s1(:,ii)=s1(:,ii)/sum(s1(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end
                s=cat(3,s,s1);
            end
             if rms(SA_2)>=rms_thresh_EEG
                [s2,w,t]=spectrogram(SA_2,1024,512,freq_vec,fs);
                s2=abs(s2);
                for ii=1:size(s2,2)
                    s2(:,ii)=s2(:,ii)/sum(s2(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end      
                s=cat(3,s,s2);
             end
             if rms(SA_3)>=rms_thresh_EEG
                [s3,w,t]=spectrogram(SA_3,1024,512,freq_vec,fs);
                s3=abs(s3);
                for ii=1:size(s3,2)
                    s3(:,ii)=s3(:,ii)/sum(s3(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end
                s=cat(3,s,s3);
             end
             if rms(SA_4)>=rms_thresh_EEG
                [s4,w,t]=spectrogram(SA_4,1024,512,freq_vec,fs);
                s4=abs(s4);
                for ii=1:size(s4,2)
                    s4(:,ii)=s4(:,ii)/sum(s4(:,ii)); % normalize each time stamp according to the sum of powers in that time stamp
                end    
                s=cat(3,s,s4);
             end
             s=nanmean(s,3);
             wake_spect=cat(3,wake_spect,s);  
        end
    end
 
    
    disp(['finished working on segment ',num2str(segments_already_staged+1),'/',num2str(total_segment_num)])

    if segments_already_staged<total_segment_num-1
        segments_already_staged=segments_already_staged+1;
        handles.segments_already_staged=segments_already_staged;
        
        % TODO FIND WHICH IS THE BETTER EEG TRACE
        EEG_trace_1=handles.EEG_trace_1;
        EEG_trace_2=handles.EEG_trace_2;
        EEG_trace_3=handles.EEG_trace_3;
        EEG_trace_4=handles.EEG_trace_4;

        EOG_trace_1=handles.EOG_trace_1;
        EOG_trace_2=handles.EOG_trace_2;

        EMG_trace_1=handles.EMG_trace_1;

        for all_plotter=1:7
            if all_plotter<5 % EEG data traces
                if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EEG_trace_',num2str(all_plotter),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EEG_trace_',num2str(all_plotter),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
                end
                trace_to_show=double(trace_to_show);
                eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);
            elseif (all_plotter<7 && all_plotter>=5)
                if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
                    eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EOG_trace_',num2str(all_plotter-4),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EOG_trace_',num2str(all_plotter-4),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
                end
                trace_to_show=double(trace_to_show);
                eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [0.1,35], [0.05,40], ''BP'');']);
            elseif all_plotter==7
                if segments_tracings(3,1+segments_already_staged-extra)==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):segments_tracings(6,1+segments_already_staged-extra));']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+1==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']);
                elseif segments_tracings(3,1+segments_already_staged-extra)+2==segments_tracings(4,1+segments_already_staged-extra)
                   eval(['trace_to_show=[polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,segments_tracings(5,1+segments_already_staged-extra):end),polysom_analog(1,segments_tracings(3,1+segments_already_staged-extra)+1).(EMG_trace_',num2str(all_plotter-6),')(1,1:end),polysom_analog(1,segments_tracings(4,1+segments_already_staged-extra)).(EMG_trace_',num2str(all_plotter-6),')(1,1:segments_tracings(6,1+segments_already_staged-extra))];']); 
                end
                trace_to_show=double(trace_to_show);
                eval(['SA_',num2str(all_plotter),'=BWfilter_designer(trace_to_show, fs, [10,100], [5,120], ''BP'');']);
            end
        end
    else
        segments_already_staged=segments_already_staged+1; % to avoid staying in the while loop forever.
        handles.segments_already_staged=segments_already_staged;        
    end
end

if wake_nrem_spect
% plot the wake and NREM matrices
nrem_spect=nanmean(nrem_spect,3);
wake_spect=nanmean(wake_spect,3);
figure;
subplot(1,2,1);
imagesc(t,w,10*log10(wake_spect)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
xlabel('Time (sec)')
ylabel('Frequency (Hz)')
title('WAKE','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)
subplot(1,2,2);
imagesc(t,w,10*log10(nrem_spect)); %plot the log spectrum
set(gca,'YDir', 'normal'); % flip the Y Axis so lower frequencies are at the bottom
hold all;
xlabel('Time (sec)')
ylabel('Frequency (Hz)')
title('NREM','FontWeight','Bold','Color',[0 0.5 0],'FontSize',14);
c=colorbar;
x1=get(gca,'position');
x=get(c,'Position');
x(3)=x(3)/2;
x(1)=x(1)+x(3)/2;
x1(3)=x1(3)+4*x(3);
set(c,'Position',x)
set(gca,'position',x1)
end


% cluster the all-night segments using the k-means algorithm
all_night_classification_matrix=measure_box(:,[1,5,6]); % 18 for ratio_3, 1 for ratio_1, 3 for SEF_50
[all_night_machine_classification,Centers,sumd,D] = kmeans(all_night_classification_matrix,num_clusters);

% next, A measure of the separation of each point in the cluster from the
% points in the other cluster and in the same, is calculated. Separation is
% calculated in very much the same way like Mati Joshua calculated his 
% isolation score 2009 manuscript. The calculation is explained below.
A=all_night_classification_matrix;
B=permute(repmat(A,[1,1,size(A,1)]),[1,3,2])-permute(repmat(A',[1,1,size(A',2)]),[3,2,1]); % find the euclidean distance between each dot in the 3-dim space and all other dots. If, for example, we want to find the distances between all points in the 1-dim vector [a,b,c,d], we create two 2-dim matrices - [a,a,a,a;b,b,b,b;c,c,c,c;d,d,d,d;] and [a,b,c,d;a,b,c,d;a,b,c,d;a,b,c,d;]. This is done by transposing A, and replicating A and A'. Then the matrices are subtracted to get [0,a-b,a-c,a-d;b-a,0,b-c,b-d;c-a,c-b,0,c-d;d-a,d-b,d-c,0], such that in the new matrix, element Aij is the difference between i and j, again, in the 1-dim space. To apply the same calculation for 3 dimensions, I use repmat in 3-d, and thenI permute the matrices to allow the right subtraction.
C=sqrt(sum(B.*B,3)); % the euclidean difference is squared, summed across the third dimension to add the x,y and z distances, and then squared.
% for each point, calculate an separation score with gamma=10
gamma=10; separation_threshold=0.9;
sep=exp(-C)*gamma; 
separation_score=nan(1,length(all_night_machine_classification));
for ii=1:length(separation_score)
    separation_score(1,ii)=sum(sep(setdiff(find(all_night_machine_classification==all_night_machine_classification(ii,1)),ii),ii))/(sum(sep(:,ii))-sep(ii,ii));
end
bold=find(separation_score>=separation_threshold);
% % next, A measure of the distance from each point in the cluster to the 
% % other centroids, is measured. The points which are in the X% of points 
% % farthest from the centroid of cluster 1 *and* in the X% of points 
% % farthest from the centroid of cluster 2 are filled. The others are not.
% bold=[];
% thresh=0.9;
% for cluster=1:3
%     this_cluster=D(find(all_night_machine_classification==cluster),~ismember([1:num_clusters],cluster));
%     [~,indx]=sort(this_cluster(:,1));
%     indx=flipud(indx);
%     a=indx(1:floor(thresh*length(indx))); % finds the thresh farthest segments from first of other clusters' centroid
%     [~,indx]=sort(this_cluster(:,2));
%     indx=flipud(indx);
%     b=indx(1:floor(thresh*length(indx))); % finds the thresh farthest segments from first of other clusters' centroid
%     bold=[bold,intersect(a,b)]; % The bolded segments are those that are exported for all clusters together.
% end
% Present Clustering Results and allow user to choose satisfactory segments
% This is done in the following way: For every sleep stage (NREM, REM and
% WAKE), a new figure is presented to the user, with three projections of
% the 3-parameter space used for clustering. In each such projection, the
% user is asked to delineate a polygon in which all the relevant N/R/W
% segments are situated. Only a segment which is included in the polygon in
% all dimensions isconsidered a genuine N/R/W point.

%draw the 3d plot to help the user orient in the 3-d
    figure;
%   filled_points=intersect(find(all_night_machine_classification==1),bold);
%   empty_points=setdiff(find(all_night_machine_classification==1),bold);
%   scatter3(all_night_classification_matrix(filled_points,1),all_night_classification_matrix(filled_points,2),all_night_classification_matrix(filled_points,3),'o','b','filled');        hold all;
%   scatter3(all_night_classification_matrix(empty_points,1),all_night_classification_matrix(empty_points,2),all_night_classification_matrix(empty_points,3),'*','b');
%   hold all;
%   filled_points=intersect(find(all_night_machine_classification==2),bold);
%   empty_points=setdiff(find(all_night_machine_classification==2),bold);
%   scatter3(all_night_classification_matrix(filled_points,1),all_night_classification_matrix(filled_points,2),all_night_classification_matrix(filled_points,3),'o','g','filled');
%   hold all;
%   scatter3(all_night_classification_matrix(empty_points,1),all_night_classification_matrix(empty_points,2),all_night_classification_matrix(empty_points,3),'*','g');
%   hold all;    
%   filled_points=intersect(find(all_night_machine_classification==3),bold);
%   empty_points=setdiff(find(all_night_machine_classification==3),bold);
%   scatter3(all_night_classification_matrix(filled_points,1),all_night_classification_matrix(filled_points,2),all_night_classification_matrix(filled_points,3),'o','r','filled');
%   hold all;
%   scatter3(all_night_classification_matrix(empty_points,1),all_night_classification_matrix(empty_points,2),all_night_classification_matrix(empty_points,3),'*','r');
    
    % modify separation scores to have a clear distribution
%     smallest=sort(separation_score);
    smallest=min(separation_score);
%     smallest=smallest(round(0.05*length(smallest))); % find the 95%-smallest separation score
    desired_colorband_width=0.8; % the separation will always be 0.8 wide, from 0.2 to 1.
    separation_score=separation_score-smallest;
    separation_score=0.8*separation_score/(1-smallest); % convert the (1-smallest) existing gap to a 0.8 gap between the ~smallest and the largest values
    separation_score=separation_score+0.2; % to move back to having the highest value at 1.
    
    zz=1:length(all_night_machine_classification);
    cluster_1=find(all_night_machine_classification==1); 
    relevant_separation_score=separation_score(cluster_1);
    rgb_colormap=[ones(length(cluster_1),1),ones(length(cluster_1),1)-relevant_separation_score',ones(length(cluster_1),1)-relevant_separation_score'];
    all_bigger=find(rgb_colormap>1);
    for yy=1:length(all_bigger)
        rgb_colormap(all_bigger)=1;
    end
    all_smaller=find(rgb_colormap<0);
    for yy=1:length(all_smaller)
        rgb_colormap(all_smaller)=0;
    end    
    scatter3(all_night_classification_matrix(cluster_1,1),all_night_classification_matrix(cluster_1,2),all_night_classification_matrix(cluster_1,3),50,zz(1:length(cluster_1)),'filled'); colormap(rgb_colormap);
    hold on;
    scatter3(Centers(1,1),Centers(1,2),Centers(1,3),100,'d','k','filled')    
    hold on;
    cluster_2=find(all_night_machine_classification==2); 
    relevant_separation_score=separation_score(cluster_2);
    rgb_colormap=[rgb_colormap;ones(length(cluster_2),1)-relevant_separation_score',ones(length(cluster_2),1),ones(length(cluster_2),1)-relevant_separation_score'];
    all_bigger=find(rgb_colormap>1);
    for yy=1:length(all_bigger)
        rgb_colormap(all_bigger)=1;
    end
    all_smaller=find(rgb_colormap<0);
    for yy=1:length(all_smaller)
        rgb_colormap(all_smaller)=0;
    end        
    scatter3(all_night_classification_matrix(cluster_2,1),all_night_classification_matrix(cluster_2,2),all_night_classification_matrix(cluster_2,3),50,zz(length(cluster_1)+1:length(cluster_1)+length(cluster_2)),'filled'); colormap(rgb_colormap);
    hold on;
    scatter3(Centers(2,1),Centers(2,2),Centers(2,3),100,'d','k','filled')
    hold on;
    cluster_3=find(all_night_machine_classification==3); 
    relevant_separation_score=separation_score(cluster_3);
    rgb_colormap=[rgb_colormap;ones(length(cluster_3),1)-relevant_separation_score',ones(length(cluster_3),1)-relevant_separation_score',ones(length(cluster_3),1)];
    all_bigger=find(rgb_colormap>1);
    for yy=1:length(all_bigger)
        rgb_colormap(all_bigger)=1;
    end
    all_smaller=find(rgb_colormap<0);
    for yy=1:length(all_smaller)
        rgb_colormap(all_smaller)=0;
    end       
    
    scatter3(all_night_classification_matrix(cluster_3,1),all_night_classification_matrix(cluster_3,2),all_night_classification_matrix(cluster_3,3),50,zz(length(cluster_1)+length(cluster_2)+1:length(cluster_1)+length(cluster_2)+length(cluster_3)),'filled'); colormap(rgb_colormap);
    hold on;
    scatter3(Centers(3,1),Centers(3,2),Centers(3,3),100,'d','k','filled')    
    hold on;
    xlabel('High to Low EEG ratio','FontSize',12,'FontWeight','Bold');
    ylabel('EMG rms','FontSize',12,'FontWeight','Bold');
    zlabel('Eye Open fraction','FontSize',12,'FontWeight','Bold');
    title('K-means Clustering Results - NREM, REM, W','FontSize',14,'FontWeight','Bold')
    
for iii=1:3 % (R -> N -> W)
    % draw three projections
    projection_matrix=[1,2;2,3;1,3]; % xy, yz, xz
    for kk=1:3
    switch iii
        case 1
            fh=figure('name','mark_R');
        case 2
            fh=figure('name','mark_N');
        case 3
            fh=figure('name','mark_W');    
    end
    set(fh,'units','normalized','position',[0 0 1 1],'outerposition',[0 0 1 1]);
    if iii==1
        suptitle('In the three following projections, Delineate the Polygon Containing all certain REM points')
    elseif iii==2
        suptitle('In the three following projections, Delineate the Polygon Containing all certain NREM points')
    elseif iii==3
        suptitle('In the three following projections, Delineate the Polygon Containing all certain WAKE points')
    end        
        subplot(3,3,[1,2,4,5,7,8])
        points_so_far=0;
        for jj=1:num_clusters
%             filled_points=intersect(find(all_night_machine_classification==jj),bold);
%             empty_points=setdiff(find(all_night_machine_classification==jj),bold);
%             if jj==1
%                 scatter(all_night_classification_matrix(filled_points,projection_matrix(kk,1)),all_night_classification_matrix(filled_points,projection_matrix(kk,2)),'o','b','filled');
%             elseif jj==2
%                 scatter(all_night_classification_matrix(filled_points,projection_matrix(kk,1)),all_night_classification_matrix(filled_points,projection_matrix(kk,2)),'o','g','filled');
%             elseif jj==3
%                 scatter(all_night_classification_matrix(filled_points,projection_matrix(kk,1)),all_night_classification_matrix(filled_points,projection_matrix(kk,2)),'o','r','filled');
%             end
%             hold all;
%             if jj==1
%                 scatter(all_night_classification_matrix(empty_points,projection_matrix(kk,1)),all_night_classification_matrix(empty_points,projection_matrix(kk,2)),'*','b');
%             elseif jj==2
%                 scatter(all_night_classification_matrix(empty_points,projection_matrix(kk,1)),all_night_classification_matrix(empty_points,projection_matrix(kk,2)),'*','g');
%             elseif jj==3
%                 scatter(all_night_classification_matrix(empty_points,projection_matrix(kk,1)),all_night_classification_matrix(empty_points,projection_matrix(kk,2)),'*','r');
%             end
            cluster_points=find(all_night_machine_classification==jj); 
             if jj==1
                scatter(all_night_classification_matrix(cluster_points,projection_matrix(kk,1)),all_night_classification_matrix(cluster_points,projection_matrix(kk,2)),50,zz(points_so_far+1:points_so_far+length(cluster_points)),'filled'); colormap(rgb_colormap);
                hold on;
                scatter(Centers(jj,projection_matrix(kk,1)),Centers(jj,projection_matrix(kk,2)),100,'d','k','filled')
                hold on;
                points_so_far=points_so_far+length(cluster_points);
            elseif jj==2
                scatter(all_night_classification_matrix(cluster_points,projection_matrix(kk,1)),all_night_classification_matrix(cluster_points,projection_matrix(kk,2)),50,zz(points_so_far+1:points_so_far+length(cluster_points)),'filled'); colormap(rgb_colormap);
                hold on;
                scatter(Centers(jj,projection_matrix(kk,1)),Centers(jj,projection_matrix(kk,2)),100,'d','k','filled')
                hold on;
                points_so_far=points_so_far+length(cluster_points);                
            elseif jj==3
                scatter(all_night_classification_matrix(cluster_points,projection_matrix(kk,1)),all_night_classification_matrix(cluster_points,projection_matrix(kk,2)),50,zz(points_so_far+1:points_so_far+length(cluster_points)),'filled'); colormap(rgb_colormap);
                hold on;
                scatter(Centers(jj,projection_matrix(kk,1)),Centers(jj,projection_matrix(kk,2)),100,'d','k','filled')
                hold on;
                points_so_far=points_so_far+length(cluster_points);                
            end
            hold all;  
            xlim([min(all_night_classification_matrix(:,projection_matrix(kk,1)))-0.2*(max(all_night_classification_matrix(:,projection_matrix(kk,1)))-min(all_night_classification_matrix(:,projection_matrix(kk,1)))),max(all_night_classification_matrix(:,projection_matrix(kk,1)))+0.2*(max(all_night_classification_matrix(:,projection_matrix(kk,1)))-min(all_night_classification_matrix(:,projection_matrix(kk,1))))]);
            ylim([min(all_night_classification_matrix(:,projection_matrix(kk,2)))-0.2*(max(all_night_classification_matrix(:,projection_matrix(kk,2)))-min(all_night_classification_matrix(:,projection_matrix(kk,2)))),max(all_night_classification_matrix(:,projection_matrix(kk,2)))+0.2*(max(all_night_classification_matrix(:,projection_matrix(kk,2)))-min(all_night_classification_matrix(:,projection_matrix(kk,2))))]);
        end
        switch kk
            case 1
                xlabel('EEG Power Ratio')
                ylabel('EMG Root Mean Square')
            case 2
                xlabel('EMG Root Mean Square')
                ylabel('Eye-Open Ratio')
            case 3
                xlabel('EEG Power Ratio')
                ylabel('Eye-Open Ratio')
        end
        if kk==1
            uicontrol(fh,'Style','pushbutton','String','Mark XY Segments','Position',[1400 800 100 40],'Callback','global temp_xy;temp_xy=impoly(gca,[]);')
            uicontrol(fh,'Style','pushbutton','String','Save and exit','Position',[1400 100 100 40],'Callback','global temp_xy iii polygon_mat; polygon_mat{iii,1}=getPosition(temp_xy); close;')
            waitfor(fh); 
        elseif kk==2
            uicontrol(fh,'Style','pushbutton','String','Mark YZ Segments','Position',[1400 500 100 40],'Callback','global temp_yz;temp_yz=impoly(gca,[]);')
            uicontrol(fh,'Style','pushbutton','String','Save and exit','Position',[1400 100 100 40],'Callback','global temp_yz iii polygon_mat; polygon_mat{iii,2}=getPosition(temp_yz); close;')
            waitfor(fh);             
        elseif kk==3
            uicontrol(fh,'Style','pushbutton','String','Mark XZ Segments','Position',[1400 200 100 40],'Callback','global temp_xz;temp_xz=impoly(gca,[]);')
            uicontrol(fh,'Style','pushbutton','String','Save and exit','Position',[1400 100 100 40],'Callback','global temp_xz iii polygon_mat; polygon_mat{iii,3}=getPosition(temp_xz); close;')
            waitfor(fh);             
        end
    end
end

% apply the polygon border data supplied by the user and create the final
% set of segments to be regarded R, N or W
all_segments=1:size(measure_box,1); % simply a list of indexes, 1 to num of segments
final_R_indices=intersect(all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,2),polygon_mat{1,1}(:,1),polygon_mat{1,1}(:,2))),intersect(all_segments(inpolygon(all_night_classification_matrix(:,2),all_night_classification_matrix(:,3),polygon_mat{1,2}(:,1),polygon_mat{1,2}(:,2))),all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,3),polygon_mat{1,3}(:,1),polygon_mat{1,3}(:,2)))));
final_N_indices=intersect(all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,2),polygon_mat{2,1}(:,1),polygon_mat{2,1}(:,2))),intersect(all_segments(inpolygon(all_night_classification_matrix(:,2),all_night_classification_matrix(:,3),polygon_mat{2,2}(:,1),polygon_mat{2,2}(:,2))),all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,3),polygon_mat{2,3}(:,1),polygon_mat{2,3}(:,2)))));
final_W_indices=intersect(all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,2),polygon_mat{3,1}(:,1),polygon_mat{3,1}(:,2))),intersect(all_segments(inpolygon(all_night_classification_matrix(:,2),all_night_classification_matrix(:,3),polygon_mat{3,2}(:,1),polygon_mat{3,2}(:,2))),all_segments(inpolygon(all_night_classification_matrix(:,1),all_night_classification_matrix(:,3),polygon_mat{3,3}(:,1),polygon_mat{3,3}(:,2)))));
final_U_indices=setdiff(setdiff(setdiff([1:size(measure_box,1)],final_R_indices),final_N_indices),final_W_indices);

for ii=1:length(final_R_indices)
    Daily_Data_Cell{final_R_indices(ii)+1,6}='R';
end
for ii=1:length(final_N_indices)
    Daily_Data_Cell{final_N_indices(ii)+1,6}='N2';
end
for ii=1:length(final_W_indices)
    Daily_Data_Cell{final_W_indices(ii)+1,6}='W';
end
for ii=1:length(final_U_indices)
    Daily_Data_Cell{final_U_indices(ii)+1,6}='U'; % unspecified segments
end

% compute post-hoc man-machine agreement

human_W=find(human_clustering==1);
human_N=find(human_clustering==2);
human_R=find(human_clustering==3);
human_U=find(human_clustering==0);
post_hoc_mm=(sum(ismember(human_W,final_W_indices))+sum(ismember(human_N,final_N_indices))+sum(ismember(human_R,final_R_indices))+sum(ismember(human_U,final_U_indices)))/length(human_clustering);

excel_file_name=['Staging for date_',recording_date]; 
xlswrite([handles.path_to_save,'\Staging Folder for Date ',recording_date,'\',excel_file_name,'.xls'],Daily_Data_Cell); % save to excel
save([handles.path_to_save,'\Staging Folder for Date ',recording_date,'\',excel_file_name], 'Daily_Data_Cell'); % save to cell
save([handles.path_to_save,'\polysom_',recording_date,'.mat'],'BeginFile','EndFile','FilesPaths','RecBegin','RecEnd','polysom_analog','recording_date','monkey_name','path_to_save','Daily_Data_Cell','segments_already_staged','-v7.3');

errordlg(['All segments are Staged. Post-hoc Man-Machine Agreement is ',num2str(post_hoc_mm)],'Error'); 
if handles.update_ddmmyy==1
    name=[handles.path_to_ddmmyy_string,'\Data_',recording_date,'.mat'];
    load(name);
    num_of_stages=6; stage_matrix={'W','N1','N2','N3','R','U'};
    for stage=1:num_of_stages
        relevant_indices=find(strcmp(stage_matrix(1,stage),Daily_Data_Cell(2:end,6)));
        ind=find(diff(relevant_indices)==1)+1;
        openers=relevant_indices;
        openers(ind)=[];
        ind=find(diff(relevant_indices)==1);
        closers=relevant_indices;
        closers(ind)=[];
        stage_openers=nan(1,length(openers));
        for pp=1:length(openers)
            stage_openers(1,pp)=datafiles(Daily_Data_Cell{openers(pp)+1,2}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{openers(pp)+1,4}/fs*(generalInfo.SAMPLING*1000));
        end
        stage_closers=nan(1,length(closers));
        for ss=1:length(closers)
            stage_closers(1,ss)=datafiles(Daily_Data_Cell{closers(ss)+1,3}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{closers(ss)+1,5}/fs*(generalInfo.SAMPLING*1000));
        end
        eval(['polysomnography.stages.',stage_matrix{1,stage},'.begin_sample=stage_openers;']);
        eval(['polysomnography.stages.',stage_matrix{1,stage},'.end_sample=stage_closers;']);
        disp(stage)
    end

     polysomnography.Staging_date=handles.date_of_staging;
     polysomnography.Stager=handles.stager_name_string;
     
     % create *all* segments openers and closers
     all_begin=Daily_Data_Cell{2:end,4}; all_end=Daily_Data_Cell{2:end,5};
     for aa=1:size(Daily_Data_Cell,1)-1
         all_begin(aa,1)=datafiles(Daily_Data_Cell{aa+1,2}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{aa+1,4}/fs*(generalInfo.SAMPLING*1000));
         all_end(aa,1)=datafiles(Daily_Data_Cell{aa+1,3}-BeginFile+1,1).start_sample-1+floor(Daily_Data_Cell{aa+1,5}/fs*(generalInfo.SAMPLING*1000));
     end
     %correction that is due to the "dead samples" that are created because
     %of /fs*(generalInfo.SAMPLING*1000)) - because
     %fs<generalInfo.SAMPLING*1000, a group of undefined "dead samples"
     %occur. added 231014 MIGHT DO PROBLEMS IN neuron profile
     all_begin_new=all_begin;
     all_end_new=all_end;
     all_begin_new(1,1)=datafiles(Daily_Data_Cell{1+1,2}-BeginFile+1,1).start_sample;
     all_end_new(end,1)=datafiles(Daily_Data_Cell{end,2}-BeginFile+1,1).end_sample;
     
     for runner=1:length(all_begin_new)-1
         midpoint=round((all_end_new(runner,1)+all_begin_new(runner+1,1))/2);
         all_end_new(runner,1)=midpoint;
         all_begin_new(runner+1,1)=midpoint+1;
     end
     
     A=polysomnography.stages.W.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.W.begin_sample=A; clear A;
     A=polysomnography.stages.W.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.W.begin_sample=A; clear A;
     A=polysomnography.stages.W.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.W.end_sample=A; clear A;
     A=polysomnography.stages.W.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.W.end_sample=A; clear A;
     
     A=polysomnography.stages.N1.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N1.begin_sample=A; clear A;
     A=polysomnography.stages.N1.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N1.begin_sample=A; clear A;
     A=polysomnography.stages.N1.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N1.end_sample=A; clear A;
     A=polysomnography.stages.N1.end_sample; A=changem(A,all_end_new,all_end);  polysomnography.stages.N1.end_sample=A; clear A;
     
     A=polysomnography.stages.N2.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N2.begin_sample=A; clear A;
     A=polysomnography.stages.N2.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N2.begin_sample=A; clear A;
     A=polysomnography.stages.N2.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N2.end_sample=A; clear A;
     A=polysomnography.stages.N2.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N2.end_sample=A; clear A;
     
     A=polysomnography.stages.N3.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N3.begin_sample=A; clear A;
     A=polysomnography.stages.N3.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.N3.begin_sample=A; clear A;
     A=polysomnography.stages.N3.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.N3.end_sample=A; clear A;
     A=polysomnography.stages.N3.end_sample; A=changem(A,all_end_new,all_end);polysomnography.stages.N3.end_sample=A; clear A;
     
     A=polysomnography.stages.R.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.R.begin_sample=A; clear A;
     A=polysomnography.stages.R.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.R.begin_sample=A; clear A;
     A=polysomnography.stages.R.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.R.end_sample=A; clear A;
     A=polysomnography.stages.R.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.R.end_sample=A; clear A;
     
     A=polysomnography.stages.U.begin_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.U.begin_sample=A; clear A;
     A=polysomnography.stages.U.begin_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.U.begin_sample=A; clear A;
     A=polysomnography.stages.U.end_sample; A=changem(A,all_begin_new,all_begin); polysomnography.stages.U.end_sample=A; clear A;
     A=polysomnography.stages.U.end_sample; A=changem(A,all_end_new,all_end); polysomnography.stages.U.end_sample=A; clear A;     
     
     all_begin=all_begin_new; clear all_begin_new
     all_end=all_end_new; clear all_end_new
     polysomnography.all_beginnings=all_begin;
     polysomnography.all_ends=all_end;

    % compute general statistics for the session
    % Total Sleep Time
    polysomnography.Stats.Total_Sleep_Time=(sum(polysomnography.stages.N1.end_sample-polysomnography.stages.N1.begin_sample)+sum(polysomnography.stages.N2.end_sample-polysomnography.stages.N2.begin_sample)+sum(polysomnography.stages.N3.end_sample-polysomnography.stages.N3.begin_sample)+sum(polysomnography.stages.R.end_sample-polysomnography.stages.R.begin_sample))/44000;
    % Sleep Efficacy
    polysomnography.Stats.Sleep_Efficacy=polysomnography.Stats.Total_Sleep_Time/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000);
    % Sleep Fragmentation
    % I'm still not sure as to how to compute this
    % REM Latency
    if (length(polysomnography.stages.R.begin_sample)>=1 && length(polysomnography.stages.N1.begin_sample)>=1 && length(polysomnography.stages.N3.begin_sample)>=1 && length(polysomnography.stages.N2.begin_sample)>=1)
        polysomnography.Stats.REM_Latency=polysomnography.stages.R.begin_sample(1,1)-min([polysomnography.stages.W.begin_sample(1,1),polysomnography.stages.N1.begin_sample(1,1),polysomnography.stages.N2.begin_sample(1,1),polysomnography.stages.N3.begin_sample(1,1),polysomnography.stages.R.begin_sample(1,1)])/44000;
    end
    % Total REM (precentage)
    polysomnography.Stats.Total_REM=(sum(polysomnography.stages.R.end_sample-polysomnography.stages.R.begin_sample)/44000)/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000);
    % Total SWS
    polysomnography.Stats.Total_SWS=(sum(polysomnography.Events.SWS.end-polysomnography.Events.SWS.start)/(polysomnography.Stats.Total_Sleep_Time+sum(polysomnography.stages.W.end_sample-polysomnography.stages.W.begin_sample)/44000));
    % And see some more in my Email
     
    save (name,'electrode','generalInfo','datafiles','polysomnography','LickBlinkData');
end

handles.extra=extra;
handles.total_segment_num=total_segment_num;
guidata(hObject,handles);

function Man_Machine_Agreement_Callback(hObject, eventdata, handles)
global measure_box
global num_clusters
global color_matrix
global human_clustering
num_clusters=3;

% compute and scatter plot the k-means clustering for the user-staged
% segments, then report a auto-human agreement measure.
classification_matrix=measure_box(:,[1,5,6]); % relevant measures are ratio1, EMG rms and eye-open ratio % 18 for ratio_3
%normalize classification matrix to have a normal distribution around 0
%with std=1
% for jj=1:size(classification_matrix,2)
%     classification_matrix(:,jj)=(classification_matrix(:,jj)-nanmean(classification_matrix(:,jj)))/std(classification_matrix(:,jj));
% end
[machine_classification,centers] = kmeans(classification_matrix,num_clusters);
% the agreement ratio is computed like this: For every *classified* segment
% in the human clustered data, the classifications of the k-means algorithm
% is inspected. The most prevalent classification in the k-means indexing
% is regarded the classification corresponding to the initial human one,
% and the others are regarded as mistakes. The total number of mistakes in
% all clusters is computed.
human_clustering=measure_box(:,7);
num_of_mistakes=0;
figure;
for ii=1:num_clusters
    cluster_to_inspect=machine_classification(find(human_clustering==ii));
    most_frequent_machine_classification=mode(cluster_to_inspect);
    if ii==1
        classifications_so_far=most_frequent_machine_classification;
    else
        a=tabulate(cluster_to_inspect); 
        [b,ind]=sort(a(:,2));
        appearance_frequencies=flipud(ind); % these 3 lines find the most frequent classifications by order of frequency.
        c=double(ismember(appearance_frequencies,classifications_so_far));
        most_frequent_machine_classification=appearance_frequencies(find(c==0,1,'first')); % the last 2 lines find the first classification which was not the "right" classification of a different cluster.
        classifications_so_far=[classifications_so_far,most_frequent_machine_classification];
    end
    num_of_mistakes=num_of_mistakes+length(find(cluster_to_inspect~=most_frequent_machine_classification));
    color_matrix=[0,0,1;0,1,0;1,0,0];
    AA=classification_matrix(find(human_clustering==ii & machine_classification==most_frequent_machine_classification),:);
    scatter3(AA(:,1),AA(:,2),AA(:,3),'o','MarkerFaceColor',color_matrix(ii,:),'MarkerEdgeColor',color_matrix(ii,:));
    hold all;
    scatter3(centers(most_frequent_machine_classification,1),centers(most_frequent_machine_classification,2),centers(most_frequent_machine_classification,3),'d','MarkerFaceColor','k','MarkerEdgeColor','k');
    hold all;
    BB=classification_matrix(find(human_clustering==ii & machine_classification~=most_frequent_machine_classification),:);
    scatter3(BB(:,1),BB(:,2),BB(:,3),'*','MarkerFaceColor',color_matrix(ii,:),'MarkerEdgeColor',color_matrix(ii,:));
    hold all;    
end
error_ratio=num_of_mistakes/length(human_clustering);
title(['Initial Clustering results. o/* = Agreement/disagreement. Error ratio is ',num2str(error_ratio)])
xlabel('High to Low EEG ratio');
ylabel('EMG rms');
zlabel('Eye Open fraction');
