function varargout = interPSvsTime(varargin)
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%     G U I   H E L P   I N F O R M A T I O N                                   % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%  ourGUI31 M-file_menu for ourGUI31.fig                                        % 
%                                                                               % 
%       The GUI is designed to produce graphs of the power between two          % 
%       frequencies as a function of time. It should be launced from the        % 
%       'Front_Panel' window of the 'LFP analysis suite'.                       % 
%                                                                               % 
%       The algorithm selects a section of the signal stored as handles.userdata% 
%       and produces a periodogram using the signal processing toolbox command  % 
%       The periodogram is stored as the columns of an array aclled STFT which  % 
%       is converted to a PSD data object so that Matlab's avgpower command     % 
%       can be used on the columns (which represent the signal at different     %
%       times)                                                                  % 
%       
%       The length window that the user selects is used to select,through
%       multiplication, a parection of th e
%
%       The GoBtn function is called when the user clicks 'go' and this calls   % 
%       ,Amongst other things, the getPSD function which gets the periodogram   % 
%       PSD to generate the plots                                               % 


%       Note taht the PSD is an ESTIMATE of the actual frequency content of
%       the signal and so can be very incorrect in magnitude, especially
%       for shorter window lengths. But longer window lengths make the time
%       accuracy of the estimation inaccurate
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %     
%                                                                               % 
%  This Software is distributed under the GNU General Public Liscence           % 
%  a copy of this liscence sould be accompany this software.                    % 
%                                                                               % 
% Copyright (C) 2006    Richard Peter                                           %
%                                                                               %  
% This program is free software; you can redistribute it and/or modify          %
% it under the terms of the GNU General Public License as published by          %
% the Free Software Foundation; either version 2 of the License, or             %
% (at your option) any later version.
%                                                                               % 
% This program is distributed in the hope that it will be useful,               %
% but WITHOUT ANY WARRANTY; without even the implied warranty of                %
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                 %
% GNU General Public License for more details.                                  %
%                                                                               %  
% You should have received a copy of the GNU General Public License             %
% along with this program; if not, write to the Free Software                   %
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301 USA %
%                                                                               % 
%                                                                               % 
%  For further information contact Richard Peter at                             % 
%  infinitybeckon@hotmail.com                                                   % 
%                                                                               %   
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @interPSvsTime_OpeningFcn, ...
                   'gui_OutputFcn',  @interPSvsTime_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


% --- Executes just before interPSvsTime is made visible.
function interPSvsTime_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to interPSvsTime (see VARARGIN)

% Choose default command line output for interPSvsTime
handles.output = hObject;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialisation Code Goes Here

%%%%%%%%%%%%%%%%
%Userdata

%Read in the variables that have been passed to it from the front panel
%window 
handles.userdata = varargin{1};
handles.user_start = varargin{2};
handles.user_end = varargin{3};
handles.tstep = varargin{4};
%%%%%%%%%%%%%%%%

handles.Numfigs = 0;
%Define required variables
handles.NotchVals = [];
handles.CustomBands = [];

% Reset all Toggle Buttons/Checkboxes
refresh_window(handles);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update handles structure
guidata(hObject, handles);

% UIWAIT makes interPSvsTime wait for user response (see UIRESUME)
% uiwait(handles.figure1);


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

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


% --- Executes when the user cheks the Select all Bands Callback
% --- This simply select all the pre-defined frequency Bands
function SelectAllBands_Callback(hObject, eventdata, handles)
%Check to see if the user has selected/deselected
if(get(hObject,'Value') == get(hObject,'Max'))
    % Then checkbox is checked- select all defalut frequency bands
    set(handles.AlphaBtn,'Value',1);
    set(handles.BetaBtn,'Value',1);
    set(handles.DeltaBtn,'Value',1);
    set(handles.ThetaBtn,'Value',1);
    set(handles.GamLBtn,'Value',1);
    set(handles.GamHBtn,'Value',1);
    set(handles.GamVHBtn,'Value',1);
else 
    % The user wants top deselect all Frequency Bands
    set(handles.AlphaBtn,'Value',0);
    set(handles.BetaBtn,'Value',0);
    set(handles.DeltaBtn,'Value',0);
    set(handles.ThetaBtn,'Value',0);
    set(handles.GamLBtn,'Value',0);
    set(handles.GamHBtn,'Value',0);
    set(handles.GamVHBtn,'Value',0);
end %End If

% save the handles structure 
guidata(hObject, handles);
% End of Function            



% --- Executes on button press in ChooseBandsTgl.
function ChooseBandsTgl_Callback(hObject, eventdata, handles)
% hObject    handle to ChooseBandsTgl (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 ChooseBandsTgl
% {1 = Pressed, 0 = Not Pressed}
if(get(hObject,'Value') == 1)
    % We first check whether the toggle button has been turned on, 
    % Then setnotchtgl the radio Button to selected (1)
    % Then get the user's bands
    [handles.CustomFilters, Y] = CustomFilters(handles);
    set(handles.ChooseBandsTgl, 'Value', Y);
    handles;
else %If its Been Un-Pressed 

end %End if

% Update handles structure
guidata(hObject, handles);

%End of Function





% --- Executes on button press in DefaultWdowBtn.
function DefaultWdowBtn_Callback(hObject, eventdata, handles)
% Hint: get(hObject,'Value') returns toggle state of DefaultFilterBtn
if get(hObject, 'Value') == 1
    %If its on
set(handles.OwnWdowBtn, 'Value', 0);          %deselect own filter radio button
set(handles.DefaultWdowList, 'Enable', 'on'); %Disable Dropdown menu
set(handles.OwnWdowList, 'Enable', 'off'); %Enable Default filter list
else
    %if its off
set(handles.OwnWdowBtn, 'Value', 1);          %select own filter button
set(handles.DefaultWdowList, 'Enable', 'off');%Enable Drop down menu
set(handles.OwnWdowList, 'Enable', 'on');  %Disable Default filter list
end


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

usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end





% --- Executes on button press in SetNotchTgl.
% --- Here the user wants to setnotchtgl some Notch filters, perhaps to remove Noise
function SetNotchTgl_Callback(hObject, eventdata, handles)
% Hint: get(hObject,'Value') returns toggle state of SetNotchTgl

if(get(hObject,'Value') == get(hObject,'Max'))
    % We first check whether the toggle button has been turned on, 
    %Next we create an empty vector of Notch limits
    %Now we launch the NotchSet GUI
    [handles.NotchVals, Y] = NotchSet(handles);
    set(handles.SetNotchTgl, 'Value', Y);
end %End if

%Save the handles structure
guidata(hObject, handles);
% End of Function

% --- Executes on button press in ViewWdowBtn. 
%The user clicks this when they want to view their window and its FFT
function ViewWdowBtn_Callback(hObject, eventdata, handles)
% hObject    handle to ViewWdowBtn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
wvtool(GenerateWdow(handles));

% --- Executes on button press in GoBtn.
% When the user wants to plot they click this.
function GoBtn_Callback(hObject, eventdata, handles)
% hObject    handle to GoBtn (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
if get(hObject, 'Value')==1 %The user has pressed the button...
if get(handles.SetNotchTgl, 'Value')


        Coeffs = fir1(1000, [handles.NotchVals]/handles.tstep*2, 'stop');
        handles.userdata = filtfilt(Coeffs, 1,  handles.userdata);
end

handles.Lims = [];
handles.Lims = GenerateLims(handles);
%Lims is now an n X 2 array of the low/high frequencies that the user has
%selected.


% %Next generate filters to allow only these frequencies
handles.PSD = getPSD(handles);
% handles.Filters = GenerateFilters(handles);

if length(handles.PSD) >0
    %Now plot it
    handles.plotme = [];
    if(size(handles.Lims, 1) ~= 0 )
        % the user wants to see some separate signals, so lets filter the
        % signal and then store it for plotting
        
            % If usng the spectrogram method :The FFt order is 4096
%               handles.Lims = ceil(handles.Lims*size(handles.PSD, 1)/(handles.tstep/2));
            for ii = 1:size(handles.Lims, 1)
%                 handles.plotme(ii, :) = sum(handles.PSD(...
%                     ceil(handles.Lims(ii, 1)):ceil(handles.Lims(ii, 2)), :), 1);
                handles.plotme(ii, :) = avgpower(handles.PSD, [handles.Lims(ii, 1) handles.Lims(ii, 2)]);
                
            end
      
        
        %now the plotting bit...
        %We ned to make a new plot and keep track of how many there are...
        %then make an x-axis
        %Then plot each setnotchtgl of data on its own subplot
        %then setnotchtgl a title and x-limits for each axis
        % finally, if its the last axis then add an x- axis label
        fighandle = figure (handles.Numfigs+1); % Generate a figure and store handle
        handles.Numfigs=handles.Numfigs+1; %Increment the number of figures
        h = newplot(fighandle); %Generate a new plot in the figure
        Axislabels = getaxislabels(handles); %Get a cell array of axis labels


        %next make an x-axis to plot against...
        %Define the increment per plot point
        STFTstep = str2num(get(handles.TimeRes, 'String'))*handles.tstep;
        %Define the number of points in a window
        hannN = str2num(get(handles.WdowLen, 'String'))*handles.tstep;
        %Define the number of user's points
        tN = length(handles.userdata);
        
        ModStart = handles.user_start+str2num(get(handles.WdowLen, 'String'))/2;
        ModEnd = (floor((tN-hannN)/STFTstep)-1)*(STFTstep)/handles.tstep+ModStart;
        
        
        f = linspace(ModStart, ModEnd, floor((tN-hannN)/STFTstep)+1);
        %If the user wanted to see the original signal then make a note...
        m = 0; %This variable will allow us to plot the original signal first
        if get(handles.PlotOrigSig, 'Value')
            m=1; %if they want the original signal add 1 onto everything else...
            %Now make an x- axis
            dataaxis = linspace(handles.user_start, handles.user_end,length(handles.userdata));
            %Next get the handle for the subplot and plot it
            h = subplot(size(handles.Lims, 1)+m,1,1); % get the handle of the subplot
            plot(h, dataaxis, handles.userdata); %Plot the graph
            xlim(h, [ModStart, ModEnd]);
            title(h, 'Original Data');
            ylabel(h,  'Signal (mV)');
        end
        for t = 1+m:size(handles.Lims, 1)+m
            h = subplot(size(handles.Lims, 1)+m,1,t); % get the handle of the subplot
            plot(h, f, handles.plotme(t-m, :)); %Plot the graph
            xlim(h, [ModStart, ModEnd]);
            title(h, Axislabels{t-m});
            ylabel(h,  'Signal (mV^2)');

            if(t -m== size(handles.Lims, 1))
                xlabel(h, 'Time (s)'); %If its the last graph, label the x-axis
            end % if
        end % for

    else
        errordlg('You must select a frequency band',...
            'Incorrect Selection','modal')
        set(hObject, 'Value', 0);
    end

end %If length PSD > 0
set(hObject, 'Value', 0);
end %If button is down

% Update handles structure
guidata(hObject, handles);
%End of Function


%This Function Simply sets the Button values back to their default values.
%It is useful to know when a Window has been refreshed whilst designing 
function refresh_window(handles)
    %First setnotchtgl the frequecny bands Box...
    %FTo start check the SelectallBands box
    set(handles.SelectAllBands, 'Value', 0);
    %Now when we call its callback it will clear the selection of all the
    %Frequency Bands Radio buttons
    SelectAllBands_Callback(handles.SelectAllBands, handles, handles)
    %next the rest of the Frequency bands Panel
    set(handles.JnAdjBands, 'Value', 0);
    set(handles.ChooseBandsTgl,'Value', 0);
    set(handles.SetNotchTgl,'Value', 0);

  
    %Now the Window Panel
    set(handles.DefaultWdowList,'Enable', 'on');

%End of function
    
    
    
% This Function Gets the Limits filters requested...
% Let us ensure that the output is a cell array of filter coefficients where 
% the top row is B coefficients and the second row is A coefficients called 
% 'Coefficients'.
%
% In addition an array (Lims) comprising a two column vector  of 
% frequency ranges. These will be the
% start and stop frequencies of each of the filters designed .
%
% The number of rows in the 2-column vector (Lims) will therefore equal the
% number of elements in the cell array of coefficients (Coefficients).

function [Lims] = GenerateLims(handles)
%First create boolean variables to describe the requirements:
Join = get(handles.JnAdjBands, 'Value');        %Join Adjacent bands?
CustomT = get(handles.ChooseBandsTgl, 'Value'); %User's Custom Bands?
NotchT = get(handles.SetNotchTgl, 'Value');     %Notch Filters eg. 50Hz?

%Next get the frequency bands that are selected...
%Bands is an array of start/stop frequencies
%BandBool is a bollean array of which bands are selected.
[Bands, BandBool] = GetFrequencyBands(handles);
Bands = Bands';
%Next we see how many filters are required (Numfilts)
LBool = bwlabel(BandBool);
Numfilts = max(LBool);


% If the user wants to join frequency bands
if Join ==1 && Numfilts>0
    %The user wants to join the frequency Bands together...
    ii = [];
    NewBands = [];
    for ii= 1:Numfilts
        F = find(LBool == ii);    %Finds all adjacent entries
        Min = min(Bands(F, 1)); %Find the minimum of the minimums 
        Max = max(Bands(F, 2)); %Find the maximum on the maximums
        NewBands(ii, (1: 2)) = [Min Max];
    end %For
    %NewBands are the uppr and lowe limits of each requested  band in a 
    % nx2 array. Where n is the number of required filters, Numfilts. 
    count = ii;
else
    ii = [];
    count =1;
    NewBands = [];
    F = find(Bands(:, 1)== 999);
    if length(F) ==0
        F=999;
    end
    %F is now the places of the unwanted frequency Bands within Bands
    for ii = 1:size(Bands, 1)
        if ii ~= F
            %If its not in F then add it to NewBands
           NewBands(count, (1: 2)) = Bands(ii, 1:2); 
           count = count+1;
        end % If 
    end %for ii
    %NewBands are the uppr and lowe limits of each band in a nx2 array.
    % Where n is the number of required filters.
    count = count-1;
end %if Join =1

%If the User has defined custom Bands
if CustomT == 1
    NumCust = length(handles.CustomFilters)/2; %NumCust=the number defined
    for ii = 1:NumCust   
        %Append each filter's limits to NewBands   
        NewBands(count+1, (1:2)) = handles.CustomFilters(ii:ii+1);
        %increment the count
        count = count+1;
    end %for
end %if

%setnotchtgl Lims to NewBands as these are the limits we wish to create filters for
Lims = NewBands; 

%End Of Function
    
    
% --- This function generates a 2xn array of frequency limits (from the
% --- frequency bands selected)
function [Bands, Bool] = GetFrequencyBands(handles)
Bands = []; %initialise the Bnads Variable 
count = 1; %Start Counting the required bands
% %Now check to see if each band is required
Bool = Bandmask(handles);
Freqs = [1 4; 4 8; 8 14; 14 30; 30 50; 50 100; 100 150]';
Bands(1, :) = Freqs(1, :).*Bool;
Bands(1, :) = (Bands(1, :)==0).*999 +Bands(1, :);
Bands(2, :) = Freqs(2, :).*Bool;
handles;

% End Of Function
    

    

%This function determines the filter that is selected from the default box
function var1 = GetSelected(handlesItem)


index_selected = get(handlesItem,'Value');
if length(index_selected) ~= 1
    errordlg('You must select only one variable',...
            'Incorrect Selection','modal')
        set(handles.GoBtn, 'Value', 0);
else
    var1 = index_selected(1);
end 

%This function generates a cell array of axis labels for use when plotting
%It is called by the GoBtn_Callback function
function AxLabels = getaxislabels(handles)
[Bands Bool] = GetFrequencyBands(handles);
Bands = Bands';
LBool = bwlabel(Bool);
ii = [];
clear('AxLabels');
if get(handles.JnAdjBands, 'Value')
    for ii = 1:max(LBool)
        %To save on memory, the following two variables have been put
        %directly into the AxLabels assignemnt line
        %Min = min(Bands( find(LBool==ii), 1));
        %Max = max(Bands( find(LBool==ii), 2));
        AxLabels{ii} = strcat(num2str(min(Bands( find(LBool==ii), 1))), '-',...
            num2str(max(Bands( find(LBool==ii), 2))), ' Hz');
    end
else
    %Find the elements that are needed
    F = find(Bool==1);
    %Now go through each elelment and add its label
    for ii = 1:sum(Bool)
        AxLabels{ii} = strcat(num2str(Bands(F(ii),1)), '-',...
            num2str(Bands(F(ii),2)), ' Hz');
    end
end

%Now add the labels of custom bands
if get(handles.ChooseBandsTgl, 'Value')
    F=handles.CustomFilters;
    iii = [];
    for iii = 1:2:length(F)
        Min = F(iii);
        Max = F(iii+1);
        ii = ii+1;
        AxLabels{ii} = strcat(num2str(Min), '-', num2str(Max), ' Hz');
    end
end
%Generates a mask of the Frequency bands selected 1 = selected, 0 = not
%selected. It is called by GetFrequencyBands and getaxislabels.
function Mask = Bandmask(handles)

Mask = [get(handles.DeltaBtn, 'Value'), get(handles.ThetaBtn, 'Value'),...
    get(handles.AlphaBtn, 'Value'), get(handles.BetaBtn, 'Value'),...
    get(handles.GamLBtn, 'Value'), get(handles.GamHBtn, 'Value'),...
    get(handles.GamVHBtn, 'Value')];


% --- Executes on button press in DontFilterCQbox.
function NULL(hObject, eventdata, handles)
% hObject    handle to DontFilterCQbox (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 DontFilterCQbox

handles;



% --- Executes during object creation, after setting all properties.
function NULL_CreateFcn(hObject, eventdata, handles)

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


%When the User edits their desired time resolution - update the overlap
function TimeRes_Callback(hObject, eventdata, handles)
% hObject    handle to TimeRes (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)




% --- Executes during object creation, after setting all properties.
function TimeRes_CreateFcn(hObject, eventdata, handles)
% hObject    handle to TimeRes (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


%This Function Generates an array for freq x time spectrogram surface.
function PSD = getPSD(handles)

handles.Wdow = GenerateWdow(handles);
STFTstep = str2num(get(handles.TimeRes, 'String'))*handles.tstep;  
hannN = str2num(get(handles.WdowLen, 'String'))*handles.tstep;%number of points in window
win = handles.Wdow;


% %%% Periodogram
handles.WelchWinLength = hannN;
Win = GenerateWdow(handles);  
 
% Next get the number of samples in the user's selection
tN = length(handles.userdata); % or (handles.user_end - handles.user_start)*handles.tstep;


NFFT = 4096;
STFT = [];

for pp = 1:floor((tN-hannN)/STFTstep)+1

    %Periodogram
        windowedSIG = handles.userdata(floor((pp-1)*(STFTstep)+1):...
        ceil((pp-1)*(STFTstep)+hannN)).*Win';
    [Pxx(:, pp), f] = periodogram(windowedSIG,[], NFFT,handles.tstep, 'onesided');

end
% if length(Pxx) ==0
if length(STFT) ==0
    errordlg('Error. Please Check that the window is shorter than the data',...
        'Incorrect Selection', 'modal');
    set(handles.GoBtn, 'Value', 0);
    PSD = [];
else
    PSD = dspdata.psd(Pxx, f, 'Fs', handles.tstep);

end %if



%This function generates the user's window.
function Wdow = GenerateWdow(handles)
%first get the selected window 
windownum = get(handles.DefaultWdowList, 'Value');

WdowLen = round(str2num(get(handles.WdowLen, 'String'))*handles.tstep);
%find out which window they want and make it
if windownum == 1
    Wdow = hann(WdowLen);
elseif windownum == 2
    Wdow = hamming(WdowLen);    
elseif windownum == 3
    Wdow = rectwin(WdowLen);    
elseif windownum == 4
    Wdow = bartlett(WdowLen);
elseif windownum == 5
    Wdow = barthannwin(WdowLen);    
elseif windownum == 6
    Wdow = blackman(WdowLen);
elseif windownum == 7
    Wdow = gausswin(WdowLen);
elseif windownum == 8
    Wdow = kaiser(WdowLen, 2.5);
elseif windownum == 9
    Wdow = tukeywin(WdowLen, 0.25);
elseif windownum == 10
    Wdow = tukeywin(WdowLen, 0.5);
elseif windownum == 11
    Wdow = tukeywin(WdowLen, 0.75);
end
%End of function