function V = pitch_variability(F0, fs, semitones)

% V = pitch_variability(F0, semitones)
%
% This function computes pitch variability, which is defined as the series
% of absolute differences between consecutive F0 values (the values define
% F0 for segments) that do not equal to zero, trimmed by the median of the
% series of the differences.
% 
% This function has the optional parameter "semitones", which can be used
% to convert the zero-trimmed F0 series to semitone scale. As repotred by
% Little et al. (pitch period entropy, DOI: 10.1109/TBME.2008.2005954) the
% conversion makes the features independent of the speakers' gender.
% 
% -----------------
% Input parameters:
% -----------------
% F0        - input F0 series [Hz]
% fs        - sampling frequency [Hz]
% semitones - conversion of F0 to semitones (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(semitones)))
    semitones = false;
end

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

%% Convert to semitones (default: false)
if (semitones)
    F0 = 12*log2(F0(2:end)./F0(1:(end - 1)));
end

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

%% Compute the pitch 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', 1.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 pitch variability series (if nargout < 1)
if (nargout < 1)
    figure;

    plot(1:length(S), S);
    title('Pitch variability');
    xlabel('segments');
    
    if (semitones)
        ylabel('semitones');
    else
        ylabel('frequency [Hz]');
    end
end