function outvec = frequencyPower(data, freq_range, srate, chan)
%
% outvec = frequencyPower(data, freq_range, srate)
%
% measures instantaneous power by measuring square of the magnitude of the 
% analytic signal according to methods of Cohen, M, J Neuro, 2014
% 
%       INPUTS: 
%               data =                 EEG data in EEGLAB format 
%                                           (channels x timepoints x trials)
%
%               freq_range =       vector of two integers (ex: [7 13])
%                                           buffer of 15% for lower and
%                                           upper bounds will be added for
%                                           edge effects
%
%               srate =                 sampling rate of the data
%
%       OUTPUTS:
%               outvec =            matrix of calculated instantaneous
%                                           power for the specified data with the mean taken over all the trials
%                                           (channels x timepoints x trials) 
%
% written by annemarie wolff, awolf037@uottawa.ca
%%   
buffer_hz = [freq_range(1)*.15 freq_range(2)*.15];

%% compute function
% tic;
if ~isempty(data)
    if nargin<4
        for chani=1:size(data,1)
            if length(size(data))==3

                    % bandpass filter the data
                    EEG_filt = bandpass(squeeze(data(chani,:,:)), [(freq_range(1)-buffer_hz(1)) (freq_range(2)+buffer_hz(2))], srate);

                   %% take the hilbert transform
                    hilb_data = hilbert(EEG_filt);

                   %%  compute power sliding
                    for triali=1:size(data,3)
                        freqslide(:,triali) = (abs(hilb_data(:,triali))).^2;
                    end
                    clear hilb_data

                   %% apply a median filter, then take the mean over all the trials
            %         outvec(chani,:,:) = mean(freqslide,2);
                    outvec(chani,:,:) = medfilt1(freqslide, 10, [], 1,'omitnan');
                    clear freqslide
            else % only 1 trial, chan * timepoints
                EEG_filt = bandpass(squeeze(data(chani,:)), [(freq_range(1)-buffer_hz(1)) (freq_range(2)+buffer_hz(2))], srate);

               %% take the hilbert tranform
                hilb_data = hilbert(EEG_filt);
               %%  compute frequency sliding
                freqslide = (abs(hilb_data)).^2;
                clear hilb_data

               %% apply a median filter
        %         outvec(chani,:) = mean(medfilt1(freqslide, 10, [], 2,'omitnan'),2);
                outvec(chani,:) = medfilt1(freqslide, 10, [], 1,'omitnan');
                clear freqslide
            end
        end
    else
        chani = chan;
        % bandpass filter the data
        EEG_filt = bandpass(squeeze(data(chani,:,:)), [(freq_range(1)-buffer_hz(1)) (freq_range(2)+buffer_hz(2))], srate);

        %% take the hilbert transform
        hilb_data = hilbert(EEG_filt);

        %%  compute power sliding
        for triali=1:size(data,3)
        freqslide(:,triali) = (abs(hilb_data(:,triali))).^2;
        end
        clear hilb_data

        %% apply a median filter, then take the mean over all the trials
        %         outvec(chani,:,:) = mean(freqslide,2);
        outvec(:,:) = medfilt1(freqslide, 10, [], 1,'omitnan');
        clear freqslide

    end
else
    
   outvec = []; 
    
end
% elapsed_time = toc;
% disp([ 'File completed frequency power in ' num2str(round(elapsed_time,1)) 's.']);
