function[output, idx_local, beginning, ending] = Extraction(file)
% Funkce slouzi k extrakci sinové aktvity
%
% Vstup: file - jeden objekt obsahujici data
%
% Vystup: output - binárni vektor, kde 1 značí výskyt sinove aktivity a 0
%                  nikoliv
%         idx_local - misto zaznamenanych piku sinove aktivity ve vzorcich
%         beginning - misto, urcene jako zacatek prubehu sinove aktivity
%         ending - misto, urcene jako konec prubehu sinove aktivity
%
%% load
    %file = jsondecode(fileread('C:\Users\zanet\Documents\BP\data_new\000043.json'))
    CS5 = file.data(13).content';
    CS4 = file.data(14).content';
    CS3 = file.data(15).content';
    CS2 = file.data(16).content';
    CS1 = file.data(17).content';
    all_CS = [CS5; CS4; CS3; CS2; CS1]';
    orig_CS = all_CS;
    
fvz = file.measurement.samplingFrequency; %sampling frequency
pom = 0;
t = [];
for k = 1 : length(all_CS) 
   pom = pom + (1/fvz);
   t = [t, pom];   %timeline
end
    
%% subtraction of the effect of the DC component (odecteni vlivu stejnosmerne slozky)

all_CS = all_CS - mean(all_CS); 

%% linear filtration HighPass 50 Hz

[n, Wn] = buttord(55/(fvz/2), 50/(fvz/2), 3, 5); %n...minimul order, Wn...cutoff freq
[b, a] = butter(n, Wn, 'high'); %designing filter
all_CS = filter(b, a, all_CS); %filtering

%% envelope of the signal - nonlinear transformation

%all_CS = abs(hilbert(all_CS));
[all_CS, ~] = envelope(all_CS, fvz);
all_CS = all_CS.^2;

%% smoothing - linear filtration LP

L = 230; 
f = 1/0.05; %50 ms (1/period)
Wn = f/(fvz/2);
bLP = fir1(L, Wn,'low');
a = 1;

all_CS = filtfilt(bLP, a, all_CS);


 %% PCA 
%[coeff, scorec ,~ ,~ ,variance ,envelope_mu] = pca(all_CS'); %  mu = mean
[coeff] = pca(all_CS');

% selecting top 1 best feature: 
envelope_CS = coeff(:,1)'; %reduced dimension - 1D

%% peak detection

L = length(envelope_CS); 
window = 50; 

% looking for first treshold
learning_phase = fvz * 2;
if learning_phase > L
    learning_phase = L;
end
threshold = max(envelope_CS) * 0.5; 
[value_local, idx_local] = findpeaks(envelope_CS(1 : learning_phase), 'MinPeakHeight', threshold);
if isempty(value_local)
    threshold = max(envelope_CS(1:learning_phase)) * 0.5; 
    [value_local, idx_local] = findpeaks(envelope_CS(1 : learning_phase), 'MinPeakHeight', threshold);
end
%number_peaks = size(idx_local, 2);

%first treshold
    %moving average from detected peaks
    threshold = movmean(value_local,8);
    threshold = threshold(1, size(threshold, 2))*0.55;
    %mean from detected peaks
    %threshold = mean(value_local)*0.55; % new treshold from detected peaks

i = 1 + learning_phase;
while i < L-window
    
    if envelope_CS(i) > threshold 
        [value, idx] = max(envelope_CS(i : i+window)); % maximum in window if amplitude above treshold
        idx_local = [idx_local, i+idx]; % index of maximum
        value_local = [value_local, value]; % value of maximum
        i = idx_local(end) + (fvz * 0.15); % refractory phase of the myocardium (0.15 sec)
        %number_peaks = number_peaks + 1;
    else
        i = i + 1; % looking for another value above threshold
    end
    %moving average from detected peaks
        threshold = movmean(value_local, 10); 
        threshold = threshold(1, size(threshold, 2))*0.4; % new threshold from detected peaks
    %mean from detected peaks
        %threshold = mean(value_local)*0.55
end

%% AUC - area uder curve

window = 100; % pair number

beginning = [];
ending = [];

for i = 1:length(idx_local)

    if idx_local(i) <= window+1
        peak = envelope_CS(1 : idx_local(i)+window);
    elseif idx_local(i) >= length(envelope_CS)-window-1
        peak = envelope_CS(idx_local(i)-window : size(envelope_CS, 2));
    else
        peak = envelope_CS(idx_local(i)-window : idx_local(i)+window); % window+1+window ... length_peak = 2*window+1
    end
    x = 1 : length(peak);
    area_o = trapz(x, abs(peak));
    peak_norm = peak./area_o; % normalized peak
    % norm_area = trapz(x, y_normalized) % normalized area, total = 1
   
    perc = 0.95;
    j = 1;
    found = 0;
    
    while found == 0
        if idx_local(i) + window + 1 + j > size(envelope_CS, 2)
            beginning = [beginning, idx_local(i) - j + 1];
            ending = [ending, size(envelope_CS, 2)];
            found = 1;
            
        elseif idx_local(i) - window - j < 1
            beginning = [beginning, 1];
            ending = [ending, idx_local(i) + j - 1];
            found = 1;
            
        else
            section = peak_norm( window+1-j : window+1+j);
            x = 1 : 2*j+1;
            area_part = trapz(x, abs(section));
            if area_part >= perc
                beginning = [beginning, idx_local(i) - j];
                ending = [ending, idx_local(i) + j];
                found = 1;
            else
                j = j + 1;
            end
        end
    end
end

%% output

output = zeros(1, size(envelope_CS, 2));
for i = 1:length(beginning)
    output(beginning(i) : ending(i)) = 1;
end
 %% figure
% figure 
% subplot 211
% plot(t, envelope_CS)
% hold on
% plot(t, output*max(envelope_CS))
% title ('Obálka s vyznačenou síňovou aktivitou')
% legend('Obálka', 'Síňová aktivita')
% xlabel('Čas(s)')
% ylabel('Amplituda(mV)')
% subplot 212
% plot(t, orig_CS(:, 3)')
% hold on
% plot(t, output*max(orig_CS(:, 3)'))
% title ('Obálka s vyznačenou síňovou aktivitou')
% legend('CS3', 'Síňová aktivita')
% xlabel('Čas(s)')
% ylabel('Amplituda(mV)')

end