% POP_REJSPEC - rejection of artifact in a dataset using 
%                 thresholding of frequencies in the data.
% Usage:
%   >>  pop_rejspec(INEEG, typerej); % pop-up interactive window mode
%   >> [OUTEEG, Indices] = pop_rejspec( INEEG, typerej, 'key', val, ...);
%
% Pop-up window options:
% "Electrode|Component" - [edit box] electrode or component indices to 
%                 take into consideration for rejection. Sets the 'elecrange'
%                 parameter in the command line call (see below).
% "Spectrum computation method" -[popupmenu] with values {'fft','multitaper'}
%                 Select method to compute spectrum.
% "Minimum power rejection threshold(s) (dB)" - [edit box] minimum rejection
%                 threshold(s) in (in dB). Sets the command line parameter 
%                 'threshold'. If more than one, apply to each electrode|
%                 component individually. If fewer than number of 
%                 electrodes|components, apply the last value to all remaining
%                 electrodes|components.
% "Maximum power  rejection threshold(s) (dB)" - Maximum rejection threshold(s) 
%                 Sets the command  line parameter 'threshold'.
% "Low frequency  limit(s)" - [edit box] Low-frequency limit(s) in Hz of 
%                 range to perform  rejection. Sets the command line parameter
%                 'freqlimits'.
% "High frequency limit(s)" - [edit box] High-frequency limit(s) in Hz of 
%                 range to perform  rejection. Sets the command line parameter
%                 'freqlimits'.
% "Display previous rejection marks: " - [Checkbox]. Sets the command line
%                 input option 'eegplotplotallrej'.
% "Reject marked trials: " - [Checkbox]  Sets the command line
%                 input option 'eegplotreject'.
%
% Command line inputs:
%   INEEG      - input dataset
%   typerej    - [1|0] data to reject on (0 = component activations; 1 = 
%                electrode data). {Default is 1}. 
%
% Optional arguments.
%   'elecrange'     - [e1 e2 ...] array of indices of electrode|component 
%                     number(s) to take into consideration during rejection.
%   'threshold'     - [lower upper] threshold limit(s) in dB.
%   'freqlimits'    - [lower upper] frequency limit(s) in Hz.
%   'method'        - ['fft'|'multitaper'] method to compute spectrum.
%   'specdata'      - [array] precomputed spectral data.
%   'eegplotcom'    - [string] EEGPLOT command to execute when pressing the
%                     reject button (see 'command' input of EEGPLOT).
%   'eegplotreject' - [0|1] 0 = Do not reject marked trials (but store the 
%                     marks. 1 = Reject marked trials. {Default: 1}.
%   'eegplotplotallrej' - [0|1] 0 = Do not superpose rejection marks on previous
%                     marks stored in the dataset. 1 = Show both previous and
%                     current marks using different colors. {Default: 0}.
%
% Outputs:
%   OUTEEG     - output dataset with updated spectrograms
%   Indices    - index of rejected trials
%   Note: When EEGPLOT is called, modifications are applied to the current 
%   dataset at the end of the call to EEGPLOT (e.g., when the user presses 
%   the 'Reject' button).
%
% Author: Arnaud Delorme, CNL / Salk Institute, 2001-
%
% See also: EEGTHRESH, EEGLAB, EEGPLOT, POP_REJEPOCH

% Copyright (C) 2001 Arnaud Delorme, Salk Institute, arno@salk.edu
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.

% 01-25-02 reformated help & license -ad 
% 03-07-02 added srate argument to eegplot call -ad
% 03-08-02 reworked spectrum to save space & add eeglab options -ad

function [EEG, Irej, com] = pop_rejspec( EEG, icacomp, varargin);
    %elecrange, negthresh, posthresh, ...
   	%startfreq, endfreq, superpose, reject);

Irej = [];
com  = '';
if nargin < 1
   help pop_rejspec;
   return;
end
if nargin < 2
   icacomp = 1;
end  

if icacomp == 0
	if isempty( EEG.icasphere )
	    ButtonName=questdlg( 'Do you want to run ICA now ?', ...
                         'Confirmation', 'NO', 'YES', 'YES');
    	switch ButtonName
        	case 'NO', disp('Operation cancelled'); return;   
        	case 'YES', [ EEG, com ] = pop_runica(EEG);
    	end % switch
	end
end	

if nargin < 3
    
    % which set to save
    % -----------------
    promptstr = { fastif(icacomp==0, 'Component indices (Ex: 2 4 6:8 10)', ...
                                     'Electrode indices (Ex: 2 4 6:8 10)'), ...
                  'Spectrum computation method',...
                  'Minimum power rejection threshold(s) (dB)',... %'Low power rejection threshold (s) (dB)', ...
                  'Maximum power rejection threshold(s) (dB)',... %'High power rejection threshold (s) (dB)', ...
                  'Low frequency limit(s) (Hz)', ...
                  'High frequency limit(s) (Hz)', ...
                  'Display previous rejection marks', ...
                  'Reject marked trial(s)' };
    
    inistr = { ['1:' int2str(EEG.nbchan)], ...
               '',...
               '-30', ...
               '30', ...
               '15', ...
               '30', ...
               '0', ...
               '0' };
    
    methodlist = {'FFT', 'Multitaper'};
    g  = [1 0.1 0.75];
    g2 = [1 0.26 0.9];
    g3 = [1 0.22 0.85];
    geometry = {g g2 g g g g 1 g3 g3};
    uilist = {...
              { 'Style', 'text', 'string', promptstr{1}} {} { 'Style','edit'      ,'string' ,inistr{1}  'tag' 'cpnum'}...
              { 'Style', 'text', 'string', promptstr{2}} {} { 'Style','popupmenu' ,'string' ,methodlist 'tag' 'specmethod' }...
              { 'Style', 'text', 'string', promptstr{3}} {} { 'Style','edit'      ,'string' ,inistr{3}  'tag' 'lowlim'}...
              { 'Style', 'text', 'string', promptstr{4}} {} { 'Style','edit'      ,'string' ,inistr{4}  'tag' 'higlim'}...
              { 'Style', 'text', 'string', promptstr{5}} {} { 'Style','edit'      ,'string' ,inistr{5}  'tag' 'lowfreq'}...
              { 'Style', 'text', 'string', promptstr{6}} {} { 'Style','edit'      ,'string' ,inistr{6}  'tag' 'higfreq'}...
              {}...
              { 'Style', 'text', 'string', promptstr{7}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value' str2double(inistr{7}) 'tag','rejmarks' }...
              { 'Style', 'text', 'string', promptstr{8}} {} { 'Style','checkbox'  ,'string'  ,' ' 'value' str2double(inistr{8}) 'tag' 'rejtrials'} ...
               };
    
    result = inputgui( geometry,uilist,'pophelp(''pop_rejspec'');', 'Reject by data spectra -- pop_rejspec()');
 
    size_result = size(result);
    if size_result(1) == 0
        return;
    end
    options = {};
    options = { options{:} 'elecrange' eval( [ '[' result{1} ']' ]  ) };
    options = { options{:} 'method'   lower(methodlist{result{2}})};
    thresholdsLow  = eval( [ '[' result{3} ']' ]  );
    thresholdsHigh = eval( [ '[' result{4} ']' ]  );
    options = { options{:} 'threshold' [thresholdsLow(:) thresholdsHigh(:) ] };
    freqLimitsLow   = eval( [ '[' result{5} ']' ]  );
    freqLimitsHigh  = eval( [ '[' result{6} ']' ]  );
    options = { options{:} 'freqlimits' [freqLimitsLow(:) freqLimitsHigh(:) ] };
    
    if result{7}, superpose=1; else superpose=0; end
    if result{8}, reject=1;    else reject=0; end
    options = { options{:} 'eegplotplotallrej' superpose };
    options = { options{:} 'eegplotreject'     reject };
else
    if isnumeric(varargin{3}) || (ischar(varargin{3}) && ~strcmpi(varargin{3}, 'threshold') && ~isempty(str2num(varargin{3})))
        options = {};
        if ischar(varargin{1}), varargin{1} = str2num(varargin{1}); end
        if ischar(varargin{2}), varargin{2} = str2num(varargin{2}); end
        if ischar(varargin{3}), varargin{3} = str2num(varargin{3}); end
        if ischar(varargin{4}), varargin{4} = str2num(varargin{4}); end
        if ischar(varargin{5}), varargin{5} = str2num(varargin{5}); end
        if nargin > 2, options = { options{:} 'elecrange'   varargin{1} }; end
        if nargin > 3, options = { options{:} 'threshold'   [ varargin{2}; varargin{3}]' }; end
        if nargin > 5, options = { options{:} 'freqlimits'  [ varargin{4}; varargin{5}]' }; end
        if nargin > 7, options = { options{:} 'eegplotplotallrej' varargin{6}  }; end
        if nargin > 8, options = { options{:} 'eegplotreject'     varargin{7}  }; end
        if nargin > 9, options = { options{:} 'eegplotcom'        varargin{8}  }; end
    else
        options = varargin;
    end
end

opt = finputcheck( options, { 'elecrange'     'integer'  []    [1:EEG.nbchan];
                              'threshold'     'real'     []    [-30 30];
                              'freqlimits'    'real'     []    [15 30];
                              'specdata'      'real'     []    EEG.specdata;
                              'eegplotcom'    'string'   []    '';
                              'method'        'string'   { 'fft';'multitaper' }    'multitaper';
                              'eegplotreject' 'integer'  []    0;
                              'eegplotplotallrej' 'integer'  []    0 }, 'pop_rejspec');
if ischar(opt), error(opt); end

sizewin = 2^nextpow2(EEG.pnts);
if icacomp == 1
    [allspec, Irej, tmprejE, freqs ] = spectrumthresh( EEG.data, opt.specdata, ...
							opt.elecrange, EEG.srate, opt.threshold(:,1)', opt.threshold(:,2)', opt.freqlimits(:,1)', opt.freqlimits(:,2)', opt.method);
    rejE = zeros(EEG.nbchan, EEG.trials);
    rejE(opt.elecrange,Irej) = tmprejE;
else
    % test if ICA was computed
    % ------------------------
    icaacttmp = eeg_getdatact(EEG, 'component', [1:size(EEG.icaweights,1)]);
    [allspec, Irej, tmprejE, freqs ] = spectrumthresh( icaacttmp, EEG.specicaact, ...
                            opt.elecrange, EEG.srate, opt.threshold(:,1)', opt.threshold(:,2)', opt.freqlimits(:,1)', opt.freqlimits(:,2)', opt.method);
    rejE = zeros(size(EEG.icaweights,1), size(icaacttmp,1));
    rejE(opt.elecrange,Irej) = tmprejE;
end

fprintf('%d channel selected\n', size(opt.elecrange(:), 1));
fprintf('%d/%d trials marked for rejection\n', length(Irej), EEG.trials);
rej = zeros( 1, EEG.trials);
rej(Irej) = 1;

calldisp = 0;
if nargin < 3 || opt.eegplotplotallrej == 2
    calldisp = 1;
    nbpnts = EEG.pnts;
    if icacomp == 1
        macrorej = 'EEG.reject.rejfreq';
        macrorejE = 'EEG.reject.rejfreqE';
    else			
        macrorej  = 'EEG.reject.icarejfreq';
        macrorejE = 'EEG.reject.icarejfreqE';
    end
	colrej = EEG.reject.rejfreqcol;
    
    elecrange  = opt.elecrange;
    superpose  = opt.eegplotplotallrej;
    reject     = opt.eegplotreject;
    topcommand = opt.eegplotcom;
	eeg_rejmacro; % script macro for generating command and old rejection arrays
    
	if icacomp == 1
		eegplot(EEG.data(opt.elecrange,:,:), 'srate', EEG.srate, 'winlength', 5, 'position', [100 550 800 500], ...
			'limits', [EEG.xmin EEG.xmax]*1000, 'command', command, eegplotoptions{:}); 
	else
		eegplot(icaacttmp(opt.elecrange,:,:), 'srate', EEG.srate, 'winlength', 5, 'position', [100 550 800 500], 'limits', ...
				[EEG.xmin EEG.xmax]*1000 , 'command', command, eegplotoptions{:}); 
    end
else
    if opt.eegplotreject
        EEG = pop_rejepoch(EEG, rej, 0);
    end
end
if ~isempty(rej)
	if icacomp	== 1
		EEG.reject.rejfreq = rej;
		EEG.reject.rejfreqE = rejE;
	else
		EEG.reject.icarejfreq = rej;
		EEG.reject.icarejfreqE = rejE;
	end
    Irej = find(rej);
end

% store variables
% ---------------
if icacomp == 1, EEG.specdata = allspec;
else,            EEG.specicaact = allspec;
end
    
com = [com sprintf('EEG = pop_rejspec( EEG, %s);', ...
   vararg2str({icacomp, 'elecrange', opt.elecrange,'method', opt.method, 'threshold', opt.threshold, 'freqlimits', opt.freqlimits, ...
     'eegplotcom', opt.eegplotcom, 'eegplotplotallrej' opt.eegplotplotallrej 'eegplotreject' opt.eegplotreject & ~calldisp })) ]; 

return;

% compute spectrum and reject artifacts
% -------------------------------------
function [specdata, Irej, Erej, freqs ] = spectrumthresh( data, specdata, elecrange, srate, negthresh, posthresh, startfreq, endfreq, method);

	% compute the fft if necessary - old version
    if isempty(specdata)
        if strcmpi(method, 'fft')
            sizewin = size(data,2);
            freqs = srate*[1, sizewin]/sizewin/2;
            specdata = fft( data-repmat(mean(data,2), [1 size(data,2) 1]), sizewin, 2);
            specdata = specdata( :, 2:sizewin/2+1, :);
            specdata = 10*log10(abs( specdata ).^2);
            specdata  = specdata - repmat( mean(specdata,3), [1 1 size(data,3)]);
        else
            if ~exist('pmtm')
                error('The signal processing toolbox needs to be installed');
            end
            [~, freqs] = pmtm( data(1,:,1), [],[],srate); % just to get the frequencies 	

            fprintf('Computing spectrum (using slepian tapers; done only once):\n');

            for index = 1:size(data,1)
                fprintf('%d ', index);    
                for indextrials = 1:size(data,3)
                    [ tmpspec(index,:,indextrials), freqs] = pmtm( data(index,:,indextrials) , [],[],srate);
                end
            end
            tmpspec  = 10*log(tmpspec);
            tmpspec  = tmpspec - repmat( mean(tmpspec,3), [1 1 size(data,3)]);
            specdata = tmpspec;
        end
    else
        if strcmpi(method, 'fft')
            sizewin = size(data,2);
            freqs = srate*[1, sizewin]/sizewin/2;
        else
            [~, freqs] = pmtm( data(1,:,1), [],[],srate); % just to get the frequencies 	
        end
    end
    
	% perform the rejection
	% ---------------------	
	[~, Irej, ~, Erej] = eegthresh( specdata(elecrange, :, :), size(specdata,2), 1:length(elecrange), negthresh, posthresh, ...
										 [freqs(1) freqs(end)], startfreq, min(freqs(end), endfreq));
	fprintf('\n');    
	
