
clear all
clc
close all
%část algoritmu převzato z "Pokročilá analýza biologických signálů (FACS) 14. Jiří Kozumplík, EMD (Empirical Mode Decomposition)"
%isimf a ismonotonic jsou pomocné funkce
 hodnoty=[];
 ts=5000/500;
% signaly=load('signaly_f_ADP9_2');
% signal=signaly.signaly_f_ADP9_2;
 signaly=load('signaly_zkraceno');
 signal=signaly.signaly_zkraceno;
for i=1:1
hodnoty=[];
puvodni=signal(i,:);
x=signal(i,:);
x = x - mean(x);
puvodni=puvodni - mean(puvodni);
sig=puvodni;

fvz=500;
n = 1:1:length(x); 
A=40;% amplituda brumu
brum =A*cos(2*pi*50*n/fvz);
x=x+brum;
zasumeny=brum;

while ~ismonotonic(x)
    q=x;
sd=Inf;
while (sd > 0.03) | ~isimf(q)  % (sd > 0.1)

[y1,x1] = findpeaks(q);

%nalezení minim
[s,Spozice] = findpeaks(-q);
s=-s;
N=length(q);
xx = 1:N;
konec=q(1,end);
p = spline([0 x1 N+1],[0 y1 0],xx);
d = spline([0 Spozice N+1],[0 s 0],xx);

rozdil = q-(p+d)/2;
sd = sum((q-rozdil).^2)/sum(q.^2);
q=rozdil;

end
 
x=x-q;
hodnoty=[hodnoty;q];
 
end
hodnoty=[hodnoty;x];
%zpětná rekonstrukce signálu s brumem
 l=size(hodnoty);
 
 %zobrazení první imf složky 
 figure(10)
 plot(hodnoty(1,:))
 title('první IMF složka')
 %spektrum první imf složky
 figure(11) 
 SPIMF1=fft(hodnoty(1,:));
 SPIMF1=SPIMF1(1,1:1000);
 f=linspace(0,500,5000);
 f=f(1,1:1000);
 plot(f,abs(SPIMF1))
 title('Spektrum první IMF složky')

 qrs=hodnoty(1,:)>(A+1);
 w=0;
 for o=2:l(1,1)
 w=w+hodnoty(o,:)+qrs;
 end
 %kresba rozlozeneho signalu
hodnoty=[zasumeny;hodnoty];
t = 1:length(x); 
 figure(4);  
nimf = size(hodnoty,1); 
 for (m=1:nimf)
  subplot(nimf,1,m); plot(t,hodnoty(m,:),'k','LineWIdth',1.5);
   ylabel('C')
 end
 subplot(nimf,1,1); title('IMFs');
 subplot(nimf,1,nimf);xlabel('počet vzorků[-]');
 
y=w;
zac=500;
kon=4500;
rozdil=[];
for d=zac:kon
    rozdil=[rozdil w(1,d)-puvodni(1,d-1)];
end
%kresba rozdílového signálu rekonstruované a původního signálu
 figure(3)
 subplot(2,1,1)
 rozdilovy=[];
 plot(y(zac:kon))
 hold on
 plot(sig(zac:kon))
 ylabel('U[\muV]')
 title('fitrovaný a původní signál')
 xlabel('n[-]')
 subplot(2,1,2)
 plot(rozdil)
 ylabel('U[\muV]')
 title('rozdílový signál')
 xlabel('n[-]')
end


