function stanovenaKvalita = networkAssessmentV4(inputData,FsInput)
x = zeros(length(inputData),1);
x(:,1) = inputData; 


% Úprava vstupních dat na vzorkovací frekvenci 250 Hz
% decFactor = FsInput/250;
% x = downsample(x,decFactor);
Fs = 250;
L = lcm(FsInput,Fs);
intFactor = L/FsInput;
decFactor = L/Fs;

x = interp(x,intFactor);
x = downsample(x,decFactor);

window = 3000;
minLen = 128*decFactor/intFactor; % Min. délka vstupních dat

% Volání klasifikační funkce
%tic
if length(inputData) >= minLen
   [stanovenaKvalita,netPreds1] = networkClassification(x,Fs,window,minLen);
else
   txt = sprintf('Příliš krátká vstupní data. Minimální délka je %d vzorků.', minLen);
   disp(txt)
   return
end
%toc

% Volání funkce pro detekci R vln
%tic
[rIndexes, QRSstart, QRSend, netPreds1] = RwaveDetection(x,netPreds1,Fs);
%toc

% Volání funkce pro hodnocení kvality EKG signálu
%tic
% Definice známek kvalit
if length(rIndexes) > 1
    stanovenaKvalita = assessment(stanovenaKvalita,rIndexes,netPreds1,QRSstart,QRSend,Fs,x);
elseif length(rIndexes) == 1
    disp('Našla se pouze jedna R vlna, signál bude označen známkou 3.')
    disp('Jde o nekvalitní nebo příliš krátká vstupní data (min. 2 periody).')
    stanovenaKvalita(1:end) = 3;
else
    disp('Nenašly se žádné R vlny, signál bude označen známkou 3.')
    stanovenaKvalita(1:end) = 3;
end    
%toc

% Výpočet procentuálního zastoupení jednotlivých známek
oneQper = (length(stanovenaKvalita(stanovenaKvalita == 1))/length(stanovenaKvalita))*100;
twoQper = (length(stanovenaKvalita(stanovenaKvalita == 2))/length(stanovenaKvalita))*100;
threeQper = (length(stanovenaKvalita(stanovenaKvalita == 3))/length(stanovenaKvalita))*100;
qPercentage = [oneQper,twoQper,threeQper];

% Interpolace, aby vektor kvalit byl o stejné délce jako vstupní signál
stanovenaKvalita = double(stanovenaKvalita);
stanovenaKvalita = interp(stanovenaKvalita, decFactor);
stanovenaKvalita = round(stanovenaKvalita);
stanovenaKvalita = downsample(stanovenaKvalita, intFactor);
stanovenaKvalita = stanovenaKvalita(1:length(inputData));
end
% -----------------------------------------------------------------------------------------------------------

function [stanovenaKvalita, netPreds1] = networkClassification(x,Fs,window,minLen)
% Segmentace neuronovou sítí a hodnocení na základě SNR
% Vstupy:
%       x - vstupní signál
%       Fs - vzorkovací frekvence (síťová)
%       window - okno s počtem vzorků vstupujících do neuronové sítě
%       minLen - minimální délka signálu
% Výstupy:
%       netPreds1 - kategorický vektor o délce vstupního signálu
%       stanovenaKvalita - vektor kvalit (ne kompletní)

% Rozsekání vstupního signálu na části o délce okna
n = numel(x);
inputCell = mat2cell(x,diff([0:window:n-1,n]));

% Prealokace vektorů
netPreds1 = [];
netPreds1 = categorical(netPreds1);
netPreds1(1:length(x)) = 'P';
stanovenaKvalita = zeros(1,length(x));
SNRlimit = 26.5;
inputCellLen = length(inputCell);
load fsstNet.mat fsstNet
% Klasifikace po jednotlivých částech vstupního signálu a výpočet SNR
for i = 1:inputCellLen-1
    input = inputCell{i};
    [mitFSST, SNR] = computeFSST(input,Fs);
    netPreds = classify(fsstNet,mitFSST,'MiniBatchSize',40);
    netPreds1(i*window-window+1:i*window) = netPreds{1};
    if SNR < SNRlimit
       stanovenaKvalita(i*window-window+1:i*window) = 2;
    end
end

% Klasifikace zbytkových dat na konci signálu a výpočet SNR
input = inputCell{end};
inputLen = length(input);
% Pro klasifikaci/SNR potřebujeme signál delší než min. délka
if inputLen >= minLen
   [mitFSST, SNR] = computeFSST(input,Fs);
   netPreds = classify(fsstNet,mitFSST,'MiniBatchSize',40);
   netPreds1(inputCellLen*window-window+1:inputCellLen*window-window+inputLen) = netPreds{1};
   if SNR < SNRlimit
      stanovenaKvalita(inputCellLen*window-window+1:inputCellLen*window-window+inputLen) = 2;
   end
end    
end   


function [rIndexes, QRSstart, QRSend, netPreds1] = RwaveDetection(x,netPreds1,Fs)
% Detekce R vln, začátků a knců QRS komplexů
% Vstupy:
%       x - vstupní signál
%       netPreds1 - kategorický vektor o délce vstupního signálu
%       Fs - vzorkovací frekvence (síťová)
% Výstupy:
%       rIndexes - vektor pozic R vln
%       QRSstart - vektor pozic začátků QRS komplexů
%       QRSend - vektor pozic konců QRS komplexů
%       netPreds1 - kategorický vektor o délce vstupního signálu

minQRSlen = round(0.06*Fs); % minimální šířka QRS komplexu

% Prealokace vektoru
prealocatedVector = zeros(1,round(length(x)/(0.75*Fs)));
QRSstart = prealocatedVector;
QRSend = prealocatedVector;
rIndexes = prealocatedVector;

% Výběr indexů vzorků z kat. vektoru, které odpovídají QRS komplexům
QRSindexes = find(netPreds1 == 'QRS');
% Rodělení jednotlivých QRS komplexů (indexy) do samostatných buněk
QRSgrouped = mat2cell( QRSindexes, 1, diff( [0, find(diff(QRSindexes) ~= 1), length(QRSindexes)] )) ;
% Výběr QRS komplexů, které splňují minimální délku
QRSgroupedTrue = QRSgrouped(cellfun('length', QRSgrouped) >= minQRSlen);
% Osamocené (falešné) QRS vzorky přepiš na 'n/a'
netPreds1(cell2mat(QRSgrouped(cellfun('length', QRSgrouped) <= 1))) = 'n/a';

% Vytvoření vektorů začátků a konců QRS komplexů a detekce R vln
for i = 1:length(QRSgroupedTrue)
    firstQRSSampleIndex = QRSgroupedTrue{i}(1);
    lastQRSSampleIndex = QRSgroupedTrue{i}(end);
    QRSstart(i) = firstQRSSampleIndex;
    QRSend(i) = lastQRSSampleIndex;
    [maxValue, maxIndex] = max(x(QRSgroupedTrue{i}(1:end)));
    [minValue, minIndex] = min(x(QRSgroupedTrue{i}(1:end)));
    % Rozhodnutí zda R vlna směřuje nahoru nebo dolů
    if abs(maxValue - x(firstQRSSampleIndex)) > abs(minValue - x(firstQRSSampleIndex))
       rIndexes(i) = QRSgroupedTrue{i}(maxIndex);
    else
       rIndexes(i) = QRSgroupedTrue{i}(minIndex); 
    end   
end

% Vymazaní zbytkových 0 pocházejících z prealokace
QRSstart = nonzeros(QRSstart);
QRSend = nonzeros(QRSend);
rIndexes = nonzeros(rIndexes);
end

function stanovenaKvalita = assessment(stanovenaKvalita,rIndexes,netPreds1,QRSstart,QRSend,Fs,x)
% Hodnocení kvality signálu
% Vstupy:
%       stanovenaKvalita - vektor kvalit (ne kompletní)       
%       rIndexes - vektor pozic R vln
%       netPreds1 - kategorický vektor o délce vstupního signálu
%       QRSstart - vektor pozic začátků QRS komplexů
%       QRSend - vektor pozic konců QRS komplexů
%       Fs - vzorkovací frekvence (síťová)
%       x - vstupní signál
% Výstupy:
%       stanovenaKvalita - stanovený vektor kvality EKG (finální)

% Vytvoření vektoru period
periods = rIndexes(2:end)-rIndexes(1:end-1);

% Výpočty limitních počtů vzorků pro jednitlivé úseky a stanovení rozměrů úseků
maxPeriod = round(2*Fs); % Nejdelší možná perioda
minPeriod = round(6*Fs/25); % Nejkratší možná perioda
minP = round(0.06*Fs); % Minimální počet vzorků P vlny
minT = round(0.1*Fs); % Minimální počet vzorků T vlny
PRinterval = round(0.3*Fs);
STinterval = round(0.4*Fs);
STsegment = round(0.15*Fs);
PRsegment = round(0.12*Fs);

% Částem nesplňujícím podmínku maximální HRV přiřaď známku 3
counter = 0;
HRVwindow = 20;
STDlimit = 30;
periodsSelected = zeros(1,HRVwindow);
j = 1;
for i = 2:length(rIndexes)
    % HRV se počítá z period o počtu HRVwindow
    periodsSelected(j) = rIndexes(i)-rIndexes(i-1);
    counter = counter + 1;
    j = j + 1;
    if counter == HRVwindow
       periodsSTD = std(periodsSelected)*100/mean(periodsSelected);
       counter = 0;
       if periodsSTD > STDlimit
          stanovenaKvalita(rIndexes(i-HRVwindow):rIndexes(i)) = 3;
       end
       periodsSelected = zeros(1,HRVwindow);
       j = 1;
    end
end

% Výpočet HRV pro konec signálu, kde je počet period již menší než HRVwindow
rIndexesLen = length(rIndexes);
if rIndexesLen > HRVwindow
   missingPartlen = rem(rIndexesLen,HRVwindow);
   if missingPartlen > 1
      periodsSelected = rIndexes(rIndexesLen-missingPartlen+1:rIndexesLen) - rIndexes(rIndexesLen-missingPartlen:rIndexesLen-1);
      periodsSTD = std(periodsSelected)*100/mean(periodsSelected);
      if periodsSTD > STDlimit
         stanovenaKvalita(rIndexes(end-missingPartlen):rIndexes(end)) = 3;
      end
   end
else
   % Pro výpočet HRV potřebujeme aspoň 2 periody
   if rIndexesLen > 1
      periodsSelected = rIndexes(2:rIndexesLen) - rIndexes(1:rIndexesLen-1);
      periodsSTD = std(periodsSelected)*100/mean(periodsSelected);
      if periodsSTD > STDlimit
         stanovenaKvalita(rIndexes(1):rIndexes(end)) = 3;
      end
   end   
end   

% Doznámkování zbytku signálu po jednotlivých periodách
for i = 1:length(QRSend)-1
    % Pokud již perioda má známku 2 nebo 3 pokračuj na další
    if all(stanovenaKvalita(rIndexes(i)+1:rIndexes(i+1)) == 3)
       continue
    end
    % Najdi falešné QRS vzorky mimo skutečný QRS komplex
    falseQRS = find(netPreds1(QRSend(i)+STsegment:QRSstart(i+1)-PRsegment)=='QRS',1);
    % Pokud je perioda delší nebo kratší než limit nebo byly nalezeny mezi
    % QRS komplexy falešné QRS vzorky, označ periodu známkou 3 a pokračuj na další
    if periods(i) > maxPeriod || periods(i) < minPeriod || ~isempty(falseQRS)
       stanovenaKvalita(rIndexes(i):rIndexes(i+1)) = 3;
       continue
    end

    if all(stanovenaKvalita(rIndexes(i)+1:rIndexes(i+1)) == 2)
       continue
    end
   
    % Spočítej počet P vzorků v PR intervalu a T vzorků v ST intervalu 
    Tsamples = find(netPreds1(QRSend(i):QRSend(i)+STinterval)=='T');
    if PRinterval > QRSstart(i+1)
       Psamples = find(netPreds1(QRSstart(i+1)-QRSstart(i+1)+1:QRSstart(i+1))=='P');
    else   
       Psamples = find(netPreds1(QRSstart(i+1)-PRinterval:QRSstart(i+1))=='P');
    end
    
    % Podle splnění minimálního počtu vzorků v ST a PR intervalu uděl známku periodě
    if length(Tsamples) > minT && length(Psamples) > minP
       stanovenaKvalita(rIndexes(i):rIndexes(i+1)) = 1;
    else
       stanovenaKvalita(rIndexes(i):rIndexes(i+1)) = 2;
    end
end

% Necelé periody před/za první/poslední R vlnou označ jako 0, aby
% neovlivnily následující výpočet
stanovenaKvalita(1:rIndexes(1)) = 0;
stanovenaKvalita(rIndexes(end):length(x)) = 0;

% Rozsekání vektoru stanovenaKvalita na části o délce Qwindow
Qwindow = 3000;
m = numel(stanovenaKvalita);
stanovenaKvalita = transpose(stanovenaKvalita);
QCell = mat2cell(stanovenaKvalita,diff([0:Qwindow:m-1,m]));

% Výpočet procentuálního zastoupení známek v daném okně 
for i = 1:length(QCell)
    Qinput = QCell{i};
    QinputLen = length(Qinput);
    oneQper = (length(Qinput(Qinput == 1))/QinputLen)*100;
    twoQper = (length(Qinput(Qinput == 2))/QinputLen)*100;
    threeQper = (length(Qinput(Qinput == 3))/QinputLen)*100;
    qPercentage = [oneQper,twoQper,threeQper];
    % Pokud je v okně převažující známka 3, označ celé ono jako 3
    if max(qPercentage) == threeQper
       stanovenaKvalita(i*Qwindow-Qwindow+1:i*Qwindow-Qwindow+QinputLen) = 3;
    end
end    

% Necelé periody před/za první/poslední R vlnou označ jako 3
stanovenaKvalita(1:rIndexes(1)) = 3;
stanovenaKvalita(rIndexes(end):length(x)) = 3;
end