function P = speech_intonation(y, fs, window, overlap, padding, ...
    min_pitch, max_pitch)

% P = speech_intonation(y, fs, window, overlap, padding, ...
%    min_pitch, max_pitch)
% 
% This function computes an "intonation function" and "intonation peaks".
% Intonation function is defined as the product of absolute differences
% of frequency and intensity components computed for intervals of selected
% length (default: 0.5 seconds).
% 
% The frequency components are computed as the median pitch divided by the
% maximum pitch per segment. The energy components are computed as median
% energy divided by the maximum energy per segment.
% 
% The intonation function therefore describes how pitch and intensity do
% change in time. The peaks of the function then determines the number of
% significant changes of intonation for a given threshold (default thr.
% is set to 1.5*mean(intonation function)).
% 
% -----------------
% Input parameters:
% -----------------
% y         - input columns vector (speech signal)
% fs        - input sampling frequency [Hz]
% window    - length of window in samples or window samples
% overlap   - window overlapp in samples
% padding   - if set to 1, the last segment will be zero-padded
% min_pitch - minimum pitch [Hz] (default: 75Hz)
% max_pitch - maximum pitch [Hz] (default: 400Hz)
% 
% ------------------
% Output parameters:
% ------------------
% P         - output intonation function peaks 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(window)))
    window = hamming(1*fs); 
end
if ((nargin < 4) || (isempty(overlap)))
    overlap = floor(length(window)/2);
end
if ((nargin < 5) || (isempty(padding)))
    padding = 0;
end
if ((nargin < 3) || (isempty(min_pitch)))
    min_pitch = 75;
end
if ((nargin < 4) || (isempty(max_pitch)))
    max_pitch = 400;
end

%% Segment the signal
Y = segmentation(y, window, overlap, padding);

%% Process the segments
F_components = nan(size(Y, 2), 1);
E_components = nan(size(Y, 2), 1);

for i = 1:size(Y, 2)
    
    % Compute the pitch components
    F0 = praat_F0(Y(:, i), fs, min_pitch, max_pitch);
    F0 = F0(F0 ~= 0);
    
    nans = isnan(F0);
    infs = isinf(F0);
    
    if ((length(F0) > 2) && (~(sum(nans) > 0)) && (~(sum(infs) > 0)))
        F_components(i) = median(F0(:))/max(abs(F0(:)));
    end
    
    % Compute the energy components
    E = energy(Y(:, i), fs);
    E = E(E ~= 0);
    
    nans = isnan(E);
    infs = isinf(E);
    
    if ((length(E) > 2) && (~(sum(nans) > 0)) && (~(sum(infs) > 0)))
        E_components(i) = median(E(:))/max(abs(E(:)));
    end
end

%% Post-process the parameters
leave = (~isnan(F_components) & ~isnan(E_components));

F_components = F_components(leave);
E_components = E_components(leave);

%% Compute the intonation function series
S = (abs(diff(F_components)).*abs(diff(E_components)))/sum(leave);

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

% Set the function series
P.intonation_series = S;

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

P.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);

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

%% Plot the results
if (nargout < 1)
    figure;
    
    subplot(2, 1, 1);
    plot(y, 'b');
    title('Waveform of the input signal');
    xlabel('samples');
    ylabel('amplitude');
    
    subplot(2, 1, 2);
    findpeaks(S, fs, 'MinPeakHeight', 1.5*mean(S));
    title('Intonation function');
    xlabel('samples');
    ylabel('intonation function (peaks)');
end