load pomiary; %ladowanie danych pomiarowych
s=length(pomiary);
nr=5; %nalezy podac numer stanowiska
randn('state',nr); %powtarzalnosc
r=0.06*randn(1,s);
p=pomiary+r;
p=p(1:10:s);r=1.0*r(1:10:s);
t=linspace(0,300,s/10+1);
figure(1);
plot(t,p); %charakterystyka skokowa zdjeta na obiekcie

%Aproksymacja Strejca
figure(2);
k=1.0672;T3=22 ;l1=1; m1=[T3 1]; %k=1.0672; T3=13 dla n=5; T3=22 dla n=3
n=3; %rzad aproksymacji
li=1;mi=[T3 1];
for i=1:n-1
   [li,mi]=series(li,mi,l1,m1);
end
printsys(k*li,mi);
y3=step(k*li,mi,t);
e3=p-y3';
plot(t,p,t,y3,'-m',t,e3,':',t,r);xlabel('Czas, s');ylabel('Wyjscie ');grid on;
title('Aproksymacja Strejca');
legend('pomiary - p','aproksymacja - y3','odchylenie - e3');
J3=sum(e3.^2), %kryterium jakosci aproksymacji
%Aproksymacja Kupfmullera (model inercyjny I-go rzedu z opoznieniem) 
figure(3);
Ti=44;Toi=25;  %Ti=44 ; Toi=25;
li=k;mi=[Ti 1];
[lp1,mp1]=pade(Toi,5);
[lo1,mo1]=series(lp1,mp1,li,mi);
y1=step(lo1,mo1,t);
e1=p-y1';
plot(t,p,t,y1,'m',t,e1);xlabel('Czas, s');ylabel('Wyjscie ');grid on;
J1=sum(e1.^2); %kryterium jakoœci aproksymacji

%Aproksymacja Kupfmullera (model inercyjny II-go rzedu z opoznieniem) 
figure(4);
mj=20;
TT2=linspace(10.0,50,mj);
for j=1:mj 
  T2=TT2(j);
  T1=22.4;To2=12.8; %T1=22.4; To2=13;T2=31.73;T2=31.73;
l2=k;m2=[T1*T2 T1+T2 1];
[lp2,mp2]=pade(To2,5);
[lo2,mo2]=series(lp2,mp2,l2,m2);
y2=step(lo2,mo2,t);
e2=p-y2';
plot(t,p,t,y2,'m',t,e2,':');xlabel('Czas, s');ylabel('Wyjscie ');grid on;hold on;
J2=sum(e2.^2); %kryterium jakosci aproksymacji
%y4=k*(1-exp(-stepfun(t,Toi).*(t-Toi)/Ti));
y4=k*(1-T1/(T1-T2)*exp(-stepfun(t,To2).*(t-To2)/T1)+T2/(T1-T2)*exp(-stepfun(t,To2).*(t-To2)/T2));
plot(t,p,t,y4);
e4=p-y4;
J4=sum(e4.^2);
[T2,J4];
df(j,1)=T2;df(j,2)=J4;
end
T22=df(:,1);J44=df(:,2);
figure(5);
plot(T22,J44);grid on;hold on;axis([10 50 0 5]);
xlabel('T2');ylabel('J4');title('Kryterium');
pp=p';
save pom.dat pp -ascii