function [seg_index, s1_end, s2_beg, s2_end, beat_len, assigned_states] = pcg_segmentation(pcg, fsig, fs, method, doPlot)
    %THRESHOLD_BASED_SEGMENTATION Summary of this function goes here
    %   Detailed explanation goes here
    % generate template signal
    if strcmp(method,'simple')
        % calculate cross-correlation
        fsig = fsig/max(fsig);
        csig = xcorr(fsig,'coeff');

        % select correct window
        tot_len = 2*ceil(60/50*fs);
        csig_win = csig((length(fsig)- tot_len):(length(fsig)+tot_len));

        % search highest max
        d_csig_win = diff(csig_win);
        d_csig_gz = d_csig_win > 0;
        zero_cross_ind = find(d_csig_gz(1:end-1) ~= d_csig_gz(2:end));
        [~, ind] = sort(csig_win(zero_cross_ind),'descend');
        max_ind = zero_cross_ind(ind(1:2));
        max_ind = sort(max_ind);
        beat_len = floor(mean(diff(max_ind)));

        % segment signal
        % global_max_ind = find(fsig == max(fsig));
        % start_ind = global_max_ind - floor(beat_len/8);
        % seg_index = (mod(start_ind,beat_len):beat_len:(length(fsig)-beat_len+1));
        % seg_start_ind = find(seg_index == start_ind);
        % [~, seg_index] = findpeaks(fsig, 'MinPeakDistance', 0.9*beat_len, 'MinPeakHeight', 0.5*max(fsig));
        [~, seg_index] = findpeaks(fsig, 'MinPeakDistance', 0.8*beat_len);
        seg_index = seg_index  - floor(0.1*beat_len);
        
        s1_end = [];
        s2_end = [];
        s2_beg = [];
        assigned_states = [];

        if doPlot
            figure;
            clf;
            hold on;
            plot(csig);
            plot((length(fsig)- tot_len):(length(fsig)+tot_len),csig_win,'y-');
            plot(max_ind+(length(fsig)- tot_len),csig_win(max_ind),'r*','Linewidth',3);
            hold off;
        end
    elseif strcmp(method, 'springer')
%         load('example_model.mat');
        load('Schmidt_B_matrix.mat');
        load('Schmidt_pi_vector.mat');

        assigned_states = runSchmidtSegmentationAlgorithm(pcg, fs, Schmidt_B_matrix, Schmidt_pi_vector, doPlot);
     
        seg_index = find(assigned_states == 4 & [assigned_states(2:end);0] == 1);
        s1_end = find(assigned_states == 1 & [assigned_states(2:end);0] == 2);
        s2_beg = find(assigned_states == 2 & [assigned_states(2:end);0] == 3) + 1;
        s2_end = find(assigned_states == 3 & [assigned_states(2:end);0] == 4);
        
        s1_end(s1_end < seg_index(1)) = [];
        s2_beg(s2_beg < seg_index(1)) = [];
        s2_end(s2_end < seg_index(1)) = [];
        s1_end(s1_end > seg_index(end)) = [];
        s2_beg(s2_beg > seg_index(end)) = [];
        s2_end(s2_end > seg_index(end)) = [];
        
        beat_len = mean(abs(diff(seg_index)));
    end
end

