function outvec = frequencySliding(data, freq_range, srate, chan)
%
% outvec = frequencySliding(data, freq_range, srate)
%
% measures instantaneous frequency by measuring the derivative of the phase
% angle timeseries 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]) 
%
%               srate =                 sampling rate of the data
%
%       OUTPUTS:
%               outvec =            matrix of calculated frequency sliding for the specified data
%                                           (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 tranform
                hilb_data = hilbert(EEG_filt);

               %%  compute frequency sliding
                for triali=1:size(data,3)
                    freqslide(:,triali) = srate*diff(unwrap(angle(hilb_data(:,triali))))/(2*pi);
                end
                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
            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 = srate*diff(unwrap(angle(hilb_data)))/(2*pi);
                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 tranform
            hilb_data = hilbert(EEG_filt);

           %%  compute frequency sliding
            for triali=1:size(data,3)
                freqslide(:,triali) = srate*diff(unwrap(angle(hilb_data(:,triali))))/(2*pi);
            end
            clear hilb_data

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

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