%^^^^^^^^^^^^^^^^^^G--P法求時間序列的最佳嵌入維數(shù)主函數(shù)^^^^^^^^^^^^^^^^^^^^^
clear all;
clc;
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
r0=5e-2;            % attentian 
rn=1;

nn=100;            %控制軌道的變化步長0-10，200
hr=(rn-r0)/nn;
r=r0:hr:(rn-hr);  

tic
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
load('D:\researchdata\S.mat');
h=1*1e-12;

xn=S(1,50001:55000);
maxS=max(xn);
xn=xn./maxS;
% xxn=xn(50001:55000);
N=length(xn);
St=1/h;
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^^^^^^^時間序列的功率譜分析^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Y=fft(xn);
Nn=length(Y);
L=round(Nn/2);
Y(1)=[];
power=abs(Y(1:L)).^2;
nyquist=1/2;
freq=nyquist*(1:L)/L*St;
period=1./freq;
[mp,index]=max(power);
avrpa=period(index)/h;    %關(guān)于軌道的平均周期不是時間
avrpa=round(avrpa)
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
J=jjautoc2(xn);                       %時序最佳時延，調(diào)用函數(shù)獲得
J
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
m0=2; mn=15;            %^^^^^^^^^^^最大嵌如維數(shù)^^^^^^^^^^^^^^^^^^^^^^^^^^
mm=mn-m0+1;
rrr=zeros(mm,nn);
dm=zeros(1,mm);
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
for m=m0:mn
    m
    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    M=N-(m-1)*J;
    xnxn=zeros(M,m);
    for i=1:M
        j=1:m;
        xnxn(i,:)=xn(i+(j-1)*J);
    end
    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    distance=zeros(M,M);   
    for i=1:M
        for j=1:M
            distance(i,j)=norm((xnxn(i,:)-xnxn(j,:)),2);
        end
    end
    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    distance2=zeros(M,M);
    ooo=zeros(M,M);
    for ij=1:nn
        distance2=r(1,ij)-distance;
        ooo=distance2>=0;
        rrr(m-m0+1,ij)=sum(sum(ooo))/(M*M);
    end
    clear distance distance2 xnxn
    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
end              
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
r=log(r);
rrr=log(rrr);

figure(3);
for m=1:(mn-m0+1)
    plot(r,rrr(m,:));hold on;
end
hold off;
xlabel('lnr');ylabel('lnrrr');
%^^^^注意:不同的潛入為數(shù)時，擬合的范圍不同，要根據(jù)圖形3的結(jié)果確定^^^^^^^^^^^^
for i=1:mm
    [hk,hj1]=min(abs(rrr(i,:)-(-3.5)));    % attentian  below limit
    [hk,hj2]=min(abs(rrr(i,:)-(-1.5)));    % attentian  up limit
    p1=polyfit(r(hj1:hj2),rrr(i,hj1:hj2),1);
    dm(i)=p1(1);
end
%^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
figure(4);
mx=m0:mn;
plot(mx,dm);xlabel('m');ylabel('dm');
dm
toc