function[idx_local] = DetectionQRS(file)
% Funkce slouzi k detekci komorove vlny - QRS komplexu.
%
% Vstup: file - jeden objekt obsahujici data
%
% Vystup: idx_local - misto zaznamenanych piku komorove aktivity 
%                     ve vzorcich
%
%% load
%file = jsondecode(fileread('C:\Users\zanet\Documents\BP\data_new\00000D.json'));

    I =  file.data(1).content';

orig_I = I;    
fvz = file.measurement.samplingFrequency; %sampling frequency
pom = 0;
t = [];
for k = 1 : length(I) 
   pom = pom + (1/fvz);
   t = [t, pom];   %timeline
end
    
%% subtraction of the effect of the DC component (odecteni vlivu stejnosmerne slozky)

I = I - mean(I); 

%% linear filtration HighPass 50 Hz
Wp = [11 21]/(fvz/2); Ws = [10 24]/(fvz/2);   % Normalizing frequency
%Rp = 3; Rs = 40;
Rp = 5; Rs = 10;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);  % gives minimum order of filter, n...minimul order, Wn...cutoff freq
[b, a] = butter(n, Wn); % designing Butterworth filter - bandpass
I = filter(b, a, I); %filtering

% plotting freq characteristic
% [z, p, k] = butter(n, Wn); %bandpass
% sos = zp2sos(z, p, k);
% fvtool(sos, 'Analysis', 'freq')

%% envelope of the signal - nonlinear transformation

[I, ~] = envelope(I, fvz);
I = I.^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;

I = filtfilt(bLP, a, I);


envelope_CS = I; %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

%if another peak wasnt found in time of 1.66*duration of peak before
%looking for another peak with different threshold value and lower fraction parametr
for i = 2 : (value_local - 1)
    active = idx_local(i+1) - idx_local(i);
    left = idx_local(i) - idx_local(i-1);
    
    if  active > (left * 1.66)
        j = idx_local(i+1);
        threshold = movmean(value_local, size(value_local, 2));
        threshold = threshold(1, size((i : size(value_local, 2)), threshold, 2))*0.55;
        while j > idx(i) + window + 1
            if envelope_CS(i) > threshold
                [value, idx] = max(envelope_CS(idx_local(i+1)-1 : -1 : idx_local(i+1)-1-window)); 
                idx_local = [idx_local, i + idx];
                value_local = [value_local, value];
            else
                j = j - 1; %step backwards
            end
        end
    end
end

%% figure
% figure
% subplot 211
% plot(t, envelope_CS)
% hold on
% plot(t(idx_local),value_local, 'r+');
% title ('Obálka s vyznačenými detekovanými píky')
% xlabel('Čas(s)')
% ylabel('Amplituda(mV)')
% subplot 212
% plot(t, orig_I)
% hold on
% plot(t(idx_local),orig_I(idx_local), 'r+');
% title ('Signál s vyznačenými detekovanými píky')
% xlabel('Čas(s)')
% ylabel('Amplituda(mV)')
end
