% Written by Péter Hajas and Márton Áron Goda
% Using some of Springer's functions
function [nf] = get_noise_features(audio_data, fs, figures)


%% 25-400Hz 4th order Butterworth band pass
audio_data = butterworth_low_pass_filter(audio_data,2,400,fs, false);
audio_data = butterworth_high_pass_filter(audio_data,2,25,fs);

%% Spike removal from the original paper:
audio_data = schmidt_spike_removal(audio_data,fs);

%% Find the homomorphic envelope
homomorphic_envelope = Homomorphic_Envelope_with_Hilbert(audio_data, fs);

%% Find the autocorrelation:
y=homomorphic_envelope-mean(homomorphic_envelope);
[c] = xcorr(y,'coeff');
signal_autocorrelation = c(length(homomorphic_envelope)+1:end);
% autocorr_null_level = -1 * min(signal_autocorrelation);
% signal_autocorrelation = signal_autocorrelation + autocorr_null_level;
% signal_autocorrelation = signal_autocorrelation / signal_autocorrelation(1);

min_index = 0.4*fs;
max_index = 1.8*fs;

[pks, i_pks, w, p] = findpeaks( signal_autocorrelation(min_index : max_index), 'MinPeakDistance', min_index);
% locs = i_pks + min_index - 1;

% [~, i_hr] = min(abs(ref_period-locs));
[~, i_max] =  max(pks);

peak = pks(i_max);
w = w(i_max);
p = p(i_max);

%% Kurtosis
win = 2 * fs;
overlay = 0.5 * fs;

n_win = floor((length(audio_data) - win) / overlay);
x_starts = zeros(1, length(n_win));
PCG_kur  = zeros(1, length(n_win));
PCG_rms  = zeros(1, length(n_win));
PCG_fftm = zeros(5, length(n_win));
for m = 1 : n_win
	x_starts(m) = (m-1)*overlay+1;
	sig = audio_data( x_starts(m) : x_starts(m) + win - 1 );
	sig = sig - mean(sig);
	PCG_kur(m) = kurtosis(sig, 0, 1);
	PCG_rms(m) = std(envelope(sig, 100, 'rms'));
	PCG_fft = abs(fft(sig, 2048));
	PCG_fft = PCG_fft(1:end/2);
	PCG_fft = PCG_fft / sum(PCG_fft);
	PCG_fftm(1, m) = mean(PCG_fft(1:32));
	PCG_fftm(2, m) = mean(PCG_fft(33:64));
	PCG_fftm(3, m) = mean(PCG_fft(65:128));
	PCG_fftm(4, m) = mean(PCG_fft(129:384));
	PCG_fftm(5, m) = mean(PCG_fft(385:1024));
end

s_kur = sort(PCG_kur, 'ascend');
s_rms = sort(PCG_rms, 'ascend');
PCG_fftm_out = trimmean(PCG_fftm, 20, 2)';
nf = [peak, w, trimmean(PCG_kur, 20), trimmean(PCG_rms, 20),...%rms: b halmazt jól kiválasztja
		PCG_fftm_out,...
		s_kur(round([end*0.25, end*0.50, end*0.75])), s_rms(round([end*0.25, end*0.50, end*0.75]))];














