function [frm_out, F1, F2, F3] = praat_formants(wav_pth, dir_pth, time_step, max_frm, win_len, preem)
%% Funkce pro vypocet formantu F1, F2, F3 a jejich pasem B1, B2, B3 za pomoci programu Praat
% frm_out = praat_formants(wav_pth, time_step, max_frm, win_len, preem)
% 
% This function extracts first three formants and their bandwidths using
% Praat. The files praatcon.exe and praat_formants.praat must be in the
% same directory.
% 
% wav_pth       - path to the *.wav file
% dir_pth       - path to the directory with Praat files with backslash at
%                 the end
% time_step     - time step in seconds (default: 0)
% max_frm       - maximum formant in Hz (default: 5500)
% win_len       - window length in s (default: 0.025)
% preem         - pre-emphasis from Hz (default: 50)
% 
% Fi represents the i-th formant [Hz]; i = {1, 2, 3}
% Fib represents the i-th formant bandwidth [Hz]; i = {1, 2, 3}
% 
% frm_out.Fi_mean = mean(Fi);
% frm_out.Fi_var = var(Fi);
% frm_out.Fi_std = std(Fi);
% frm_out.Fi_max = max(Fi);
% frm_out.Fi_min = min(Fi);
% frm_out.Fi_med = median(Fi);
% 
% frm_out.Fib_mean = mean(Fib);
% frm_out.Fib_var = var(Fib);
% frm_out.Fib_std = std(Fib);
% frm_out.Fib_max = max(Fib);
% frm_out.Fib_min = min(Fib);
% frm_out.Fib_med = median(Fib);

%% Paths and variables
if((nargin < 2) || isempty(dir_pth))
    dir_pth = '';
end
if((nargin < 3) || isempty(time_step))
    time_step = 0;
end
if((nargin < 4) || isempty(max_frm))
    max_frm = 5500;
end
if((nargin < 5) || isempty(win_len))
    win_len = 0.025;
end
if((nargin < 6) || isempty(preem))
    preem = 50;
end

%% Prepare and run Praat command
tmp_pth = [pwd '\' num2str(round(rand(1,1)*1e5)) '.tmp'];
while(exist(tmp_pth, 'file'))
    tmp_pth = [pwd '\' num2str(round(rand(1,1)*1e5)) '.tmp'];
end

% tmp_pth = 'a4.txt';
% cesta_praat = 'D:\Matlab\Funkce\';
% cesta = [cesta_praat 'praatcon.exe ' cesta_praat 'praat_formants.praat'];
% cmd = [cesta ' ' wav_pth ' ' tmp_pth ' ' time_step ' 3 ' max_frm ' ' win_len ' ' preem];
% stat = system(cmd);

cmd = [dir_pth 'praatcon.exe' ' ' dir_pth 'praat_formants.praat ' ...
    wav_pth ' ' tmp_pth ' ' num2str(time_step) ' 3 ' num2str(max_frm) ...
    ' ' num2str(win_len) ' ' num2str(preem)];
stat = system(cmd);

%% Load data to the environment
if(~stat)
    fid = fopen(tmp_pth, 'r');
    frm = textscan(fid, '%f %f %f %f %f %f', 'HeaderLines', 1, ...
        'treatAsEmpty', '--undefined--');
    fclose(fid);

    frm{1,1}(~(frm{1,1} > 0)) = [];
    frm{1,2}(~(frm{1,2} > 0)) = [];
    frm{1,3}(~(frm{1,3} > 0)) = [];
    frm{1,4}(~(frm{1,4} > 0)) = [];
    frm{1,5}(~(frm{1,5} > 0)) = [];
    frm{1,6}(~(frm{1,6} > 0)) = [];

    delete(tmp_pth);
    
    % Extract the formant info
    F1 = frm{1,1};
    frm_out.F1_mean = mean(frm{1,1});
    frm_out.F1_var = var(frm{1,1});
    frm_out.F1_std = std(frm{1,1});
    frm_out.F1_max = max(frm{1,1});
    frm_out.F1_min = min(frm{1,1});
    frm_out.F1_med = median(frm{1,1});
    frm_out.F1_1p = prctile(frm{1,1},1);
    frm_out.F1_99p = prctile(frm{1,1},99);
    frm_out.F1_ir = prctile(frm{1,1},99) - prctile(frm{1,1},1);
    
    frm_out.F1b_mean = mean(frm{1,2});
    frm_out.F1b_var = var(frm{1,2});
    frm_out.F1b_std = std(frm{1,2});
    frm_out.F1b_max = max(frm{1,2});
    frm_out.F1b_min = min(frm{1,2});
    frm_out.F1b_med = median(frm{1,2});
    frm_out.F1b_1p = prctile(frm{1,2},1);
    frm_out.F1b_99p = prctile(frm{1,2},99);
    frm_out.F1b_ir = prctile(frm{1,2},99) - prctile(frm{1,2},1);
    
    F2 = frm{1,3};
    frm_out.F2_mean = mean(frm{1,3});
    frm_out.F2_var = var(frm{1,3});
    frm_out.F2_std = std(frm{1,3});
    frm_out.F2_max = max(frm{1,3});
    frm_out.F2_min = min(frm{1,3});
    frm_out.F2_med = median(frm{1,3});
    frm_out.F2_1p = prctile(frm{1,3},1);
    frm_out.F2_99p = prctile(frm{1,3},99);
    frm_out.F2_ir = prctile(frm{1,3},99) - prctile(frm{1,3},1);
    
    frm_out.F2b_mean = mean(frm{1,4});
    frm_out.F2b_var = var(frm{1,4});
    frm_out.F2b_std = std(frm{1,4});
    frm_out.F2b_max = max(frm{1,4});
    frm_out.F2b_min = min(frm{1,4});
    frm_out.F2b_med = median(frm{1,4});
    frm_out.F2b_1p = prctile(frm{1,4},1);
    frm_out.F2b_99p = prctile(frm{1,4},99);
    frm_out.F2b_ir = prctile(frm{1,4},99) - prctile(frm{1,4},1);
    
    F3 = frm{1,5};
    frm_out.F3_mean = mean(frm{1,5});
    frm_out.F3_var = var(frm{1,5});
    frm_out.F3_std = std(frm{1,5});
    frm_out.F3_max = max(frm{1,5});
    frm_out.F3_min = min(frm{1,5});
    frm_out.F3_med = median(frm{1,5});
    frm_out.F3_1p = prctile(frm{1,5},1);
    frm_out.F3_99p = prctile(frm{1,5},99);
    frm_out.F3_ir = prctile(frm{1,5},99) - prctile(frm{1,5},1);
    
    frm_out.F3b_mean = mean(frm{1,6});
    frm_out.F3b_var = var(frm{1,6});
    frm_out.F3b_std = std(frm{1,6});
    frm_out.F3b_max = max(frm{1,6});
    frm_out.F3b_min = min(frm{1,6});
    frm_out.F3b_med = median(frm{1,6});
    frm_out.F3b_1p = prctile(frm{1,6},1);
    frm_out.F3b_99p = prctile(frm{1,6},99);
    frm_out.F3b_ir = prctile(frm{1,6},99) - prctile(frm{1,6},1);
end

