function M = speech_monopitch(y, fs, min_pitch, max_pitch)

% M = speech_monopitch(y, fs, min_pitch, max_pitch)
% 
% This function computes speech monopitch series defined as the cumulative 
% sum of absolute F0 differences divided by the length of F0 computed for
% speech segments with duration of one second.
% 
% -----------------
% Input parameters:
% -----------------
% y         - input columns vector (speech signal)
% fs        - input sampling frequency [Hz]
% min_pitch - minimum pitch [Hz] (default: 75Hz)
% max_pitch - maximum pitch [Hz] (default: 400Hz)
% 
% ------------------
% Output parameters:
% ------------------
% M         - output speech monopitch 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(min_pitch)))
    min_pitch = 75;
end
if ((nargin < 4) || (isempty(max_pitch)))
    max_pitch = 400;
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 pitch
S = nan(size(Y, 2), 1);

for i = 1:size(Y, 2)
    
    % Exract the fundamental frequency (F0)
    F0 = praat_F0(Y(:, i), fs, min_pitch, max_pitch);
    
    % Trim the F0 series (by zero)
    F0 = F0(F0 ~= 0);
    
    % Process the F0 series
    if ((length(F0) > 2) && ...
        (~(sum(isnan(F0)) > 0)) && ...
        (~(sum(isinf(F0)) > 0)))
        
        % Convert to semitones
        F0 = 12*log2(F0(2:end)./F0(1:(end - 1)));
        
        % Compute cumulative sum of absolute F0 differences divided by
        % the length of F0 (normalization by the number of F0 segments)
        c = (std(F0)/mean(F0))*(cumsum(abs(diff(F0)))/length(F0));
        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.monopitch_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 monopitch');
    xlabel('segments of speech');
    ylabel('relative cumulative pitch variability');
end