% G-P 算法同時(shí)求關(guān)聯(lián)維和Kolmogorov熵 (輸入時(shí)間序列數(shù)據(jù))% 使用平臺(tái) - Matlab7.0% 作者：陸振波，海軍工程大學(xué)% 歡迎同行來信交流與合作，更多文章與程序下載請(qǐng)?jiān)L問我的個(gè)人主頁% 電子郵件：luzhenbo@yahoo.com.cn% 個(gè)人主頁：http://luzhenbo.88uu.com.cn% 參考文獻(xiàn): 趙貴兵,石炎福,段文峰等.從混沌序列同時(shí)計(jì)算關(guān)聯(lián)維和Kolmogorov熵[J].計(jì)算物理,1999;16(5):309~315clcclear allclose all%--------------------------------------------------------------------------% 產(chǎn)生 Lorenz 時(shí)間序列% dx/dt = sigma*(y-x)% dy/dt = r*x - y - x*z% dz/dt = -b*z + x*ysigma = 16;             % Lorenz方程參數(shù)r = 45.92;               b = 4;          y = [-1;0;1];           % 起始點(diǎn) (3x1 的列向量)h = 0.01;               % 積分時(shí)間步長(zhǎng)k1 = 30000;             % 前面的迭代點(diǎn)數(shù)k2 = 5000;              % 后面的迭代點(diǎn)數(shù)X = LorenzData(y,h,k1+k2,sigma,r,b);X = X(k1+1:end,1);      % 時(shí)間序列(列向量)%--------------------------------------------------------------------------% G-P算法計(jì)算關(guān)聯(lián)維rr = 0.5;Log2R = -6:rr:0;        % log2(r)R = 2.^(Log2R);fs = 1/h;               % 信號(hào)采樣頻率t = 10;                 % 時(shí)延dd = 4;                 % 嵌入維間隔D = 4:dd:36;            % 嵌入維    p = 50;                 % 限制短暫分離，大于序列平均周期(不考慮該因素時(shí) p = 1)  tic Log2Cr = log2(CorrelationIntegral(X,t,D,R,p));   % 輸出每一行對(duì)應(yīng)一個(gè)嵌入維toc%--------------------------------------------------------------------------% 結(jié)果作圖figureplot(Log2R,Log2Cr','k.-'); axis tight; hold on; grid on;xlabel('log2(r)'); ylabel('log2(C(r))');title(['Lorenz, length = ',num2str(k2)]);%--------------------------------------------------------------------------% 最小二乘擬合Linear = [6:10];                            % 線性似合區(qū)域[a,B] = LM2(Log2R,Log2Cr,Linear);           % 最小二乘求斜率和截距disp(sprintf('Correlation Dimension = %.4f',a));for i = 1:length(D)    Y = polyval([a,B(i)],Log2R(Linear),1);    plot(Log2R(Linear),Y,'r');endhold off;%--------------------------------------------------------------------------% 求梯度Slope = diff(Log2Cr,1,2)/rr;                % 求梯度xSlope = Log2R(1:end-1);                    % 梯度所對(duì)應(yīng)的log2(r)figure;plot(xSlope,Slope','k.-'); axis tight; grid on;xlabel('log2(r)'); ylabel('Slope');title(['Lorenz, length = ',num2str(k2)]);%--------------------------------------------------------------------------% 差分求K熵KE = -diff(B)/(dd*t)*fs*log(2);             % 用采樣頻率 fs 和公式 log(x) = log2(x)*log(2) 將單位轉(zhuǎn)化成 nats/sD_KE = D(1:end-1);                          % K熵所對(duì)應(yīng)的嵌入維figure;plot(D_KE,KE,'k.-'); grid on; hold on;xlabel('m'); ylabel('Kolmogorov Entropy (nats/s)');title(['Lorenz, length = ',num2str(k2)]);%--------------------------------------------------------------------------% 輸出顯示disp(sprintf('Kolmogorov Entropy = %.4f',min(KE)));