%      sigma = 10;          % Lorenz方程參數
%      r = 28;              
%      b = 8/3;          
% 
%      y = [-1;0;1];         % 起始點 (3x1 的列向量)
%      h = 0.01;             % 積分時間步長
% 
%      k1 =8000;           % 前面的迭代點數10000
%      k2 =5000;             % 后面的迭代點數5000
% 
%      z = LorenzData(y,h,k1+k2,sigma,r,b);
%      X = z(k1+1:end,1);

    clc
%     clear
    load data13;
    X=data13;
%     X=denosie1(X);%去除噪聲
%     X=denosie2(X);%再次去噪
%     tau=Mutual_Information_main(X);%互信息法求延遲時間tau
% plot(X)
tau=6;
    X = normalize_1(X);  % 歸一化
    data=X';
    logdelt = 0.2;
      % delt=std(X);       %計算標準差,求r的范圍
       ln_r = [-7:logdelt:0];
       %r=[delt./2:0.004:2*delt];
       delt = exp(ln_r);
       K1=[];
       K1(1)=0;
       K13=[];
      K13(1)=0;
       GP=[];
       GP(1)=0;
m=1;
        for k=1:length(ln_r)
           r=delt(k); 
           C(k)=correlation_interal(m,data,r,tau);%  輸出變量為關聯積分
           m,k
          if (C(k)<0.0001)
             C(k)=0.0001;
          end
         ln_C(k)=log(C(k));%lnC(r)
        % ln_r(k)=log(r(k));%lnr
       end   %k循環完畢
       
       ln_Cr=ln_C;
       ln_r=ln_r;
       LinearZone = [10:25];
       
       Xj=ln_r(LinearZone);
       Yj=ln_Cr(LinearZone);
       X_mean=mean(Xj);
       Y_mean=mean(Yj);
%        for j=1:length(LinearZone)
       l_xy=sum((Xj-X_mean).*(Yj-Y_mean));
       l_xx=sum((Xj-X_mean).*(Xj-X_mean));
       a=l_xy./l_xx;
       b1=Y_mean-a*X_mean;
       K13(m)=b1./(m*tau);
       
for m=2:20
       % for k=1:length(r)
       for k=1:length(ln_r)
           r=delt(k); 
           C(k)=correlation_interal(m,data,r,tau);%  輸出變量為關聯積分
           m,k
          if (C(k)<0.0001)
             C(k)=0.0001;
          end
         ln_C(k)=log(C(k));%lnC(r)
        % ln_r(k)=log(r(k));%lnr
       end   %k循環完畢
         C;
          %------------------------------------------------------
           % 擬合線性區域
       ln_Cr=ln_C;
       ln_r=ln_r;
       LinearZone = [10:25];
       F = polyfit(ln_r(LinearZone),ln_Cr(LinearZone),1);
       GP(m)=F(1);  K1(m)=F(2);
       K1(m)=K1(m)/(m*tau);
      
      
       Xj=ln_r(LinearZone);
       Yj=ln_Cr(LinearZone);
       X_mean=mean(Xj);
       Y_mean=mean(Yj);
%        for j=1:length(LinearZone)
       l_xy=sum((Xj-X_mean).*(Yj-Y_mean));
       l_xx=sum((Xj-X_mean).*(Xj-X_mean));
       a=l_xy./l_xx;
       b=Y_mean-a*X_mean;
       K13(m)=(b1-b)./(tau);
       b1=b;
       
%       if (abs((GP(m)-GP(m-1)))<0.01)
%             break;
%       end  
    
%        subplot(411)
%      plot(ln_r,ln_C,'+:');grid on;
%      xlabel('ln r'); ylabel('ln C(r)');
%      hold on;
end

       subplot(211)%畫出K熵隨嵌入維數m變化的圖,
       X=[1:1:m];
       Y1=K1;
       plot(X,Y1,':m*'); 
       axis([2 20 0 0.9])
       hold on;grid on;
       xlabel('m'); ylabel('K1');
       
%      subplot(421)%畫出關聯積分隨嵌入維數m變化的圖,
%      X=[1:1:m];
%      Y2=GP;
%      plot(X,Y2,':m*'); 
%       title('The plot lnC(r) vs.lnr');
%     %legend('關聯維數隨 m 變化最終趨于不變',2)
%     hold on;
%     grid on;
%     xlabel('m'); ylabel('關聯維數 D(m) ');
       
     subplot(212)%畫出K熵隨嵌入維數m變化的圖,
     X=[1:1:m];
     Y3=abs(K13);
     plot(X,Y3,':m*'); 
     axis([6 20 0 0.2])
     hold on;grid on;
     xlabel('m'); ylabel('K2');