%skyhawk&flyinghawk
clear all;
data=load('bk1.txt');%讀取數(shù)據(jù)
whl=data(:,4);
[whsl,lllll]=size(whl);
m=6;
P=26;  %平均循環(huán)周期

N=80;%最多可預(yù)測(cè)步數(shù)估計(jì)值

for i=1:whsl
    whlsj(i)=whl(i);
end

[lmd_m,idx,min_d,idx1,min_d1]=lyapunov(m,whlsj,whsl,P);%求lyapunov指數(shù)
lmd_m
% t_m=fix(1/lmd_m)+1       %最大預(yù)測(cè)步數(shù)
t_m=1
% for j=1:(whsl-t_m)
%     whlsj(j)=whl(j);
% end   
for i=(whsl-N+2-t_m):(whsl-N+1)         %預(yù)測(cè)后t_m步
    [yc,y1(i),y2(i)]=jiaquanyijie(m,whlsj,i-1);%計(jì)算第i步預(yù)測(cè)值 
    whlsj(i)=yc;  %將第i步預(yù)測(cè)值作為完好率數(shù)據(jù)的第i個(gè)值進(jìn)行下一步預(yù)測(cè)
end

y(whsl-N+1)=yc;
fch=(y(whsl-N+1)-whl(whsl-N+1))*(y(whsl-N+1)-whl(whsl-N+1));
shuliang=1;

for i=(whsl-N+2):(whsl)
    whlsj(i-t_m)=whl(i-t_m);  %換為實(shí)際值
    [y(i),y1(i),y2(i)]=jiaquanyijie(m,whlsj,i-1);
    whlsj(i)=y(i);           %換為預(yù)測(cè)值
    
    fch=fch+(y(i)-whl(i))*(y(i)-whl(i));
    shuliang=shuliang+1;
end

fch=sqrt(fch)/shuliang

yyy=[whl,y'];
save('bkycqush.txt','yyy','-ASCII');

kk=1:whsl;
plot(kk,whl,'b',kk,y,'r')