function M = speech_monoloudness(y, fs)

% M = speech_monoloudness(y, fs)
% 
% This function computes speech monopitch series defined as the cumulative 
% sum of absolute E differences divided by the length of E computed for
% speech segments with duration of one second.
% 
% -----------------
% Input parameters:
% -----------------
% y         - input columns vector (speech signal)
% fs        - input sampling frequency [Hz]
% 
% ------------------
% Output parameters:
% ------------------
% M         - output speech monoloudness structure
%
%
%
% --
% ing. Zoltán Galáž
% xgalaz00@stud.feec.vutbr.cz       
% 
% Department of Telecommunications
% Faculty of Electrical Engineering and Communication
% Brno University of Technology

%% Normalize the signal
if (any(abs(y) > 1))
    y = y/max(abs(y(:)));
end

%% Segment the signal
window  = hamming(1.0*fs);
overlap = floor(length(window)/2.0);
padding = 0;

Y = segmentation(y, window, overlap, padding);

%% Compute the monotony of energy
S = nan(size(Y, 2), 1);

for i = 1:size(Y, 2)
    
    % Exract the short-time energy
    E_window  = hamming(0.02*fs);
    E_overlap = floor(length(E_window)/2.0);
    E_padding = 0;
    
    E = energy(Y(:, i), fs, E_window, E_overlap, E_padding);
    
    % Trim the E series (by zero)
    E = E(E ~= 0);
    
    % Process the E series
    if ((length(E) > 2) && ...
        (~(sum(isnan(E)) > 0)) && ...
        (~(sum(isinf(E)) > 0)))
        
        % Normalize the energy
        E = E/max(E(:));
        
        % Compute cumulative sum of absolute E differences divided by
        % the length of E (normalization by the number of E segments)
        c = (std(E)/mean(E))*(cumsum(abs(diff(E)))/length(E));
        x = (0:(length(c) - 1));
        
        c = c(:);
        x = x(:);
        
        m_c = mean(c);
        m_x = mean(x);
        
        % Compute the slope
        S(i) = sum((x - m_x).*(c - m_c))/sum((x - m_x).^2);
    end
end

%% Trim the NaN values
S = S(~isnan(S));

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

% Set the function series
M.monoloudness_series = S;

% 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);

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

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

    plot(1:length(S), S);
    title('Speech monoloudness');
    xlabel('segments of speech');
    ylabel('relative cumulative energy variability');
end