%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%never change a running system%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
clc

%Eingabe der Parameter
Abtastfrequenz=1080;                   %Abtastfrequenz in Hz
h=108;                                  %Glättungsfenster ca. 100ms
Prokraft=0.1;                          %10 Prozent der maximalen Kraft
tp=20;                                 %Tiefpass Grenzfrequenz
hp=500;                                %Hochpass Grenzfrequenz
Wn=[tp/Abtastfrequenz hp/Abtastfrequenz];%Bandpassfilter 20Hz und 500 Hz
%Ende der Eingabe der Parameter

%Lädt die Bodenreaktionskräfte Ein 
p1=load('Platte2.txt');
F1=abs(p1);

%Lädt alle EMG datein Ein
lade=load('emg12.txt');

%Ende Ladevorgang

%%%%%%%%%%%%%%%
[b,a]=butter(3,Wn,'bandpass');
%

%Bestimmung der Zeitskala%%
%begin=1/Abtastfrequenz;
begin=0;
l=length(lade);                        %Alle EMG dateien habe die gleiche Länge
f=l/Abtastfrequenz;
t=linspace(begin,f,l);                 %Definition eines Vektors für die Zeitskala

%%Signalverarbeitung%%
emg_1=detrend(lade);                  %Entfernt jegliche DC offset vom EMG Signal
emg1=abs(emg_1);                       %Gleichrichten des EMG Signals

%Filtern
d_1=filtfilt(b,a,emg1);                %Anwendung des Bandpassfilter auf emg1
e_1=abs(d_1);

%RMS Glättung
EMGrms_1= sqrt( e_smooth(e_1.*e_1, h, 'SumFilter') );%Ruft ein Funktion auf,durch der dass EMG Signal glättet wird
%ISEMG_2= ( e_smooth(e_1, h, 'SumFilter') );  %IEMG Glättung Gleitender
                                             %Mittelwert
mean(EMGrms_1) ;
max(EMGrms_1);
%Ende Gläattung                                      

%%Amplituden-Normalisierung%%
maxi=max(EMGrms_1);                    %Such das Maximum im geglätteten EMG Signal
EMG_1=linspace(1,1,l);                 %Definition eines Vektor mit der Länge L
k=0;
for i=1:l
    k=k+1;
    EMG_1(k)=(EMGrms_1(i)/maxi)*100;   %Die normalisierte Werte in einen neuen Vektor schreiben
end

%%Stützphase definieren%%
fmax=max(F1);                          %Maximum der Kraft definieren
fl=fmax*Prokraft;                      %Prozent der maximalen Kraft
p=find(F1>=fl);                        %Finde 10% der maximalen Kraft
p1=p(1);                               %Erster Punkt
p2=p(end);                             %Letzter Punkt


G=F1((p1-216):p2);
fmax1=max(G);                           %Maximum der Kraft definieren
fl1=fmax*Prokraft;                      %Prozent der maximalen Kraft
pk=find(G>=fl);                        %Finde 10% der maximalen Kraft
pk1=pk(1);                               %Erster Punkt
pk2=pk(end);                             %Letzter Punkt

%Integral
% z=integrate(EMGrms_1,Abtastfrequenz,p1,p2);
% s=sum(z);
% round(s)

%figure,plot(s)
 
%Ende Integral

%Ausschnitt der Stützphase
ab=EMG_1((p1-216):p2);
l1=length(ab);
f1=l1/Abtastfrequenz;
t1=linspace(0,f1,l1);
%figure,plot(t1,ab)
m=max(ab);
% hold on
% plot([t1(pk1) t1(pk1)],[0 m],'k','LineWidth',2)
% plot([t1(pk2) t1(pk2)],[0 m],'k','LineWidth',5)
%axis('tight')

maxiaus=max(ab);                   %Such das Maximum im geglätteten EMG Signal
EMGaus_1=linspace(1,1,l1);                 %Definition eines Vektor mit der Länge L
k=0;
for i=1:l1
    k=k+1;
    EMGaus_1(k)=(ab(i)/maxiaus)*100;   %Die normalisierte Werte in einen neuen Vektor schreiben
end
%figure,plot(t1,EMGaus_1)
%%%%%
X=EMGaus_1;
%-6dB und -20dB Zeichen
pro30=100/33;
pro60=100/66;

dreizig=ones(1,length(X))*max(X)/pro30;
sixty=ones(1,length(X))*max(X)/pro60;
%Schnittpunkte suchen
k=0; 
n=0;
for i=1:length(X)-1,
    if X(i) <= max(X)/pro30 && X(i+1) >= max(X)/pro30,
        k=k+1;
        cut_sixdb(k)=i;
    elseif X(i) >= max(X)/pro30 && X(i+1) <= max(X)/pro30,
        k=k+1;
        cut_sixdb(k)=i;
    end
    if X(i) <= max(X)/pro60 && X(i+1) >= max(X)/pro60,
        n=n+1;
        cut_twentydb(n)=i;
    elseif X(i) >= max(X)/pro60&& X(i+1) <= max(X)/pro60,
        n=n+1;
        cut_twentydb(n)=i;
    end
end
%END mit suchen
%%%%Case abfrage für 66% bis 100%
l33=length(cut_sixdb);
l66=length(cut_twentydb);
[pos,wert] = max(X);
%j=0;
%Azu=0;
% switch l33
%     case l33==2
%        cut_sixdb(1);
%        cut_sixdb(end);
%     case l33>=3
%         
%   for i=1:l66
%    if l66==l66;
%        
%    elseif cut_twentydb(i)< wert && cut_twentydb(i+1) > wert
%         Azu=cut_twentydb(i)
%    end
%    
%   end
% end

for i=1:l66
    if i==l66
    elseif cut_twentydb(i)< wert && cut_twentydb(i+1) > wert
  Azu=cut_twentydb(i);
    end
end
%%%Von 100 bis 66
for i=1:l66
    if i==l66
    else cut_twentydb(i)< wert && cut_twentydb(i+1)> wert
  Bab=cut_twentydb(i+1);
    end
end
Aab=NaN;
% 30 bis 60
for i=1:l33
    if i==l33
    elseif cut_sixdb(i) >= cut_twentydb(2) && cut_sixdb(i+1) <= cut_twentydb(2)
  Aab=cut_sixdb(i)
    end
end
%Von 60 bis 30
Bzu=NaN;
for i=1:l33
    if i==l33
    else cut_sixdb(i)<= cut_twentydb(2) && cut_sixdb(i+1)>= cut_twentydb(2)
  Bzu=cut_sixdb(i+1)
    end
end

%gewünschte länge erzeugen ohne Inhalt
i_sixdb(1:length(X))=NaN;
i_twentydb(1:length(X))=NaN;
%schnittpunkte im Vektor schreiben
i_sixdb(cut_sixdb)=dreizig(cut_sixdb);
i_twentydb(cut_twentydb)=sixty(cut_twentydb);
figure,plot(t1,EMGaus_1,'k')
hold on
plot(t1,dreizig,'r')
plot(t1,sixty,'g')
plot(t1,i_sixdb,'r.','Markersize',15)
plot(t1,i_twentydb,'r.','Markersize',15)
%%%

for xx=1:cut_sixdb(1)
    plot([t1(xx) t1(xx)],[0 X(xx)],'y')
    hold on
end
hold on
for xx=cut_sixdb(end):length(X)
    %ui=X(xx)
    plot([t1(xx) t1(xx)],[0 X(xx)],'y')
    hold on
end
hold on
for xx=cut_sixdb:cut_twentydb
    %ui=X(xx)
    plot([t1(xx) t1(xx)],[0 X(xx)])
    hold on
end

for xx=cut_twentydb(2):cut_sixdb(end)
    %ui=X(xx)
    plot([t1(xx) t1(xx)],[0 X(xx)],'g')
    hold on
end
hold on
for xx=Azu:wert
    %ui=X(xx)
    plot([t1(xx) t1(xx)],[0 X(xx)],'r')
    hold on
end
hold on
for xx=wert:Bab
    %ui=X(xx)
    plot([t1(xx) t1(xx)],[0 X(xx)],' m')
    hold on
end

hold on

m=max(ab);
plot([t1(pk1) t1(pk1)],[0 m],'k','LineWidth',2)
plot([t1(pk2) t1(pk2)],[0 m],'k','LineWidth',1)
hold on
plot(t1(wert),X(wert),'r.','Markersize',15);
%text(t1(wert),X(wert),['Mittenfreq. =',Mittenfrequenz1,'Hz']);
%realtive Brandbreite ist Quotient aus Bandbreite und der Mittelfreq.
%%%%%TEST OnSet Bestimmung
% figure(3)
% [d,c]=butter(5,10/1080,'low');
% 
% d_3=filtfilt(d,c,emg1);                %Anwendung des Bandpassfilter auf emg1
% d_2=abs(d_3);
% onset_time=onset(d_2,600,900,900,1500)
%figure,plot(d_2)
%%%%
% for xx=cut_sixdb:1:cut_twentydb
%     plot([X(xx) X(xx)],[0 60])
%     hold on
% end
    
%%%%
%%Plotten der EMG Analyse%%

% figure, plot(lade,'k')
% title ('Rohsignal')
% axis('tight')
% figure,plot(emg1)
% axis('tight')
% title('Gleichgerichtet')

% figure, subplot(2,1,1)
% plot(emg1,'g')
% axis('tight')
% title('Gleichgerichtet')
% subplot(2,1,2)
% plot(e_1,'r')
% axis('tight')
% title('Gefiltert')
% legend('Bandpassfilter 20-500Hz')
% saveas(gcf,'Bandpassfilter','bmp')
%on=mean(emg1)+std(emg1)
% on=6.4719;
% figure,plot(t,emg1,'g')
% hold on
% plot(on,'r.','Markersize',15)
% axis('tight')
%figure,plot(t,e_1,'g')
%title('smoothing')
% plot(on,'r.','Markersize',15)
% % hold on
% plot(t,EMGrms_1,'g')
% legend('Gefiltert','RMS Glättung')
% saveas(gcf,'RMS','bmp')

% figure,subplot(2,1,1)
% plot(t,e_1,'g')
% title('smoothing')
% hold on
% plot(t,ISEMG_2,'r')
% axis('tight')
% legend('Gefiltert','ISEMG Glättung')
% subplot(2,1,2)
% plot(t,e_1,'g')
% title('smoothing')
% hold on
% plot(t,EMGrms_1,'r')
% axis('tight')
% legend('Gefiltert','RMS Glättung')
% saveas(gcf,'RMS und ISEMG','bmp')

figure,subplot(2,1,1)
plot(t,F1)
hold on
plot([t(p1) t(p1)],[0 800],'k')
plot([t(p2) t(p2)],[0 800],'k')
axis('tight')
title('Bodenreaktionskraft in Z-Richtung')
hold off
subplot(2,1,2)
plot(t,EMG_1,'r')
hold on
plot([t(p1) t(p1)],[0 80],'k','LineWidth',2)
plot([t(p2) t(p2)],[0 80],'k','LineWidth',2)
axis('tight')
title('Normalisiert und Zeitpunkte definiert')
hold off

figure,plot(t1,EMGaus_1,'k')
hold on
plot(t1,dreizig,'r')
plot(t1,sixty,'g')
plot(t1,i_sixdb,'r.','Markersize',15)
plot(t1,i_twentydb,'r.','Markersize',15)
saveas(gcf,'Schnittpunkte','jpg')