function syllables_positions = BACD_alg(speechFile, fs)
% BACD_alg.m
% Input: speechFile - analysed speech file
%        fs - its sampling frequency
% Output: sample positions of detected syllables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Master thesis
% winter semester 2019/20
% Vojtech Illner
%
% Algorithm for syllables detection based on BACD approach [1]
% and pitch verification (via SWIPE method) [2].
% version 1 (22. 5. 2020)
% 
% References:
% [1] R. Cmejla, J. Rusz, P. Bergl, J. Vokral, Bayesian changepoint detection
% for the automatic assessment of fluency and articulatory disorders. Spee.
% Comm., 55(1) (2013), 178-189.
% [2] A. Camacho, G. Harris, A sawtooth waveform inspired pitch estimator
% for speech and music. J. Ac. Soc. Am. 124 (2008), 1638.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fs_re = 16000; % [Hz]
% Thresholds
% Minimal length of continous pitch segment to count it as a syllable.
thr.f0_int_len = 0.025*fs_re; % [samples, fs = 16 kHz]
% Minimal distance from BACD peak to count F0 local minima as syllable
% boundary.
thr.f0_min_dist = 0.125*fs_re; % [samples, fs = 16 kHz]
% Minimal interval length to allow counting of F0 local minimas.
thr.int_len = thr.f0_int_len; % [samples, fs = 16 kHz]
% Minimal distance between the syllables
thr.minSyllDist = 0.09*fs_re; %[samples, fs = 16 kHz]

% Analyzed length -testing puposes
len = size(speechFile, 1); % [samples]
signal_speech = speechFile(1 : len);

% Resample to fs = 16 000 Hz
signal_speech = resample(signal_speech, fs_re, fs);

% equalize the frame signal to match BACD requirements (input cant be of too low amplitude)
low_output_threshold = 10^(-6);
low_output_map = (abs(signal_speech) < low_output_threshold);
signal_speech(low_output_map) = [];

% Segment the signal into frames of length frame_len
frame_len = fs_re * 2; % [samples]
frame_step = frame_len * 0.95; % [samples] - calculated as percentage of the frame length
N = ceil((length(signal_speech)-frame_len) / frame_step) + 1; % number of frames
peaks = cell(1, N);
locs = cell(1, N);
syllables_positions = zeros(N*15, 1);
syll_count = 0;

for j = 1 : N
    syll_count_temp = 0;
    syll_pos_temp = zeros(15, 1);
    
    frame_start = (j-1)*frame_step +1;
    % extract frame
    if j == N % situation of the last frame (might be different length)
        signal_speech_frame = signal_speech(frame_start : end);
    else
        signal_speech_frame = signal_speech(frame_start : (j-1)*frame_step + frame_len);
    end
    
    % Run the BACD algorithm [1]
    win_length_bacd = 0.03 * fs_re; % [samples]
    bacd_out = RBD(signal_speech_frame, win_length_bacd);

    % Lowpass the BACD output
    break_freq = 20; % break frequency [Hz]
    bacd_out_filt = lowpass(bacd_out, break_freq, fs_re);

    % Smooth the lowpassed BACD output using its envelope
    win_length_envelope = 0.025*fs_re; % [samples]
    [env_upper, ~] = envelope(bacd_out_filt, win_length_envelope, 'peak');

    % Locate the peaks in the envelope
    MinPeakDistance = 0.075*fs_re; % [samples]
    MinPeakProminence = 0.4; % minimal height from either side
    [peaks{j}, locs{j}] = findpeaks(env_upper, 'MinPeakDistance', MinPeakDistance, 'MinPeakProminence', MinPeakProminence);
    
    % Perform pitch analysis
    Fmin = 60;
    Fmax = 400;
    STHR = 0.22; % def 0.22
    timestep = 0.01;
    % Run the SWIPE algorithm [2]
    [SWIPE.f0,SWIPE.time] = swipep(signal_speech_frame, fs_re, [Fmin Fmax], timestep, [], 1/20, 0.5, STHR);
    % locate significant local minimas in F0 contour
    MinPeakDistance_F0 = 10; % [samples]
    MinPeakProminence_F0 = 5;
    [~, locs1_f0] = findpeaks(-SWIPE.f0, 'MinPeakDistance', MinPeakDistance_F0, 'MinPeakProminence', MinPeakProminence_F0);
    [~, locs2_f0] = findpeaks(SWIPE.f0, 'MinPeakDistance', MinPeakDistance_F0, 'MinPeakProminence', MinPeakProminence_F0);
    locs_f0 = sort([locs1_f0; locs2_f0]);
    
    % Localize the intervals
    if isempty(peaks(j)) % number of intervals
        N_intervals = 0;
    else
        N_intervals = size(peaks{j}, 2);
    end
    locs{j} = [locs{j}(:)' length(signal_speech_frame)]; % add end of the frame as interval boundary
    % Define the intervals bounds and pitch contour inside the interval
    if N_intervals ~= 0
        intervals_bounds = zeros(N_intervals, 2);
        for iii = 1 : N_intervals
            intervals_bounds(iii, 1) = locs{j}(iii);
            intervals_bounds(iii, 2) = locs{j}(iii + 1);

            % pitch boundary and contour
            [~,id_start] = min(abs(SWIPE.time*fs_re - intervals_bounds(iii, 1)));
            [~,id_end] = min(abs(SWIPE.time*fs_re - intervals_bounds(iii, 2)));
            interval.f0 = SWIPE.f0(id_start:id_end);
            interval.f0_time = SWIPE.time(id_start:id_end);

            % evaluate pitch segments
            idx_f0seg = all(isnan(interval.f0),2);
            idr_f0seg = diff(find([1;diff(idx_f0seg);1]));
            f0_seg = mat2cell(interval.f0,idr_f0seg(:),size(interval.f0,2));
            f0_seg(cellfun(@(f0_seg) any(isnan(f0_seg)),f0_seg)) = [];
            interval.f0_segments = f0_seg;

            % length of the segments
            if ~isempty(interval.f0_segments)
                interval.f0_segments_length = zeros(size(interval.f0_segments, 1), 1);
                for iiii = 1 : size(interval.f0_segments, 1)
                    interval.f0_segments_length(iiii) = size(interval.f0_segments{iiii}, 1) * timestep * fs_re; % [samples]
                end
            else
                interval.f0_segments_length = 0; 
            end

            % localize the presence of potencial F0 local minimas in the interval
            if ~isempty(locs_f0(locs_f0 >= id_start & locs_f0 <= id_end))
                interval.f0_minimas = locs_f0;
                interval.f0_minimas = interval.f0_minimas(interval.f0_minimas >= id_start & interval.f0_minimas <= id_end);
                interval.f0_minimas = SWIPE.time(interval.f0_minimas) * fs_re;

                % calculate its distance from closest interval boundary
                if ~isempty(interval.f0_minimas)
                    interval.f0_minimas_distFromBoundary = zeros(size(interval.f0_minimas));
                    for iiiii = 1 : size(interval.f0_minimas, 1)
                        interval.f0_minimas_distFromBoundary(iiiii) = min([abs(interval.f0_minimas(iiiii) - intervals_bounds(iii,1));...
                            abs(interval.f0_minimas(iiiii) - intervals_bounds(iii,2))]);
                        interval.f0_minimas_distFromBoundary(iiiii) = unique(interval.f0_minimas_distFromBoundary(iiiii)); % case where the location is precisely in the middle of the interval
                    end
                else
                    interval.f0_minimas_distFromBoundary = 0;
                end
            else
                interval.f0_minimas = 0;
                interval.f0_minimas_distFromBoundary = 0;
            end

            % Count the number of syllables in the segment based on the
            % info above and given thresholds contained in thr struct.
            % Number of pitch intervals longer then f0_int_len - N1
            interval.N_syll = size(interval.f0_segments_length(interval.f0_segments_length > thr.f0_int_len), 1);
            interval.N_syll_pos = zeros(interval.N_syll, 1);
            for iiiiii = 1 : interval.N_syll
                pos_times = interval.f0_time(interval.f0 == interval.f0_segments{iiiiii}(1)) * fs_re;
                interval.N_syll_pos(iiiiii) = pos_times(1);
            end
            % Add potential influency of f0 local minima
            if (intervals_bounds(iii,2) - intervals_bounds(iii,1)) > thr.int_len
                  N_syll2 = size(interval.f0_minimas_distFromBoundary(interval.f0_minimas_distFromBoundary > thr.f0_min_dist), 1);
                  interval.f0_minimas = interval.f0_minimas(interval.f0_minimas_distFromBoundary > thr.f0_min_dist);
                interval.N_syll = interval.N_syll + N_syll2;
                  for iiiiiii = 1 : N_syll2
                      interval.N_syll_pos = [interval.N_syll_pos; interval.f0_minimas(iiiiiii)];
                  end
            end
           if interval.N_syll ~= 0
               for iiiiiiii = 1 : interval.N_syll
                   syll_count_temp = syll_count_temp +1;
                   syll_pos_temp(syll_count_temp) = frame_start + interval.N_syll_pos(iiiiiiii);
               end
           end
        end
    end
    % Check distances between estimated syllables
    syll_pos_temp(syll_pos_temp == 0) = [];
    if syll_count_temp ~= 0
        for i9 = 1 : syll_count_temp
            if i9 == 1
                syll_count = syll_count +1;
                syllables_positions(syll_count) = syll_pos_temp(i9);
            elseif syll_pos_temp(i9) > syll_pos_temp(i9-1) + thr.minSyllDist
                syll_count = syll_count +1;
                syllables_positions(syll_count) = syll_pos_temp(i9);
            end
            syllables_positions(syll_count) = syllables_positions(syll_count) + sum(low_output_map(1:syllables_positions(syll_count))); % shift because of primal signal equalization for BACD
        end
    end
end
syllables_positions(syllables_positions == 0) = [];
syllables_positions = round(syllables_positions * (fs/fs_re));
end




