function classifyResult = challenge(recordName)
%
% Sample entry for the 2016 PhysioNet/CinC Challenge.
%
% INPUTS:
% recordName: string specifying the record name to process
%
% OUTPUTS:
% classifyResult: integer value where
%                     1 = abnormal recording
%                    -1 = normal recording
%                     0 = unsure (too noisy)
%
% To run your entry on the entire training set in a format that is
% compatible with PhysioNet's scoring enviroment, run the script
% generateValidationSet.m
%
% The challenge function requires that you have downloaded the challenge
% data 'training_set' in a subdirectory of the current directory.
%    http://physionet.org/physiobank/database/challenge/2016/
%
% This dataset is used by the generateValidationSet.m script to create
% the annotations on your training set that will be used to verify that
% your entry works properly in the PhysioNet testing environment.
%
%
% Version 1.0
%
%
% Written by: Chengyu Liu, Fubruary 21 2016
%             chengyu.liu@emory.edu
%
% Last modified by: Chengyu Liu, April 3 2016
%             Note: using 'audioread' to replace the 'wavread' function to load the .wav data.
%
% Last modified by: Péter Hajas and Márton Áron Goda

%% 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');

%% Load data and resample data
springer_options   = default_Springer_HSMM_options;
[PCG, Fs1]         = audioread([recordName '.wav']);  % load data
PCG_resampled      = resample(PCG,springer_options.audio_Fs,Fs1); % resample to springer_options.audio_Fs (1000 Hz)

fs = springer_options.audio_Fs;
[startd, stopd] = select_nonsilent_longest(PCG_resampled, fs, 0);
PCG_resampled = PCG_resampled(startd:stopd);

if((stopd-startd) < 1.8*fs)
	classifyResult = 0;
	return;
end

%% Noisy?
load('svm_model_noise7.mat');
audio_data = PCG_resampled;
audio_data = normalise_signal(audio_data);
nf = get_noise_features(audio_data, fs, 0);
[~, score] = predict(ScoreSVMModel, nf);
if(score(2)<0.3)
	classifyResult = 0;
	return;
end

PCG_resampled = PCG_resampled( 1 : min(length(PCG_resampled), 80 * fs) );

audio_data = PCG_resampled;
[b, a] = butter(3, 400/(fs/2), 'low');
audio_data = filtfilt(b, a, audio_data);
[b, a] = butter(3, 20/(fs/2), 'high');
audio_data = filtfilt(b, a, audio_data);
audio_data = normalise_signal(audio_data);

%% 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

%% Running extractFeaturesFromHsIntervals.m to obtain the features for normal/abnormal heart sound classificaiton
% features  = extractFeaturesFromHsIntervals(assigned_states,PCG_resampled);

%% Running classifyFromHsIntervals.m to obtain the final classification result for the current recording
% classifyResult = classifyFromHsIntervals(features);

%% Features
[hs_boundaries] = get_hs_boundaries_VS(assigned_states, fs);
% hs_boundaries = hs_boundaries(2:end, :); 
[m, ~] = size(hs_boundaries);
if(m<4)
	classifyResult = 0;
	return;
end

width_S1 = trimmean( (hs_boundaries(:, 2) - hs_boundaries(:, 1)), 10);
width_S2 = trimmean( (hs_boundaries(:, 4) - hs_boundaries(:, 3)), 10);

len_fft_sys = 512;
len_fft_dia = 1024;
dia_freq11 = round([25, 35] * len_fft_dia / (fs/2));
dia_freq12 = round([80, 120] * len_fft_dia / (fs/2));
sys_freq11 = round([28, 40] * len_fft_sys / (fs/2));
sys_freq12 = round([80, 120] * len_fft_sys / (fs/2));
dia_freq21= round([3, 7] * len_fft_dia / (fs/2));
sys_freq21= round([5, 12] * len_fft_sys / (fs/2));
dia_freq31= round([5, 25] * len_fft_dia / (fs/2));
sys_freq31= round([8, 25] * len_fft_sys / (fs/2));
sys_ratio_freq = round([20, 300] * len_fft_sys / (fs/2));
dia_ratio_freq = round([20, 300] * len_fft_dia / (fs/2));

[~, fft_sys_mat] = get_seg_fft_amp_mat_v01(audio_data, fs, hs_boundaries(:, 2), hs_boundaries(:, 3), len_fft_sys);
fft_sys_norm = fft_amp_mat2fft_mean_v01(fft_sys_mat, 1);
fft_sys_orig = fft_amp_mat2fft_mean_v01(fft_sys_mat, 0);
fft_sys_mat2 = abs(fft(fft_sys_norm,len_fft_sys*2));
fft_sys_mat2 = fft_sys_mat2(1:end/2);
fft_sys_mat3 = abs(fft(fft_sys_mat2,len_fft_sys*2));
fft_sys_mat3 = fft_sys_mat3(1:end/2);

[~, fft_dia_mat] = get_seg_fft_amp_mat_v01(audio_data, fs, hs_boundaries(:, 4), hs_boundaries(:, 5), len_fft_dia);
fft_dia_norm = fft_amp_mat2fft_mean_v01(fft_dia_mat, 1);
fft_dia_orig = fft_amp_mat2fft_mean_v01(fft_dia_mat, 0);
fft_dia_mat2 = abs(fft(fft_dia_norm,len_fft_dia*2));
fft_dia_mat2 = fft_dia_mat2(1:end/2);
fft_dia_mat3 = abs(fft(fft_dia_mat2,len_fft_dia*2));
fft_dia_mat3 = fft_dia_mat3(1:end/2);

dia_fft_sum11 = sum(fft_dia_norm(dia_freq11));
dia_fft_sum12 = sum(fft_dia_norm(dia_freq12));
sys_fft_sum11 = sum(fft_sys_norm(sys_freq11));
sys_fft_sum12 = sum(fft_sys_norm(sys_freq12));

dia_fft_sum21 = sum(fft_dia_mat2(dia_freq21));
sys_fft_sum21 = sum(fft_sys_mat2(sys_freq21));

dia_fft_sum31 = sum(fft_dia_mat3(dia_freq31));
sys_fft_sum31 = sum(fft_sys_mat3(sys_freq31));

sys_dia_ratio = sum(fft_sys_orig(sys_ratio_freq)) / sum(fft_dia_orig(dia_ratio_freq));

n_levels = 3;
wavelet_type = 'db5';
decomp_levels = n_levels;
audio_data = audio_data(1 : floor(length(audio_data) / (2^n_levels)) * (2^n_levels) );
if(length(audio_data) > hs_boundaries(end, 5))
	hs_boundaries = hs_boundaries(1:end-1, :);
end
%envelope_type = 'Hilbert';
%envelope_type = 'Shannon';
envelope_type = 'Homomorphic';
to_normalize = 'Max';
%to_normalize = 'Std';
%to_normalize = 'None';
figures = 0;
n_hs_kurt = [1, 1, 1, 1; 2, 3, 2, 5];%, kurtosis_column_indices
n_hs_portions = [12, 24, 12, 48];
hs_start_corr = ones(1,4)*-0.020;
hs_stop_corr  = ones(1,4)*-0.020;
hs_corr = [hs_start_corr; hs_stop_corr];
%[wavelet_morphology_features, column_indices] = get_wavelet_morphology_vector_v01(audio_data, fs, wavelet_type, decomp_levels, hs_boundaries, envelope_type, to_normalize, figures);
% [wavelet_morphology_features, wave_col_inds, kurtosis_features, kurt_col_inds] = get_wavelet_morphology_vector_v04(audio_data, fs, wavelet_type, decomp_levels, hs_boundaries, envelope_type, to_normalize, n_hs_portions, hs_corr, n_hs_kurt, figures);
[wavelet_morphology_features, wave_col_inds, kurtosis_features, kurt_col_inds] = get_wavelet_morphology_vector_v04(audio_data, fs, wavelet_type, decomp_levels, hs_boundaries, envelope_type, to_normalize, n_hs_portions, hs_corr, n_hs_kurt, figures);
%load('dim_red_aux_data.mat');
load('dim_red_wav_data.mat');	%save('dim_red_wav_data.mat', 'daub_column_indices_new', 'daub_base_vectors', 'daub_orig_feature_means');
load('dim_red_kur_data.mat');	%save('dim_red_kur_data.mat', 'kurtosis_column_indices_new', 'kurtosis_base_vectors', 'kurtosis_orig_feature_means');

%% Morphology vector transformation
wf = wavelet_morphology_features(1, :);
wfo = zeros(1, daub_column_indices_new(end)-1);
n_levels = length(wave_col_inds) - 2;
for i = 1 : n_levels+1
	feature_part = wf(1, wave_col_inds(i) : wave_col_inds(i+1)-1);
	part_mean = daub_orig_feature_means(1, wave_col_inds(i) : wave_col_inds(i+1)-1);
	feature_part = feature_part - part_mean;
	base_vectors_of_level = daub_base_vectors{1, i};
	wfo(1, daub_column_indices_new(i) : daub_column_indices_new(i+1)-1) = feature_part * base_vectors_of_level;
end
daub_shape_features_reduced = wfo;

%% Kurtosis transformation
kf = kurtosis_features(1, :);
kfo = zeros(1, kurtosis_column_indices_new(end)-1);
n_levels = length(kurt_col_inds) - 2;
for i = 1 : n_levels+1
	feature_part = kf(1, kurt_col_inds(i) : kurt_col_inds(i+1)-1);
	part_mean = kurtosis_orig_feature_means(1, kurt_col_inds(i) : kurt_col_inds(i+1)-1);
	feature_part = feature_part - part_mean;
	base_vectors_of_level = kurtosis_base_vectors{1, i};
	kfo(1, kurtosis_column_indices_new(i) : kurtosis_column_indices_new(i+1)-1) = feature_part * base_vectors_of_level;
end
daub_kurtosis_features_reduced = kfo;



pcg_old_props = [width_S1; width_S2];
% sys_dia_fft_energy_props = [sys_fft_sum1; sys_fft_sum2; dia_fft_sum1; dia_fft_sum2; sys_dia_ratio];
sys_dia_fft_energy_props = [sys_fft_sum11; sys_fft_sum12; dia_fft_sum11; dia_fft_sum12; sys_dia_ratio];
sys_dia_new = [sys_fft_sum21; dia_fft_sum21; sys_fft_sum31; dia_fft_sum31];
% feature_vector = [pcg_old_props; sys_dia_fft_energy_props; daub_shape_features_reduced'];
% feature_vector = [pcg_old_props; sys_dia_fft_energy_props; sys_dia_new; daub_shape_features_reduced'];
feature_vector = [pcg_old_props; sys_dia_fft_energy_props; sys_dia_new; daub_shape_features_reduced'; daub_kurtosis_features_reduced'];

%% Calssification
load('svm_model.mat');
[classifyResult,~] = predict(SVMModel, feature_vector');

