% Written by Péter Hajas and Márton Áron Goda
% References:
% Adaptive neuro-fuzzy inference system for diagnosis of the heart valve diseases using wavelet transform with entropy (Harun Uguz)
% Support Vector Machine-based identification of heart valve diseases using heart sounds (Ilias Maglogiannis et al)

function [wavelet_morphology_features, wavelet_column_indices, kurtosis_features, kurtosis_column_indices] = get_wavelet_morphology_vector_v04(audio_data, fs, wavelet_type, decomp_levels, states_starts, envelope_type, to_normalize, n_hs_portions, hs_corr, n_hs_kurt, figures)

S1_start = states_starts(:, 1);
dia_end  = states_starts(:, 5);

[n_heartcycles, ~] = size(states_starts);

% Daubechies dekompozíciók
[decomps, decomp_lengths] = wavedec(audio_data, decomp_levels, wavelet_type);

decomp_starts = cumsum([1; decomp_lengths]);
decomp_starts = decomp_starts(1:end-2);
decomp_lengths= decomp_lengths(1:end-1);

n_hs_portions = n_hs_portions .* 2;
hs_start_corr = hs_corr(1, :);
hs_stop_corr  = hs_corr(2, :);

wavelet_morphology_features = zeros(1, sum(n_hs_portions) * sum(1./pow2([1:decomp_levels, decomp_levels])));
n_kur_feat = sum(sum(n_hs_kurt));
kurtosis_features = zeros(1, length(decomp_starts) * n_kur_feat);

wavelet_column_indices = sum(n_hs_portions) .* fliplr((1./pow2([1:decomp_levels, decomp_levels])));
wavelet_column_indices = cumsum([1, wavelet_column_indices]);

kur_col_lvl_offset = cumsum([1, sum(n_hs_kurt(1,:))]);
kurtosis_column_indices = [kur_col_lvl_offset, n_kur_feat+kur_col_lvl_offset, 2*n_kur_feat+kur_col_lvl_offset, 3*n_kur_feat+kur_col_lvl_offset];
kurtosis_column_indices = [kurtosis_column_indices, 4*n_kur_feat + 1];

plotting = cell(2, decomp_levels);

if(figures)
	figure(55);
	clf;
end

for k = 1 : length(decomp_starts)
	
	cD = decomps(decomp_starts(k) : decomp_starts(k) + decomp_lengths(k) - 1);
	
	decomp_div = length(audio_data) / length(cD);
	decomp_div = round(decomp_div);
	n_excess = length(cD) - length(audio_data) / decomp_div;
	cD = cD( 1 + floor(n_excess/2) : length(cD) - ceil(n_excess/2));
	
	fs_decomp = fs / decomp_div;
	
	if(figures)
		plotting{1, k} = cD;
		plotting{2, k} = linspace(1, length(audio_data), length(cD));
	end
	
	[max_vector_of_heartcycles, std_vector_of_heartcycles] = get_max_std_vector_of_heartcycles(cD, S1_start/decomp_div, dia_end/decomp_div);
	
	wavelet_features_of_level = zeros(1, sum(n_hs_portions) / decomp_div);
	kurtosis_features_of_level = zeros(1, n_kur_feat);
	for i = 1 : 4
		n_wave_parts = round(n_hs_portions(i) / decomp_div);
		n_kurt_parts = n_hs_kurt(:, i)';
		hs_start = max(states_starts(:, i)   / decomp_div + fs_decomp * hs_start_corr(i) * [0;ones(n_heartcycles-1, 1)], 1);
		hs_stop  = max(states_starts(:, i+1) / decomp_div + fs_decomp * hs_stop_corr(i)  * [0;ones(n_heartcycles-1, 1)], 1);
		hsf_stop = sum(n_hs_portions(1:i)) / decomp_div;
		[wavelet_level_features_of_segment, kurtosis_features_of_segment] = get_level_features_of_segment(cD, fs_decomp, hs_start, hs_stop, n_wave_parts, n_kurt_parts, max_vector_of_heartcycles, std_vector_of_heartcycles, envelope_type, to_normalize);
		wavelet_features_of_level(1, hsf_stop - n_wave_parts + 1 : hsf_stop) = wavelet_level_features_of_segment;
		kurtosis_features_of_level(1, i) = kurtosis_features_of_segment(1);
		position_shift = [0, n_hs_kurt(2,:)];
		kurt_strt = 5 + sum(position_shift(1:i));
		kurtosis_features_of_level(1, kurt_strt : kurt_strt+n_hs_kurt(2,i)-1) = kurtosis_features_of_segment(2:end);
	end
	lkl = length(kurtosis_features_of_level);
	kurtosis_features(1, (k-1)*lkl+1 : k*lkl) = kurtosis_features_of_level;
	wavelet_morphology_features(wavelet_column_indices(k):wavelet_column_indices(k+1)-1) = wavelet_features_of_level;
	
end

if(figures)
	for k = 1 : length(decomp_starts)
		y = plotting{1, k};
		x = plotting{2, k};
		figure(55);
		hold on;
		plot(x, y);
	end
end

%% Intervallum maximuma
function [max_vector_of_heartcycles, std_vector_of_heartcycles] = get_max_std_vector_of_heartcycles(audio_data, heartcycle_starts, heartcycle_ends)
heartcycle_starts = floor(heartcycle_starts(:)') + 1;
heartcycle_ends = floor(heartcycle_ends(:)');
max_vector_of_heartcycles = zeros(1, length(heartcycle_starts));
std_vector_of_heartcycles = zeros(1, length(heartcycle_starts));
for k = 1 : length(heartcycle_starts)
	audio_portion = audio_data( heartcycle_starts(k) : heartcycle_ends(k) );
	max_vector_of_heartcycles(k) = max( abs( audio_portion ) );
	std_vector_of_heartcycles(k) = std( audio_portion );
end

%% Szint paraméter vektor
function [level_features_of_segment, kurtosis_feature_of_segment] = get_level_features_of_segment(audio_decomposition, fs_decomp, hs_start, hs_stop, n_parts, n_kurt_parts, max_vector_of_heartcycles, std_vector_of_heartcycles, envelope_type, to_normalize)

hs_stop = hs_stop(hs_stop < length(audio_decomposition));
hs_stop = round(hs_stop(hs_stop > hs_start(1)));
hs_start = round(hs_start(hs_start < hs_stop(end)));

if( length(hs_stop) ~= length(hs_start) )
	disp('error')
	return
end

audio_decomposition = audio_decomposition(:)';
homomorphic_envelope = Homomorphic_Envelope_with_Hilbert(audio_decomposition, fs_decomp, 20);
hilbert_transform = abs(hilbert(audio_decomposition));

wavelet_features = zeros(length(hs_start), n_parts);
kurtosis_features = zeros(length(hs_start), sum(n_kurt_parts));

for i = 1 : length(hs_start)
	len_seg = hs_stop(i) - hs_start(i) + 1;
	audio_portion = audio_decomposition( hs_start(i) : hs_stop(i) );
	kurtosis_features(i, 1) = kurtosis(audio_portion, 1);
% 	len_kur_subseg = floor(len_seg / n_kurt_parts(2)) + 1;
	seg_kur_starts = 1 + round([0 : n_kurt_parts(2)-1] .* len_seg ./ n_kurt_parts(2));
	seg_kur_limits = [seg_kur_starts, len_seg];
	for k = 1 : n_kurt_parts(2)
		kurtosis_features(i, 1+k) = kurtosis(audio_portion(seg_kur_limits(k) : seg_kur_limits(k+1)), 1);
	end
	
	switch to_normalize
		case 'Max'
			audio_portion = audio_portion / max_vector_of_heartcycles(i);
		case 'Std'
			audio_portion = audio_portion / std_vector_of_heartcycles(i);
		case 'None'
% 			
	end
	switch envelope_type
		case 'Shannon'
			envelope = audio_portion;
			envelope = power(envelope, 2);
			envelope = envelope .* (-log(envelope));
			envelope(isnan(envelope)) = 0;
			envelope(envelope<0) = 0;
		case 'Hilbert'
			envelope = hilbert_transform( hs_start(i) : hs_stop(i) );
		case 'Homomorphic'
			envelope = homomorphic_envelope( hs_start(i) : hs_stop(i) );
	end
	
	len_subseg = floor(len_seg / n_parts) + 1;
	envel_seg_starts = 1 + round([0 : n_parts-1] .* len_seg ./ n_parts)';
	envel_seg_starts = repmat(envel_seg_starts, 1, len_subseg);
	envel_offsets = repmat([0 : len_subseg-1], n_parts, 1);
	envel_index = envel_offsets + envel_seg_starts;
	i_skip = envel_index > length(envelope);
	envel_index(i_skip) = length(envelope);
	envel_index = reshape(envel_index', 1, len_subseg*n_parts);
	envelope_rs = reshape(envelope(envel_index), [len_subseg, n_parts])';
	
% 	features(i, :) = mean(envelope_rs, 2);	0821
	wavelet_features(i, :) = trimmean(envelope_rs, 20, 2);	%0821
end

level_features_of_segment = mean(wavelet_features);
kurtosis_feature_of_segment = trimmean(kurtosis_features, 20);


%% Szint kurtosis
% function kurtosis_level_features_of_segment = get_kurtosis_level_features_of_segment(audio_decomposition, fs_decomp, hs_start, hs_stop, n_parts, max_vector_of_heartcycles, std_vector_of_heartcycles, envelope_type, to_normalize)



