%% metoda EMD
% Petr Procházka
close all

%% výběr signálu, svodu

svod=1;
cas= 0.002:0.002:10;               % Prepocet casove osy
signal=xT(svod,1:5000);           % Výběr svodu a úseku signálu
puvodni=signal;                         % Záloha původního signálu


%% výpočet EMD
%matice, do kterých jsou v každém cyklu ukládány IMF, spektra IMF a zbytkový signál
matice_IMF=[puvodni];
matice_signal=[puvodni];
matice_spektrum=abs(fft(puvodni));

[nic, extremy]=findpeaks(signal);   % zjištění počtu extrémů v signále
while length(extremy)>1

odchylka=1;                                     % Pomocná hodnota, aby pro 1. průběh cyklu while
while odchylka>0.3

% OBÁLKA SIGNÁLU

%obalka maxim
[peak,pozice]=findpeaks(signal);                             % Nalezení souřadnicvšech lokálních maxim
pozice=pozice/500;                                                   % Přepočet jednotlivých vzorků na sekundy
obalka_max=spline([cas(1) pozice cas(end)],[signal(1) peak signal(end) ],cas);   % Vytvoření obálky

%obalka minim
signal1=0-signal;                                           % Přetočení podél x-ové osy pro výpočet obálky minim
[peak1,pozice1]=findpeaks(signal1);         % Nalezení souřadnicvšech lokálních minim
peak2=0-peak1;                                           % Zpětné přetočení hodnoty lokálních minim podél osy x
pozice1=pozice1/500;                                 % Přepočet jednotlivých vzorků na sekundy
obalka_min=spline([cas(1) pozice1 cas(end)],[signal(1) peak2 signal(end) ],cas);  % Vytvoření obálky

%pruměr obálek
prumer=(obalka_max+obalka_min)/2;

%odečtení průměru od signálu
IMF=signal-prumer;

odchylka=sum(((signal-IMF).^2)/((signal).^2));      % výpočet  odchylky

signal=IMF;         %výstup vnitřního cyklu je přejmenován na signál protože je přiven na vstup vnitřního cyklu v následujícím kroku
end

% rozhodovací kritérium, kdy je ukončen proces empirického rozkladu. 
%kritérium je počítáno ze spekter jednotlivých vnitřních funkcí (IMF).
spektrum=abs(fft(IMF)); %počítání spektra danné IMF

%% podmínka pro ukončení rozkladu metodou IMF
if (max(spektrum(1:10))/10)>(max(spektrum(20:2500)))
break
end

signal=matice_signal(end,:) - IMF;    
    
matice_IMF=[matice_IMF;IMF];
matice_spektrum=[matice_spektrum; abs(fft(IMF))];
matice_signal=[matice_signal;signal];
[nic, extremy]=findpeaks(signal);
end


%% VYKRESLENÍ
% obecné vykreslení jednotlivých IMF (vnitřních funkcí)
[radky,sloupce]=size(matice_IMF);           % zjištění počtu IMF

%pokud je počet IMF větší než zadaná hodnota (nastaveno na 7), dojde k vykreslení vnitřních
%funkcí a zbytku do více oken
if radky>7
radky1=ceil((radky+1)/2)            % dojde k rozdělení subplotů do jednotlivých oken rovným dílem 
vykresleni=radky1*100+11;        % výpočet obecného označení prvního subplotu                  

   figure(1)
     subplot(vykresleni)
     plot(cas,matice_IMF(1,:))      %vykreslení původního signálu
     ylabel('signal')
     vykresleni=vykresleni+1;         % hodnota (označení) subplotu se zvýší o 1
     hold on
for j=2:radky1
     subplot(vykresleni)
     plot(cas,matice_IMF(j,:))      %vykreslení danné vnitřní funkce
     hodnota=[ 'IMF' num2str(j)]
     ylabel(hodnota)
     hold on
    vykresleni=vykresleni+1;        % hodnota (označení) subplotu se zvýší o 1
end
vykresleni=vykresleni-radky1;    % hodnota (označení) subplotu se zvýší o 1

for j=radky1+1:radky
    figure(2)
    subplot(vykresleni)
    plot(cas,matice_IMF(j,:))       %vykreslení danné vnitřní funkce
     ylabel('IMF')
    hold on
    vykresleni=vykresleni+1;
end
    subplot(vykresleni)                 % hodnota (označení) subplotu se zvýší o 1
    plot(cas,matice_signal(end,:),'r')              % vykreslení zbytku
    ylabel('zbytek')
    xlabel('čas [s]')
else
  vykresleni=(radky+1)*100+11;   % výpočet obecného označení prvního subplotu
     figure(1)
     subplot(vykresleni)
     plot(cas,matice_IMF(1,:))      %vykreslení původního signálu
     ylabel('signal')
     title('Rozklad signálu na jednotlivé vlastní modální funkce IMFa zbytek, pomocí metody EMD')
     vykresleni=vykresleni+1;         % hodnota (označení) subplotu se zvýší o 1
     hold on
    
for j=2:radky
     subplot(vykresleni)
     plot(cas,matice_IMF(j,:))      %vykreslení danné vnitřní funkce
     hodnota=[ 'IMF' num2str(j+1)]
     ylabel(hodnota)
     hold on
     vykresleni=vykresleni+1;         % hodnota (označení) subplotu se zvýší o 1
end
    subplot(vykresleni)
    plot(cas,matice_signal(end,:),'r')              % vykreslení zbytku
    ylabel('zbytek')
    xlabel('čas [s]')
end

vyfiltrovany=sum(matice_IMF(2:end,:));

% vykreslení původního signálu bez navázaného trendu a vyfiltrovaného signálu
figure
subplot(311)
plot(cas, xZ(svod,:))
title('originální signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(312)
plot(cas, puvodni)
title('původní signál s navázaným trendem')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(313)
plot(cas, vyfiltrovany,'r')
title('"vyfiltrovaný" signál vzniklý součtem IMF')
xlabel('Čas [s]')
ylabel('napětí [uV]')

figure
subplot(211);
plot(cas, xZ(svod,:))
hold on
plot(cas, vyfiltrovany,'r')
title('originální a "vyfiltrovaný" signál vzniklý součtem IMF')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(212);
plot(cas,(xZ(svod,:))-vyfiltrovany );
title('rozdílový signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')

%% vykreslení spekter

spekt_puv=abs(fft(puvodni));
spekt_vyf=abs(fft(vyfiltrovany));
spekt_orig=abs(fft(xZ(svod,:)));

delka_vyrezu_spektra=50;
figure
subplot(131)
stem((1:delka_vyrezu_spektra)/10,spekt_puv(1:delka_vyrezu_spektra))
title('                  výřez spektra původního signálu')
xlabel('frekvence [Hz]')
ylabel('amplituda[uV/Hz]')
axis([0 5 0 max(spekt_puv(1:delka_vyrezu_spektra))+10000])
subplot(132)
stem((1:delka_vyrezu_spektra)/10,spekt_puv(1:delka_vyrezu_spektra))
title('                  výřez spektra původního signálu')
xlabel('frekvence [Hz]')
ylabel('amplituda[uV/Hz]')
axis([0 5 0 max(spekt_vyf(1:delka_vyrezu_spektra))+10000])
subplot(133)
stem((1:delka_vyrezu_spektra)/10,spekt_vyf(1:delka_vyrezu_spektra),'r')
title('                  výřez spektra vyfiltrovaného signálu')
xlabel('frekvence [Hz]')
ylabel('amplituda[uV/Hz]')
axis([0 5 0 max(spekt_vyf(1:delka_vyrezu_spektra))+10000])


delka_vyrezu_spektra1=100;
figure
subplot(211)
stem((1:delka_vyrezu_spektra1)/10,spekt_puv(1:delka_vyrezu_spektra1))
hold on
stem((1:delka_vyrezu_spektra1)/10,spekt_vyf(1:delka_vyrezu_spektra1),'r')
title('výřez spektra vyfiltrovaného a původního signálu s trendem')
xlabel('frekvence [Hz]')
ylabel('amplituda[uV/Hz]')
subplot(212)
stem((1:delka_vyrezu_spektra1)/10,spekt_orig(1:delka_vyrezu_spektra1))
hold on
stem((1:delka_vyrezu_spektra1)/10,spekt_vyf(1:delka_vyrezu_spektra1),'r')
title('výřez spektra originálního a vyfiltrovaného signálu')
xlabel('frekvence [Hz]')
ylabel('amplituda[uV/Hz]')

%% vypočet SNR
orig=xZ(svod,:);
VykonSig=sum(orig.^2);
VykonSum=sum((vyfiltrovany-orig).^2);
SNR= 10*log( VykonSig/VykonSum );

SNR=SNR-5
