function FrequencyDetect(varargin)
%FreqDetect GUI which comes as part of Auto_Detect.
%User runs seizure detection algorithm based on changes of the frequecy
%components of an LFP trace.


global Freq_channel reference_channel KA_injection_time sample_rate
global t1 t2 Freq_channel_time Filter_Lower_Bound Filter_Upper_Bound



%initialising saving variables
Freq_channel = []; KA_injection_time = []; t1 = []; t2 = [];
sample_rate = []; Freq_channel_time =[];


inputVar = varargin{1};

Freq_channel = inputVar{1}(:);
KA_injection_time = inputVar{2};
t1 = inputVar{3}(1);
t2 = inputVar{3}(2);
sample_rate = inputVar{4};
reference_channel = inputVar{5};


%%%%%%%%%%%%%%%%%%%%%%% Default Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

WindowSize_Default = 1;
WindowOverlap_Default = 0.5;
Seizure_Duration_Default = 10;
InterSeizure_duration_Default = 5;
Standard_Deviation_Threshold_Default = 1;
Filter_lower_limit_Default = 49;
Filter_Upper_limit_Default = 51;

%Delta = 1, Theta = 2, Alpha = 3, Beta = 4, Gamma = 5, Broad = 6, Manual =
%7
Default_wave_band = 1; %alpha band


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%name of wave bands to filter for analysis
WaveBandName = {'Delta (0-4) Hz','Theta (4-8) Hz','Alpha (8-12) Hz','Beta (12-30) Hz','Gamma (30-100) Hz',...
     'Broad (0.4-100) Hz','Manual'};

%gets the screen size of user monitor
screensize=get(0,'ScreenSize');

% Generating the Freq Detect figure window
FreqFigure = figure('Units','characters','MenuBar','none','Name',...
    'Frequency Detect','NumberTitle','off','Position',[10 0 230 70],'ReSize',...
    'on','Visible','on');

%Generating panels 
FreqInputPanel = uipanel('Parent',FreqFigure,'Units','normalized',...
    'Position',[0.001 0.51 0.20 0.48],'Title', 'Parameters', 'FontWeight', 'bold');
FreqAllPlotsPanel =uipanel('Parent',FreqFigure,'Units','normalized',...
    'Position',[0.205 0.01 0.79 0.97],'Title', 'Frequency Detect Results', 'FontWeight', 'bold');


% Loading icons for the plot manipulation buttons
undock = load('undock.mat');



%Create axis to display LFP 
FreqChannelAxes = axes('Parent',FreqAllPlotsPanel,'Units','normalized',...
    'OuterPosition',[0 0.6 1 0.4]);
%creating plot manipulation toggle buttons
FreqChannelUnDock = uicontrol('Parent',FreqAllPlotsPanel,'Units',...
    'normalized','Callback',{@UnDock,FreqChannelAxes},'CData',...
    undock.cdata,'Position',[0.93 0.69375 0.03 0.05],'TooltipString',...
    'Click to Undock Plot into a new Figure');

%Create axis to display seizure detection graph
FreqDetectAxes = axes('Parent',FreqAllPlotsPanel,'Units','normalized',...
    'OuterPosition',[0 0.4 1 0.2]);
%creating plot manipulation toggle buttons
FreqDetectUnDock = uicontrol('Parent',FreqAllPlotsPanel,'Units',...
    'normalized','Callback',{@UnDock,FreqDetectAxes},'CData',...
    undock.cdata,'Position',[0.93 0.39375 0.02 0.03],'TooltipString',...
    'Click to Undock Plot into a new Figure');


%create axis to display spectrogram
FreqSpectrogramAxes = axes('Parent',FreqAllPlotsPanel,'Units','normalized',...
    'OuterPosition',[0 0 1 0.4]);
%creating plot manipulation toggle buttons
FreqSpectrogramUnDock = uicontrol('Parent',FreqAllPlotsPanel,'Units',...
    'normalized','Callback',{@UnDock,FreqSpectrogramAxes},'CData',...
    undock.cdata,'Position',[0.93 0.0375 0.03 0.05],'TooltipString',...
    'Click to Undock Plot into a new Figure');

% Locating all plot manipulation buttons
allToggleBtns = findall(gcf,'Style','togglebutton');


%setting string for input titles
FreqInputText = {'Sample Rate:','','KA Injection Time:','', 'Filter:',...
    '','Sliding Window Size:','','Window Overlap:','',...
    'Min Inter-Seizure Period:','', 'Min Seizure Duration:','',...
    'Frequency Band for Analysis:','','','',...
    'Or, manually set band:','','Standard Deviation Threshold:'};


%setting user input title text
FreqInputTitleText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','left',...
    'Position',[0.05 0.04 0.75 0.92],'String',FreqInputText,'Style',...
    'text');


% Setting the information section values - {''} gets filled in later func
FreqInputVarText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','right',...
    'Position',[0.65 0.78 0.3 0.18],'String',{''},'Style',...
    'text');

% Setting butterworth Filter  lower bound text field
 FreqFilterLowerBoundTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.455 0.785 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter upper limit of Butterworth filter');


% Setting butterworth Filter upper bound text field
 FreqFilterUpperBoundTextEdit= uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.785 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter lower limit of Butterworth filter');

% Setting Filter  text field units
 FreqFilterUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.77 0.07 0.05],'String','Hz',...
    'Style','text');

% Setting Sliding Window text field
FreqSlidingWindowTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.72 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter the size of the sliding window');

%  Setting Sliding Window text field units
 FreqSlidingWindowUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.71 0.07 0.05],'String','s',...
    'Style','text');

% Setting text field
FreqWindowOverlapTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.65 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter the overlap of the sliding window');

%  Setting WindowOverlap text field units
 FreqWindowOverlapUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.64 0.07 0.05],'String','s',...
    'Style','text');

% Setting minimum InterSeizure period text field
FreqMinInterSeizTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.58 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter the minimum duration of time between seizure events');

%  Setting minimum Inter Seizure text field units
 FreqMinInterSeizUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.57 0.07 0.05],'String','s',...
    'Style','text');

% Setting minimum SeizDuration text field
FreqMinSeizDurationTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.51 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter the minimum seizure duration');

%  Setting minimum SeizDuration text field units
FreqMinSeizDurationUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.50 0.07 0.05],'String','s',...
    'Style','text');


% Band Selection drop-down list
FreqBandSelectDrop = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Enable',...
    'on','Position',[0.55 0.38 0.37 0.05],'String',WaveBandName,...
    'Callback',@DropListManual,'Style','popupmenu','BackgroundColor','w','Value',...
    Default_wave_band);


% Setting manual frequency bound text field
FreqManualUpBoundTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.76 0.31 0.16 0.05],'String',[],...
    'Callback',@ManualBoundsCallback,'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter upper frequency bounds for analysis');

%  Setting ManualBounds field units
FreqManualLowBoundTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.55 0.31 0.16 0.05],'String',[],...
    'Callback',@ManualBoundsCallback,'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter lower frequency bounds for analysis');

%  Setting ManualBounds field units
FreqManualBoundsUnitsText = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','HorizontalAlignment','Left','Position',[0.923 0.29 0.07 0.05],'String','Hz',...
    'Style','text');

% Setting Standard Deviation Threshold text field
FreqStandardDevThreshTextEdit = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.715 0.24 0.21 0.05],'String',[],...
    'Style','edit','BackgroundColor','w','TooltipString',...
    'Enter the standard deviation threshold above the chosen frequency band');

% algorithm 'Run' button
FreqRunDetectButton = uicontrol('Parent',FreqInputPanel,'Units',...
    'normalized','Position',[0.2 0.01 0.6 0.09],'Callback',@FreqRunDetectButtonCallback,'Enable','on',...
    'String','Run','TooltipString','Runs algorithm with the set parameters');

% Export data to excel button
FreqExportButton = uicontrol('Parent',FreqFigure,'Units',...
    'normalized','Position',[0.04 0.04 0.13 0.034],'Callback',@ExcelExport,'Enable','on',...
    'String','Export Data to Excel','FontWeight','bold','TooltipString','Exports Seizure information to excel');



% Locating all objects on the Main Window
allObjects = findall(FreqFigure);



%Plot segment chosen from previous window
FreqChannelPlot();


set(FreqInputVarText,'String',{[num2str(sample_rate),' Hz'],'',...
    [num2str(KA_injection_time),'  s']});


% Make the Freq Detect window visible now that initialisation has been completed
set(FreqFigure,'Visible','on');





function  ButtonDown(hObject,eventdata)
        
        % Clear the curent string in the parent text field if there is a
        % button down function
        set(hObject,'Enable','on','String',[])
        uicontrol(hObject);
        
end

 function UnDock(hObject,eventdata,axes_handles)
        
        % Locating all plot manipulation buttons except for the currently
        % selected one
        otherToggleBtns = allToggleBtns;
        otherToggleBtns(allToggleBtns==hObject) = [];

        % Set the properties for all other toggle buttons
        for jj = 1:length(otherToggleBtns)
            
            % Toggle off all other plot manipulation buttons
            set(otherToggleBtns(jj),'Value',0)
            
        end
        
        % Create a new figure window on which to export the parent plot
        SeparatePlot = figure('ReSize','on');
        
        % Copy all the objects from the parent plot axes to this new figure
        % window
        copyobj(axes_handles,SeparatePlot);
        
        % Resize the plot axes to fill the new figure window
        set(gca,'OuterPosition',[0 0 1 1])
        
    end

    function DropListManual(hObject,eventdata)
        
        val = get(hObject, 'Value')
        
        if val > 6 
            
            set(FreqManualUpBoundTextEdit,'String',[],'ButtondownFcn',@ButtonDown,...
            'Enable','on'); 
            set(FreqManualLowBoundTextEdit,'String',[],'ButtondownFcn',@ButtonDown,...
            'Enable','on'); 
        end
        
    end

  function ChannelNotchFilter(Filter_Lower_Bound,Filter_Upper_Bound)
            %finds length of segment chosen in Main Window
        L=length(Freq_channel); 
        
               
        % Next power of 2 from length of y
        NFFT = 2^nextpow2(L); 
        freq = sample_rate/2*linspace(0,1,NFFT/2+1);
        nyquist=(sample_rate/2)+1;
             
        w1=Filter_Lower_Bound/nyquist;
        w2=Filter_Upper_Bound/nyquist;

        %Butterworth filter - stops frequencies between bounds w1 and w2
        %defined by the user.Filters mains frequency noise by default between 48 and 52Hz
        %mains frequency noise = 50Hz
        [b,a]=butter(5,[w1,w2],'stop');
        
        %applying Buttorworth filter and redefining global variable
        Freq_channel=filter(b,a,Freq_channel);
        
  end

    function ManualBoundsCallback(hObject,eventdata)
        
       %IF manual bounds is selected then the manual option is enabled       
        if get(FreqBandSelectDrop,'Value') <= 6;
            
            set(FreqManualBoundsTextEdit,'String',...
            Standard_Deviation_Threshold_Default,'ButtondownFcn',@ButtonDown,...
            'Enable','off');
       
        end
        
    end
    

    function FreqChannelPlot
        
        %Set and dispay default parameters in Paramter panel
        set(FreqFilterLowerBoundTextEdit,'String',Filter_lower_limit_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqFilterUpperBoundTextEdit,'String',Filter_Upper_limit_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqSlidingWindowTextEdit,'String',WindowSize_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqWindowOverlapTextEdit,'String',WindowOverlap_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqMinSeizDurationTextEdit,'String',Seizure_Duration_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqMinInterSeizTextEdit,'String',InterSeizure_duration_Default,...
            'ButtondownFcn',@ButtonDown,'Enable','on');
        set(FreqManualLowBoundTextEdit,'String',[],'ButtondownFcn',@ButtonDown,...
            'Enable','off'); 
        set(FreqManualUpBoundTextEdit,'String',[],'ButtondownFcn',@ButtonDown,...
            'Enable','off'); 
        set(FreqStandardDevThreshTextEdit,'String',...
            Standard_Deviation_Threshold_Default,'ButtondownFcn',@ButtonDown,...
            'Enable','on');
        
        


        % Determines the time scale of segment chosen in main window
        Freq_channel_time = t1:1/sample_rate:t2;

        % Plots the segment signal to the segment plot axes
        plot(FreqChannelAxes,Freq_channel_time,Freq_channel);
        grid(FreqChannelAxes);
        xlabel(FreqChannelAxes,'Time (s)');
        ylabel(FreqChannelAxes,'LFP amplitude(mV)');
        title(FreqChannelAxes,['Selected Epoch for ',reference_channel,' (',num2str(t1,...
            '%.0f'),'s - ',num2str(t2,'%.0f'),'s)']);
        xlim(FreqChannelAxes,[min(Freq_channel_time) max(Freq_channel_time)]);


         %highlights basal period in green if segment chosen in Main_Prog_Window
         %includes it
        if t1 < KA_injection_time == 1 && t2 > KA_injection_time == 1 ;
           
            
            ka_sig_length = KA_injection_time - t1;
            KA_plot_sig = Freq_channel(1:ka_sig_length*sample_rate+1);
            KA_plot_time = t1:1/sample_rate:KA_injection_time;
            
            %plots the chosen segment with the basal period displayed in a
            %different colour
            plot(FreqChannelAxes,Freq_channel_time,...
                Freq_channel,KA_plot_time,KA_plot_sig);
            grid(FreqChannelAxes);
            xlabel(FreqChannelAxes,'Time (s)');
            ylabel(FreqChannelAxes,'LFP amplitude (mV)');
            title(FreqChannelAxes,[' ',...
                reference_channel,' Epoch (',num2str(t1,...
                '%.0f'),'s - ',num2str(t2,'%.0f'),'s)']);
            legend(FreqChannelAxes,'Post-KA administration',...
                'Pre-KA administration', 'Location', 'North');
            xlim(FreqChannelAxes,...
                [min(Freq_channel_time) max(Freq_channel_time)]);
            
        elseif  t2 <= KA_injection_time == 1 ;
            
            
            %plots the basal segment in green
            plot(FreqChannelAxes,Freq_channel,...
                Freq_channel_time,'color',[0 0.5 0]);
            grid(FreqChannelAxes);
            xlabel(FreqChannelAxes,'Time (s)');
            ylabel(FreqChannelAxes,'LFP amplitude (mV)');
            title(FreqChannelAxes,['Segment of ',...
                reference_channel,' (',num2str(t1,...
                '%.0f'),'s - ',num2str(t2,'%.0f'),'s)']);
            legend(FreqChannelAxes,'Pre-KA administration',...
                'Location', 'North','Left');
            xlim(FreqChannelAxes,...
                [min(Freq_channel_time) max(Freq_channel_time)]);
            
            
            
        end
        
        
    end

    function FreqRunDetectButtonCallback(hObject,eventdata,handles)
        
         %Gets values of the text fields of parameters panel and assigns
        %values to variables
        Filter_Lower_Bound = str2double(get(FreqFilterLowerBoundTextEdit,...
            'String'));
        Filter_Upper_Bound = str2double(get(FreqFilterUpperBoundTextEdit,...
            'String'));
        Sliding_Window_Size = str2double(get(FreqSlidingWindowTextEdit,...
            'String'));
        Window_Overlap = str2double(get(FreqWindowOverlapTextEdit,...
            'String'));
        MinInter_Seiz = str2double(get(FreqMinInterSeizTextEdit,...
            'String'));
        MinSeiz_Duration = str2double(get(FreqMinSeizDurationTextEdit,...
            'String'));
        Standard_Dev_Thresh = str2double(get(FreqStandardDevThreshTextEdit,...
            'String'));
        Drop_band_select = get(FreqBandSelectDrop,'Value');
        
        Manual_band_select_Low = str2double(get(FreqManualLowBoundTextEdit,...
            'String'));
        Manual_band_select_Up = str2double(get(FreqManualUpBoundTextEdit,...
            'String'));
        
       
     
        %applies notch filter to channel and changes the mean_channel global variable
        %value which is the signal
        
        ChannelNotchFilter(Filter_Lower_Bound,Filter_Upper_Bound);
        
        %the value of I saved from this line is the nth datapoint of the 'time' vector 
        [min_dist, pre_KA_datapoints]=min(abs(Freq_channel_time-KA_injection_time));

        % this one is very slow with 'rand' data, and crashes with actual lfp data.
        % The only solution is to either vectorise, or use parpool function with
        % multiple processors.
        window_size = Sliding_Window_Size*sample_rate;
        lag = Window_Overlap*sample_rate;
        num_segs=length(Freq_channel)/lag;

        %finds last window
        a = find((1:lag:((lag*num_segs)-window_size)),1,'last');
        a = a+1;


        % create empty matrix to fill each column with fft of each window
        windows = zeros((1 + window_size),a);
        
        
        
        tic
        a=1;
        for n = 1:lag:((lag*num_segs)-window_size);
        windows(:,a) = (Freq_channel(n:(n+window_size)));
        time1(a) = Freq_channel_time(n);
        a = a+1;
        end
        
        
        % this accounts for the final extra count of 'a' which needs to be
        % reversed to keep vector lengths the same
        a = a-1;
        windows_fft = fft(windows);
        toc

        % %% 
        % %Brief code to plot a column (1st in this case) of fourier transformed data
        % 
        L=length(windows_fft(:,1)); 
        NFFT = 2^nextpow2(L); % Next power of 2 from length of y
        freq = sample_rate/2*linspace(0,1,NFFT/2+1);
        % window=windows_fft(:,1);
        % figure
        % plot(freq,2*abs(window(1:NFFT/2+1)));% Fig 2: Discrete absolute fourier transform of column of fft data.
        % title('Example plot of FFT of 1st extracted column of LFP data');
        % xlabel('Frequency (Hz)');
        % ylabel('Power (mV)');
        % xlim ([0 200]);

        % This section of code integrates only under the user selected frequency band
        if Drop_band_select == 1; %delta
            b1 = 0.4;
            b2 = 4;
        elseif Drop_band_select == 2%theta
            b1 = 4;
            b2 = 8;
            elseif Drop_band_select == 3%alpha
                b1 = 8;
                b2 = 12;
                elseif Drop_band_select == 4 %beta
                    b1 = 12;
                    b2 = 30;
                    elseif Drop_band_select == 5 %gamma
                        b1 = 30;
                        b2 = 100;
                        elseif Drop_band_select == 6 %broad
                            b1 = 0.4;
                            b2 = 100;           

        
        %if manual bounds entered, set lower and upper bounds based on the
        %users selection
        else
            b1 = Manual_band_select_Low;
            b2 = Manual_band_select_Up;
        end
    

        %finds the index of the frequency closest to the frequency chosen
        %by the user or the chosen band.
        [~, index] = min(abs(freq-b1));
        [~, index1] = min(abs(freq-b2));

         
        for n = 1:a
        window = windows_fft(:,n);
        fft_column_data = 2*abs(window(1:NFFT/2+1));
        power = (fft_column_data(index:index1)).^2;%gives power
        fft_win_sum(n) = sum(power);
        end
        
        for l=1:length(KA_injection_time)
            
             [min_dist, pre_KA_datapoints_1]=min(abs(time1-KA_injection_time));%the value of I saved from this line is the nth datapoint of the 'time' vector 
             %for this particular channel when it is closest to the inputted 'KA_injection_time'
        end
        
       
        mean_val = sum(fft_win_sum(1:pre_KA_datapoints_1)/pre_KA_datapoints_1);
        deviation = std(fft_win_sum(1:pre_KA_datapoints_1));

        threshold = mean_val + (Standard_Dev_Thresh*deviation); 
        threshold_lines = zeros(length(time1),1);


        %This is if threshold pop out is wanted!!!
        % for a= 1:length(time1)
        %     threshold_lines(a)=threshold;
        % end
        % 
        % figure; %Fig 3: FFT integral plot
        % plot(time1,fft_win_sum);
        % title('FFT integral values for each data window vs. time')
        % xlabel('Time (s)');
        % ylabel('Integral values of data window frequency spectra (mV)');
        % 
        % hold on
        % 
        % 
        % plot(time1,threshold_lines,'r');

        
        data.seiz=zeros(length(Freq_channel),1);

        trigger = find((fft_win_sum(1:end))>threshold); %finds indices of when window_average is above the set threshold

        x =(trigger*lag)-(lag-1);%corresponding indices of data.seiz(Freq_channel) where window starts
        for k=1:length(x)
            if (x(k)+window_size)<length(data.seiz)% sets the condition so that the length of data.seiz isnt altered by this condition
         data.seiz(x(k):(x(k)+ window_size))=1;% sets the data points in data.seiz to 1, corresponding to those contained between each window above threshold from 'trigger'.
            else
                end_diff=(x(k)+window_size)-length(data.seiz);%obtains number of data points from 'triggered' above threshold window start point in data.seiz to end of data.seiz
                data.seiz(x(k):(x(k)+ end_diff))=1;%sets these values to 1 also without altering the length of data.seiz.
            end
        end


      on_off=find(abs(diff(data.seiz))>0);
      data_dist=diff(on_off);
       
      
      
         
      data.seiz(1:3)=0;%  this ensures that the ordering of the algorithms selection of seizure start and end points holds true for data starting with a seizure.
      data.seiz(end)=0; %this ensures that even if seizue continues beyond dataset, the algorithm can detect and measure its duration within data selected.
   
      
%    if isempty(find(data_dist(2:2:end)<MinInter_Seiz*sample_rate,1)) == 0 % i.e. do this if it's not empty (if inequality is met)
           
      % extracts indices of array data_dist of interictal events above user defined min duration.
        interictal=find(data_dist(2:2:end)<MinInter_Seiz*sample_rate);
        
     
                %This next lines extract the start point of each interictal
                %interval that exceeds the time threshold individually using the
                %command 'interictal(int_num)' acting upon 'on_off'' which
                %contains the start/end locations of each event. The
                %algebra '*2' ensures selection of the start of the 'off'
                %location of seizures in the 'on_off' array, which is the
                %start of the interictal period. All of these array indices
                %are subseqently selected from the data.seiz array and set
                %back to 0 from 1, giving an array of 1's and zero's for
                %selected seizure events. 
        for int_num = 1:length(interictal) 
  data.seiz(on_off(interictal(int_num)*2):(on_off(interictal(int_num)*2+1))) = 1; 
        end
%    else
%         disp('Equality in line 140 not met; interictal period above selected threshold duration not found within data set');
%        end
    
  % NOW REMOVE ALL FLASE SEIZURES FROM DATA.SEIZ ARRAY
  
    %This extracts indices of where in the data the threshold is crossed,indicating either the start or end of a seizure.
      on_off_1=find(abs(diff(data.seiz))>0);
      data_dist_1=diff(on_off_1);
      
      
      %DO WE NEED THE NEXT LINE OF CODE?? THE SUBSEQUENT OPERATIONS ARENT
      %DEPENDANT ON THIS RELATION HOLDING TRUE
%       if isempty(find(data_dist_1(1:2:end)>sample_rate*MinSeiz_Duration,1)) == 0 % i.e. do this if it's not empty (if inequality is met)    
      
        seizure=find(data_dist_1(1:2:end)<sample_rate*MinSeiz_Duration);% extracts indices of array data_dist of seizure events BELOW user defined min duration.
        
        %The next two lines extract the start point of each seizure event
        %that is less than the user defined minimum MinSeiz_Duration
        %individually using the command 'seizure(sez_num)' acting upon 'on_off_1' which
        %contains the start/end locations of each event. The
        %algebra '*2-1' ensures selection of the start of the 'on'
        %location of seizures in the 'on_off_1' array, which is the
        %start of each seizure period. All of these array indices
        %are subseqently selected from the data.seiz array and set
        %back to 0 from 1, giving an array of 1's and zero's for
        %remaining thresholded seizure events.
        for sez_num = 1:length(seizure) 
  data.seiz(on_off_1(seizure(sez_num)*2-1):(on_off_1(seizure(sez_num)*2))) = 0; 
        end
%       else
%           disp('Equality in line 172 not met, no seizures above selected duration found within data set');
%   
%       end
    
    Data_Seizure = data.seiz;
    
    SeizureEventDisplay(Freq_channel_time,Data_Seizure);
    
    %% Rewrite another method for seizure detection here, based on a certain number of consecutive windows above the threshold value.


% %% Output of final seizure data to user
% 
%  data.seiz_real(:,i) = data.seiz(1:length(time));% this one works better with multiple channel analysis to store each final output.
%  
%     data.seiz_n(i) = length(find((diff(data.seiz))>0));% stores the final number of detected seizures present in 'data.seiz' 
%     %(note that 'abs' isnt used as a prequel, therefore this will only select start points of seizures.
%     
%     data.seiz_total_dur(i) = sum(data.seiz)/sample_rate; % this works as elements of data.seiz are either 1 or 0, so dividing total sum of array by sample_rate gives total seizure time
%     
% 
%  if data.seiz_n > 1
%         seiz_start = find(abs(diff(data.seiz))>0); %outputs indices of each seizure start and end point
%         for seiz = 1:data.seiz_n(i)% cycles through each seizure
%             data.seiz_dur(seiz) = (seiz_start(seiz*2)-seiz_start(seiz*2-1))/sample_rate;
%         end
%     else
%         seiz_start = find(abs(diff(data.seiz))>0);
%  end
%     
%     mean_dur=mean(data.seiz_dur);
%     standard_dev=std(data.seiz_dur);
%     standard_error=standard_dev/(sqrt(data.seiz_n));
%     
%     
%     disp(' ')
%     disp(['Variable ' nth_trace(i,:) ', Overall Seizure Duration = ' num2str(data.seiz_total_dur(i),'%.4f') 's'])
%      disp(' ')
%     disp(['Variable ' nth_trace(i,:) ', Number of Seizure Events = ' num2str(data.seiz_n(i))])
%      disp(' ')
%     disp(['Variable ' nth_trace(i,:) ', Mean seizure duration = ' num2str(mean_dur,'%.4f') 's'])
%      disp(' ')
%     disp(['Variable ' nth_trace(i,:) ', Standard Error of Mean = ' num2str(standard_error,'%.4f') 's'])
%      disp(' ')
%     if isempty(seiz_start) == 1
%         disp(' '); disp('*** NO SEIZURES DETECTED WITH CURRENT THRESHOLD SETTING ***')
%     else
%         
%         k=1;
%         for l=1:2:length(seiz_start)
%             
%         disp(['Variable ' nth_trace(i,:) ', start time of seizure number ' num2str(k) ' = ' num2str((time(seiz_start(l))),'%.4f') 's'])
%         disp(['Variable ' nth_trace(i,:) ', duration of seizure number ' num2str(k) ' = ' num2str(data.seiz_dur(k),'%.4f') 's'])
% %         format short g;
% %         if data.seiz_n == 1
% %             data.seiz_total_dur
% %             
% %         else
% %             data.seiz_dur
% %       
% %         end
% start_time(k)=(time(seiz_start(l)));
% start_time=start_time';
% duration(k)=data.seiz_dur(k);
% duration=duration';
%         k=k+1;
%         end
%         disp(['NB Calculations utilise only seizures > ' num2str(MinSeiz_Duration) 's in duration; seizures occurring less than ' num2str(MinInter_Seiz) 's apart are considered as one seizure'])
%         format short g
%         
%     end
%     
%     if export == 'y'
%         
% start={'Seizure Start Time'};
% xlswrite('seizure_events.xls',start,'seizure_events','B1')
% xlswrite('seizure_events.xls',start_time,'seizure_events','B2')
% 
% dur={'Seizure Duration'};
% xlswrite('seizure_events.xls',dur,'seizure_events','C1')
% xlswrite('seizure_events.xls',duration,'seizure_events','C2')
%     else
%     end


       %plots spectrogram 
        SpectroPlot();
    end


   function SpectroPlot

        F=linspace(1,sample_rate/4,sample_rate/4);      


        [S,F,T,P] = spectrogram(Freq_channel,sample_rate,0,F,sample_rate);
        surf(FreqSpectrogramAxes,T,F,10*log10(P),'edgecolor','none'); axis tight;
        view(FreqSpectrogramAxes,0,90);
        xlabel(FreqSpectrogramAxes,'Time (s)'); ylabel('Frequency (Hz)');
        shading interp
        chandle= colorbar('YTickLabel', [0:10]);
        ylabel(chandle, 'Frequency Intensity (mV s ^-1)');
        title(FreqSpectrogramAxes,['Spectrogram of Epoch ',...
                        reference_channel,' (',num2str(t1,...
                        '%.0f'),'s - ',num2str(t2,'%.0f'),'s)']);
        
        
    end

    function SeizureEventDisplay(Freq_channel_time,Data_Seizure)
        
             
    %seizure detect Display
    plot(FreqDetectAxes,Freq_channel_time,Data_Seizure,'r','LineWidth',1.6);
    title(FreqDetectAxes,'Seizure Events & Duration');
    xlabel(FreqDetectAxes,'Time (s)');
    grid(FreqDetectAxes);
    ylabel(FreqDetectAxes,'1 = seizure event');
    ylim (FreqDetectAxes,[ -0.2 1.2]);
    xlim(FreqDetectAxes,[min(Freq_channel_time) max(Freq_channel_time)]);
    
    end

 function ExcelExport(hObject,eventdata)

        start={'Seizure Start Time'};
        xlswrite('seizure_events.xls',start,'seizure_events','B1')
        xlswrite('seizure_events.xls',start_time,'seizure_events','B2')

        dur={'Seizure Duration'};
        xlswrite('seizure_events.xls',dur,'seizure_events','C1')
        xlswrite('seizure_events.xls',duration,'seizure_events','C2')

            
            
 end
    
    end