%ucitavanje EMG signala
load gk002;
s1=gk002(:,1);      %1. EMG signal
s2=gk002(:,2);      %2. EMG signal
s3=gk002(:,3);      %3. EMG signal
kut=gk002(:,4);     %4. signal kut
[N,d]=size(gk002);

%prikaz EMG signala
figure;
Y=[1:N];
subplot(4,1,1); plot(Y,s1);title('Prvi kanal');
hold;
%stem(tren,ones(1,length(tren)),'-xw');
subplot(4,1,2); plot(Y,s2);title('Drugi kanal');
hold
%stem(tren,ones(1,length(tren)),'-xw');
subplot(4,1,3); plot(Y,s3);title('Treci kanal');
hold
%stem(tren,ones(1,length(tren)),'-xw');
subplot(4,1,4); plot(Y,kut);title('Kut');



%odreðivanje vektora trenutaka koji odvajaju pojedine intervale;
M=700;                                  
num=ones(1,M);                          
kut=abs(kut);
y=filter(num,1,kut);                    %primjena filtra usrednjavanja duljine M=700 uzoraka na signal kuta
y=y/450;        
H=abs(fftshift(fft(y)));                %FFT signala kuta
H2=H(round(length(H)/2):length(H));     %uzimamo samo polovicu simetricnog spektra
y2=y(350:length(y));                    %filtar usrednjavanja unosi kasnjenje od M/2, pa pomicemo signal
y3=y2-min(y2(5000:65000))+1;            %y3 je razlika pocetnog signala kuta te pomaknute aproksimacije
for i=1:length(y3)
if (y3(i)>150)||(y3(i)<0)               %ogranicavanje signala na 0<y3<150
y3(i)=0;
end
end
tren=[];                                %vektor koji sadrzi vremenske trenutke razdvajanja segmenata
podrucja=zeros(length(y3),1);                  %vektor duljine pocetnog signala u kojem se oznacavaju trenutci koji   
for i=2:length(y3)                      %razdvajaju pojedine segmente.
    if(y3(i)>0)&&(y3(i-1)==0)
      j=0;
      while (y3(i+j)>0)&&(i+j<length(y3)-100)
          j=j+1;
      end
      if j>100                          %ako je detektirani interval dovoljno dug, zabiljezi novi trenutak
          podrucja(i+round(j/2))=1;            %u vektorima t i tren
          tren=[tren i+round(j/2)];
      end
    end
end
t2=zeros(length(y3),1);                  
for i=1:length(tren)-1
    int=round((tren(i+1)-tren(i))*0.1);
    t2(tren(i)+int:tren(i+1)-int)=1;        %t2 se sastoji od niza vremenskih prozora koji cine 80% od ukupnih intervala kontrakcija
end;
figure;
plot(kut);title('KUT - izvorni signal i usrednjeni signal, te prikaz trenutaka koji odvajaju pojedine odsjecke');hold on;plot(y2,'r');hold on; plot(podrucja*3000,'k');
%plot(kut); hold; plot(t*3000,'k');

save('podatci.mat','s1','s2','s3','kut','t2','tren');
pause
clear all;
close all;
%vektor trenutaka koji razdvajaju pojedine kontrakcije - APROKSIMACIJA "OD OKA"
%tren=[3000 6000 9250 12700 16350 20100 23300 26500 29600 33000 36600 40400 44400 47800 51600 55500 59200 62200 65500];


load Podatci.mat;
emg(:,1)=s1-mean(s1); emg(:,2)=s2-mean(s2); emg(:,3)=s3-mean(s3);
Nsr=3500;               %prosjecna perioda kontrakcije (=(ukupno trajanje)/(broj perioda))
L=ceil(0.1*Nsr);        %definiranje duljine vremenskog prozora (10% od periode)
if mod(L,2)~=1          %ako je parna, ucini ju neparnom
    L=L+1;
end
W=hamming(L);           %Hammingov vremenski otvor
R=20;                   
prek=L-R;               %koliko uzoraka se poklapaju dva susjedna vremenska otvora
MA=ceil(0.1*Nsr/R);     %moving average, duljina filtra;
N=1024;                 %broj frekvencija FT
fs=1000;                %frekvencija uzorkovanja 
n=length(emg(:,1));                         %duljina signala
brstup=fix((n-prek)/(length(W)-prek));      %broj stupaca u matrici B, duljini vektora vremena
brred=N/2+1;                                %broj redova u matrici B, odgovaraju pojedinim frekvencijama
U=sum(W.^2);                                


for j=1:3
    [B, f, t]=specgram(emg(:,j),N,fs,W,prek);
    for i=1:brstup
        PSD(:,i)=B(:,i).*conj(B(:,i));
        PSD(:,i)=PSD(:,i)./U;
        k=1;
        suma=sum(PSD(1,i));
        pola=sum(PSD(:,i))/2;
        while (suma<pola)
            k=k+1;
            suma=sum(PSD(1:k,i));
        end
        median(i,j)=f(k);
    end
    fmed(:,j)=filtfilt(ones(1,MA)/MA,1,median(:,j));
end

for j=1:3 
    for k=1:length(t)-1
        if t2(k*20) 
            o(k*20)=fmed(k,j);
        else
            o(20*k)=0;
        end
    end
    for i=1:length(tren)-1
        medmax(i,j,1:2)=0;
        for k=tren(i):tren(i+1)
            if o(k)>medmax(i,j,1)
                medmax(i,j,1)=o(k);           
            	medmax(i,j,2)=k;
            end
        end
    end
    yj=medmax(:,j,1);
    xj=medmax(:,j,2);
    n=length(medmax(:,j,1));
    a(j)=(sum(xj.*yj)/n-mean(xj)*mean(yj))/(mean(xj.^2)-mean(xj)^2);
    b(j)=mean(yj)-a(j)*mean(xj);
end
x=[1:length(emg(:,1))];
for i=1:3
    figure;
    plot(a(i)*x+b(i)); hold; plot(medmax(:,i,2),medmax(:,i,1),'r');
    axis([0 length(emg(:,i)) 0 round(max(medmax(:,i,1)))+10]);
    text(10000,20, ['a(' num2str(i) ')= ' num2str(a(i)*60*10000, '%6.2f') ' Hz/min       MSE= ' num2str(sum((a(i)*medmax(:,i,2)+b(i)-medmax(:,i,1)).^2)/(length(medmax(:,i,1))))]);
end

