% Skript predstavuje nosnu cast pre spracovanie umelych signalov a nasledne podklad pre vlnkovy a wienerovsky filter, obsahuje:
% 1. generovanie umelych signalov s roznymi parametrami SNR
% 2. filtraciu s vyuzitim vlnkovej transformacie a prahovania s roznymi parametrami vlnky, prahu, prahovania a  poctu filtrovnych pasiem
% 3. odhad SNR vypoctom z vyfiltrovaneho signalu po prahovani (uzitocna zlozka) a odfiltrovaneho rusenia zo signalu
% Autor: Martina Babicova

% Prametry, ktore sa menia:
% 1. typy vlniek su 3: ''bior1.1'', ''bior1.3'', ''bior1.5''
parametreVlnky = ['bior1.1'; 'bior1.3'; 'bior1.5'];
% 2. hodnota prahu moze byt v dvoch variantach Univerzalny a Empiricky prah
parametrePrah = ['U', 'E'];
% 3. hodnota SNR (predstavuje 17 rozne zarusenych signalov)
parametreSNR = [-20, -15, -10, -7, -5, -3, -1, 0, 1, 3, 5, 7, 10, 15, 20, 25, 30];
% 4. prahovanie moze byt v troch variantach: tvrde, makke, hybridne
% 5. pocetFiltrovanychPasiem (meni sa od 1 do 3 pre umele signaly)

% Hlavna cast - vytvorenie .txt suborov s obsahom vsetkych vypoctov:
% vytvor 3 textove subory (pre kazdu vlnku z moznosti), do ktoreho sa budu postupne vkladat jednotlive kombinacie parametrov   
for vlnkaParameter = 1:size(parametreVlnky, 1)           % vyber vlnku
    obsahSuboru = "";
    vlnka = parametreVlnky(vlnkaParameter,:);            % dvojbodka preto, aby to vzalo cely 'bior1.1' string a nie len 'b'
    fid = fopen(vlnka + ".txt", 'wt');
for prahParameter = 1:length(parametrePrah)              % vyber typ prahu (U - univerzalny, E - empiricky)
    typPrahu = parametrePrah(prahParameter);
    obsahSuboru = obsahSuboru + "THRESHOLD " + typPrahu + "\n";
    obsahSuboru = obsahSuboru + "Vypocitane SNR\tT prahovanie 3\tT prahovanie 2\tT prahovanie 1\tM prahovanie 3\tM prahovanie 2\tM prahovanie 1\tH prahovanie 3\tH prahovanie 2\tH prahovanie 1";
for SNRParameter = 1:length(parametreSNR)                  % vyber skusanu hodnotu SNR zo 17 moznosti
    SNRp = parametreSNR(SNRParameter);                     % [-20, -15, -10, -7, -5, -3, -1,0, 1, 3, 5, 7, 10, 15, 20, 25, 30]
    obsahSuboru = obsahSuboru + "\n" + SNRp + "\t";
for prahovanie = 0:2                                       % vyber typ prahovania (0 - tvrde, 1 - makke, 2 - hybridne)
for pocetFiltrovanychPasiem = 1:3                          % vyber pocet filtrovanych pasiem (1 až 3)
%% Generovanie signalu EKG a rusenia v pomere podla zadaneho SNR
[signalEKGSRusenim] = UpravaSignaluPodlaSNR(40, 250, 10000, SNRp); % premenna signalEKGSRusenim obsahuje signal EKG s EMG podla zadaneho SNR
% na upravu generovaneho signalu je pouzita funkcia UpravaSignaluPodlaSNR
%% Vlnkova transformacia - SWT
stupenDekompozicie = 4;                                                      % stupen dekompozicie (parameter pre swt)
[signalPoRozkladeWT] = swt(signalEKGSRusenim, stupenDekompozicie, vlnka);    % funkcia swt rozdeli signal do 5 frekvencnych pasiem
vykreslenie = 0;                                                             % pre vykreslenie je potrebne nastavit premennu vykreslenie = 1;
if vykreslenie == 1  % Vykreslenie pasiem po SWT
    figure()
    subplot 511
    plot(signalPoRozkladeWT(1,:))
    title('Frekvencne pasmo 62.5 - 125 Hz')
    subplot 512
    plot(signalPoRozkladeWT(2,:))
    title('Frekvencne pasmo 31.25 - 62.5 Hz')
    subplot 513
    plot(signalPoRozkladeWT(3,:))
    title('Frekvencne pasmo 15.625 - 31.125 Hz')
    subplot 514
    plot(signalPoRozkladeWT(4,:))
    title('Frekvencne pasmo 7.78125 - 15.5625 Hz')
    subplot 515
    plot(signalPoRozkladeWT(5,:))
    title('Frekvencne pasmo 0 - 7.78125 Hz')
end

%% Prahovanie
dlzkaOkna = 250;
% pri dlzke signalu = 10 000 okne = 250 je empirickyPrah, univerzalnyPrah, smerodajnaOdchylka aj SNR s poctom slpcov 40 (10 000: 250 = 40)
% pocet riadkov = maximalny pocet filrovanych pasiem
% vytvor premenne pre empiricky, univerzalny prah, smerodajnu odchlylku a signal po prahovani vyplnene nulami
empirickyPrah = zeros(stupenDekompozicie -1,((length(signalPoRozkladeWT)))/dlzkaOkna);
univerzalnyPrah = zeros(stupenDekompozicie -1,((length(signalPoRozkladeWT)))/dlzkaOkna);
smerodajnaOdchylka = zeros(stupenDekompozicie -1,((length(signalPoRozkladeWT)))/dlzkaOkna);
signalPoPrahovani = zeros(stupenDekompozicie + 1,length(signalPoRozkladeWT)); 
% do premennych postupne v cykle pridavaj aktualne hodnoty pre konkretne plavajuce okno
% v cykle prechadzaj s krokom dlzkaOkna cely signal a v ramci jedneho okna
% spocitaj smerodajnu odychylku, empricky prah, univerzalny prah
for j = 1:dlzkaOkna:(length(signalPoRozkladeWT(pocetFiltrovanychPasiem,:)) - dlzkaOkna) % j = 1, 251, 501, 751 ,1001
    % zrataj smerodajnu odchylku
    smerodajnaOdchylka(pocetFiltrovanychPasiem,(((j-1)/dlzkaOkna)+1)) = (median(abs(signalPoRozkladeWT(pocetFiltrovanychPasiem,j:j+(dlzkaOkna-1)))-median(signalPoRozkladeWT(pocetFiltrovanychPasiem,:))))/0.6745;
    % zrataj empricky prah
    empirickyPrah(pocetFiltrovanychPasiem,(((j-1)/dlzkaOkna)+1)) = smerodajnaOdchylka(pocetFiltrovanychPasiem,(((j-1)/dlzkaOkna)+1)) .* 3; 
    % zrataj univerzalny prah
    univerzalnyPrah(pocetFiltrovanychPasiem,(((j-1)/dlzkaOkna)+1)) = smerodajnaOdchylka(pocetFiltrovanychPasiem,(((j-1)/dlzkaOkna)+1))*sqrt(2*log(length(signalPoRozkladeWT)));
    % podla toho comu je typPrahu rovny pouzi konkretny typ prahu
    if (typPrahu == 'U')
                prah = univerzalnyPrah(pocetFiltrovanychPasiem, (((j-1)/dlzkaOkna)+1));
    else
                prah = empirickyPrah(pocetFiltrovanychPasiem, (((j-1)/dlzkaOkna)+1));
    end
        
    for k = j:(j + dlzkaOkna - 1) % iteracia v ramci okna 
            % porovnaj absolutnu hodnotu cisla v povodnom
            % vektore(signalPoRozkladeWT) na pozicii k s hodnotou prahu
            % a ak je vascia pouzi prahovanie a uloz do noveho vektora na poziciu k (signalPoPrahovani)   
        if abs(signalPoRozkladeWT(pocetFiltrovanychPasiem,k)) > prah
            if (prahovanie == 0) % tvrde
                signalPoPrahovani(pocetFiltrovanychPasiem,k) = signalPoRozkladeWT(pocetFiltrovanychPasiem,k);                                                    
            elseif (prahovanie == 1) % makke
                signalPoPrahovani(pocetFiltrovanychPasiem,k) = sign(signalPoRozkladeWT(pocetFiltrovanychPasiem,k))*(abs(signalPoRozkladeWT(pocetFiltrovanychPasiem,k)- prah));      
            else % hybridne
                signalPoPrahovani(pocetFiltrovanychPasiem,k) = signalPoRozkladeWT(pocetFiltrovanychPasiem,k)- (prah^2/signalPoRozkladeWT(pocetFiltrovanychPasiem,k));                
            end
        else
            signalPoPrahovani(pocetFiltrovanychPasiem,k) = 0;
        end

    end
end
% pridanie nedotknutych pasiem signalPoRozkladeWT do premennej signalPoPrahovani podľa toho koľko pásiem bolo prahovaných
if (pocetFiltrovanychPasiem < 2)
    signalPoPrahovani(2,:) = signalPoRozkladeWT(2,:);
end

if (pocetFiltrovanychPasiem < 3)
    signalPoPrahovani(3,:) = signalPoRozkladeWT(3,:);
end
signalPoPrahovani(4,:) = signalPoRozkladeWT(4,:);
signalPoPrahovani(5,:) = signalPoRozkladeWT(5,:);

if vykreslenie == 1  % Vykreslenie pasiem po PRAHOVANI
    figure()
    subplot 511
    plot(signalPoPrahovani(1,:))
    title('Frekvencne pasmo 62.5 - 125 Hz')
    subplot 512
    plot(signalPoPrahovani(2,:))
    title('Frekvencne pasmo 31.25 - 62.5 Hz')
    subplot 513
    plot(signalPoPrahovani(3,:))
    title('Frekvencne pasmo 15.625 - 31.125 Hz')
    subplot 514
    plot(signalPoPrahovani(4,:))
    title('Frekvencne pasmo 7.78125 - 15.5625 Hz')
    subplot 515
    plot(signalPoPrahovani(5,:))
    title('Frekvencne pasmo 0 - 7.78125 Hz')
end

% pre vypocet SNR priamo pouzi ako uzitocny signal vystup zo spätnej vlnkovej transformacie
uzitocnySignal = iswt(signalPoPrahovani, vlnka);    % vystup zo spatnej vlnkovej transformacie
rusenie = signalEKGSRusenim - uzitocnySignal;       % hodnota rusenia pre vypocet SNR

hodnotySNR = [];
% prechadzaj v cykle s krokom dlzkaOkna vzorky uzitocneho signalu a sumu a
% vypocitaj z nich SNR pomocou funkcie VypocetSNR
indexacia = 0; % premenna potrebna pre indexaciu, potreba brat napr. pri dlzkaOkna = 250 : 1 až 250 vzorkov, 251 až 500 vzorkov, atd.
aktualneOkno = 1;                           % zacni na prvom okne
for t = 1:dlzkaOkna:length(uzitocnySignal)
    hodnotySNR = [hodnotySNR, VypocetSNR(uzitocnySignal(1, indexacia+1:aktualneOkno*dlzkaOkna), rusenie(1, indexacia+1:aktualneOkno*dlzkaOkna))];
    indexacia = indexacia+dlzkaOkna;
    aktualneOkno = aktualneOkno + 1;
end
% z pola SNR vytvor median, ktory bude dalej porovnavany s vygenerovanym SNR
medianSNR = median(hodnotySNR);
  
  if vykreslenie == 1 % Vykreslenie vstupneho a uzitocnej zlozky pre porovnanie
    figure()
    subplot 211
    plot(signalEKGSRusenim)
    title('Vstupny vygenerovany signal EKG so sumom')
    subplot 212
    plot(uzitocnySignal)
    title('Uzitocna zlozka z celkoveho vygenerovaneho signalu ')
  end 
% uloz a uzavri txt subor 
obsahSuboru = obsahSuboru + num2str(medianSNR) + "\t";
end
end
end
obsahSuboru = obsahSuboru + "\n\n";
end
fprintf(fid, strrep(obsahSuboru,".",","));
fclose(fid);
end