% Apply the springer segmentation onto one audio file.
% path_to_wav = Path to the audio file.
% output_file = Optionally write the result to this .mat file.
% Returns: A list like [[s1, s2, sys, dia], [s1, s2, sys, dia], ...] .
function [heartbeat_list, segment_lengths, summary_features] = do_single_segmentation(path_to_wav, output_file)

    %% Load the trained parameter matrices for Springer's HSMM model.
    % The parameters were trained using 409 heart sounds from MIT heart
    % sound database, i.e., recordings a0001-a0409.
    load('Springer_B_matrix.mat');
    load('Springer_pi_vector.mat');
    load('Springer_total_obs_distribution.mat');
    
    if(~strcmp(path_to_wav(max(1,(end-3)):end),'.wav'))
        path_to_wav = [path_to_wav '.wav'];
    end

    %% Load data and resample data
    springer_options   = default_Springer_HSMM_options;
    [PCG, Fs1]         = audioread(path_to_wav);  % load data
    
    psd_win_len = 256;
    psd_stride = 1;
    
    % crop PCG to max. 200 seconds
    crop = Fs1*200;
    
    if crop < length(PCG)
        fprintf('crop to %d samples\n', crop);
        PCG = PCG(1:crop);
    end
    
    PCG_resampled      = resample(PCG,springer_options.audio_Fs,Fs1); % resample to springer_options.audio_Fs (1000 Hz)

    %% Running runSpringerSegmentationAlgorithm.m to obtain the assigned_states
    assigned_states = runSpringerSegmentationAlgorithm(PCG_resampled, springer_options.audio_Fs, Springer_B_matrix, Springer_pi_vector, Springer_total_obs_distribution, false); % obtain the locations for S1, systole, s2 and diastole

    tokenized_states = round(tokenize_states(assigned_states));
    num_beats = size(tokenized_states,1);

    heartbeat_list = cell(0); % This is a list containing lists of segments.
    segment_lengths = zeros(num_beats-1,5);

    for k = 1:num_beats-1

        % Split the heartbeat into its segments.
        seg1 = PCG_resampled(tokenized_states(k,1) : (tokenized_states(k,2)-1))';
        seg2 = PCG_resampled(tokenized_states(k,2) : (tokenized_states(k,3)-1))';
        seg3 = PCG_resampled(tokenized_states(k,3) : (tokenized_states(k,4)-1))';
        seg4 = PCG_resampled(tokenized_states(k,4) : (tokenized_states(k+1,1)-1))';
        
        segments = {seg1 seg2 seg3 seg4};
        heartbeat_list{k} = segments;
        
        lengths = [length(seg1) length(seg2) length(seg3) length(seg4)];
        lengths = [lengths sum(lengths)];
        segment_lengths(k,:) = lengths;
    end
    
    summary_features = extractFeaturesFromHsIntervals(assigned_states, PCG_resampled)';
    
    pcg_psd = pwelch(PCG_resampled,psd_win_len,[],psd_win_len);
%   pcg_psd = pwelch(PCG_resampled-mean(PCG_resampled),psd_win_len,[],psd_win_len);
    pcg_psd = pcg_psd(1:psd_stride:end);
%   pcg_psd = pcg_psd/sum(pcg_psd);
    
   summary_features = [summary_features; pcg_psd];
    
    if(nargin > 1)
        % Write heartbeat_list to file.
        patientlist = {heartbeat_list};
        summary_features = {summary_features};
        segment_lengths = {segment_lengths};
        save(output_file, 'patientlist', 'segment_lengths', 'summary_features');
    end
