function V = energy_variability(E, fs, decibels)

% V = energy_variability(E, decibels)
% 
% This function computes energy variability, which is defined as the series
% of absolute differences between consecutive E values (the values define
% E for segments) that do not equal to zero, trimmed by the median of the
% series of the differences.
% 
% This function has the optional parameter "decibels", which can be used
% to convert the zero-trimmed E series to decibels scale.
% 
% -----------------
% Input parameters:
% -----------------
% E         - input energy series
% fs        - sampling frequency [Hz]
% decibels  - conversion of amplitude to decibels (default: false)
% 
% ------------------
% Output parameters:
% ------------------
% V         - output pitch variability structure
%
%
%
% --
% ing. Zoltán Galáž
% xgalaz00@stud.feec.vutbr.cz       
% 
% Department of Telecommunications
% Faculty of Electrical Engineering and Communication
% Brno University of Technology

%% Paths and variables
if ((nargin < 3) || (isempty(decibels)))
    decibels = false;
end

%% Trim the E series (by zero)
E = E(E ~= 0);

%% Convert to decibels (default: false)
if (decibels)
    E = mag2db(E);
end

%% Compute the absolute differences
D = abs(diff(E));

%% Compute the energy variability series
S = D(D >= median(D));

%% Compute the output features
V = struct();

% Set the function series
V.variability_series = S;

% Set the peak features
[peak_values, ~, ~, ~] = findpeaks(S, fs, 'MinPeakHeight', 2.5*iqr(S));

V.num_peaks = length(peak_values);

% Set the slope of cumulative function
amps = S/max(S);
cums = cumsum(amps)/length(amps);
samp = 0:(length(cums) - 1);

cums = cums(:);
samp = samp(:);

m_c = mean(cums);
m_x = mean(samp);

V.cumsum_slope = sum((samp - m_x).*(cums - m_c))/sum((samp - m_x).^2);

%% Plot the energy variability series (if nargout < 1)
if (nargout < 1)
    figure;

    plot(1:length(S), S);
    title('Energy variability');
    xlabel('segments');
    
    if (semitones)
        ylabel('decibels');
    else
        ylabel('energy');
    end
end