疑问：关于CC 方法计算时间延迟和嵌入维数

winterdij

我使用自相关法计算时间延迟 计算得出Lorenz 时间序列的时间延迟为13.
Logistic 混沌时间序列的时间延迟为1.
其程序如下：
clc
clear
close all

%---------------------------------------------------
% 产生 Lorenz 时间序列

sigma = 16;             % Lorenz 方程参数 a = 16 | 10
b = 4;                  %                 b = 4 | 8/3
r = 45.92;              %                 c = 45.92 | 28        

y = [-1,0,1];       % 起始点 (1 x 3 的行向量)
h = 0.01;           % 积分时间步长

k1 = 30000;          % 前面的迭代点数
k2 = 10000;          % 后面的迭代点数

z = LorenzData(y,h,k1+k2,sigma,r,b);
X = z(k1+1:end,1);

%---------------------------------------------------
% 自相关法 (直接求 tau)
% 自相关函数下降到初始值的 1-1/e 时的 tau 即为所求 (tau 从 1 开始)

maxLags = 100;
IsPlot = 1;
tau_AutoCorrelation = AutoCorrelation(X,maxLags,IsPlot);
tau_AutoCorrelation

winterdij

是否可以验证以上程序正确呢？

winterdij

其中AutoCorrelation(X,maxLags,IsPlot)函数如下：
function [tau] = AutoCorrelation(X,maxLags,IsPlot)
% 自相关法求混沌时间序列重构的时间延迟
% 输入参数：X         混沌时间序列
%           maxLags　 最大时间延迟
% 输出参数：tau　　   时间延迟　

ACF = autocorr(X,maxLags);   % 求自相关函数
ACF = ACF(2:end);            % 将第一个 tau=0 的情况去掉   

% 自相关函数下降到初始值的 1-1/e 时的 tau 即为所求 (tau 从 1 开始)
gate = (1-exp(-1))*ACF(1);
temp = find(ACF<=gate);
if (isempty(temp))
    disp('err: max delay time is too small!')
    tau = [];
else
    tau = temp(1);   
end

if IsPlot
    figure;
    plot(1:length(ACF),ACF)
    xlabel('Lag');
    title('自相关法求时延');
end

sssssxxxxx921

大家都是同命中人啊，俺也得在十月分左右把先前做的东西总结分析，想把这部分东西加进去，就看出来出不来了   革命尚未成功  同志仍需努力啊  ^_^
不管用什么方法   时间延迟是最终要的  也是必须的啊   G-P算法计算时就用到了时间延迟 这样肯定行不通了

sssssxxxxx921

你的程序我算得Lorenz的时间延迟是12
单纯地从像空间重构来说，自相关法也是值得肯定的，也比较成熟，它认为时间延迟和嵌入维数的选取是互相独立的（相空间重构奠基人Takens认为是独立的）。吕金虎论文及书、相关期刊都说此方法使用于线性线性相关性，对于高维系统不太适合，对此我也没有理解太深

winterdij

我今天这么做了一下，你看看是否可以。
1.使用自相关法计算时间延迟。
2.使用Takens方法计算关联维数，这里需要嵌入维数m，我是从1开始给m赋值的，直到计算出的关联维数d不满足下面关系为止。
m>=2d+1
这样就找到了嵌入维数m，你试试，看看是否可以。
Takens方程如下：
clc
clear
close all

%---------------------------------------------------
% 产生 Lorenz 时间序列

sigma = 10;          % Lorenz方程参数
r = 28;               
b = 8/3;         

y = [-1;0;1];        % 起始点 (3x1 的列向量)
h = 0.01;            % 积分时间步长

k1 = 6000;           % 前面的过渡点数
k2 = 2000;           % 后面的迭代点数

z = LorenzData(y,h,k1+k2,sigma,r,b);    % 用四阶 Runge-Kutta 法产生 k1+k2 个点
X = z(k1+1:end,1);                      % 去掉前面 k1 个过渡点

X = normalize_1(X);  % 归一化到均值为 0，振幅为 1

%---------------------------------------------------
% 相空间重构

tau = 14;
m = 3;
[xn] = PhaSpaRecon(X,tau,m);

xn = xn';       % 重构相空间点集，每一行为一个点
r = 0.04;       % r 相对于吸引子直径的大小,(0,1)之间

%---------------------------------------------------
% 调用 mex 函数

Dim = CorrelationDimension(xn,r)


————————————————————————————————————————
其中normalize_1函数如下：
function [sig_output] = normalize_1(sig_input)
% 信号归一化到均值为 0,振幅为 1
% [sig_output] = normalize_sig(sig_input)
% 输入参数：sig_input  输入信号(可以批处理)
% 输出参数：sig_output 标准化的信号

[rows,cols] = size(sig_input);
if (rows ==1)
    sig_input = sig_input';
    len = cols;
    num = 1;
else
    len = rows;
    num = cols;
end

mean_sig = mean(sig_input);
sig_input = sig_input - repmat(mean_sig,len,1);  % 0 均值

dis = max(sig_input) - min(sig_input);
sig_output = sig_input./repmat(dis,len,1);       % 振幅为 1
————————————————————————————————————————
其中PhaSpaRecon函数如下：
function [xn,dn] = PhaSpaRecon(s,tau,m)
% 混沌序列的相空间重构 (phase space reconstruction)
% [xn, dn, xn_cols] = PhaSpaRecon(s, tau, m)
% 输入参数：    s          混沌序列(列向量)
%               tau        重构时延
%               m          重构维数
% 输出参数：    xn         相空间中的点序列(每一列为一个点)
%               dn         一步预测的目标(行向量)

[rows,cols] = size(s);
if (rows>cols)
    len = rows;
    s = s';
else
    len = cols;
end

if (nargout==2)
   
    if (len-1-(m-1)*tau < 1)
        disp('err: delay time or the embedding dimension is too large!')
        xn = [];
        dn = [];
    else
        xn = zeros(m,len-1-(m-1)*tau);
        for i = 1:m
            xn(i,:) = s(1+(i-1)*tau : len-1-(m-i)*tau);   % 相空间重构，每一行为一个点
        end
        dn = s(2+(m-1)*tau : end);    % 预测的目标
    end
   
elseif (nargout==1)
   
    if (len-1-(m-1)*tau < 0)
        disp('err: delay time or the embedding dimension is too large!')
        xn = [];
    else
        xn = zeros(m,len-(m-1)*tau);
        for i = 1:m
            xn(i,:) = s(1+(i-1)*tau : len-(m-i)*tau);   % 相空间重构，每一行为一个点
        end
    end
   
end

winterdij

我给你个我的qq号26156847

sssssxxxxx921

这几天有事，你的程序验证没，最大李雅普诺夫指数出来和权威结果对不

nzddqr

时间延迟是指delt_s的第一个极小值所对应的t再乘以你的采样时间
嵌入维数是由s_cor的最小值所对应的t乘以采样时间得到的延时时间窗口再除以时间延迟再加1所得

JulianChin

tau应该是一个正整数么？

平均周期p应该是一个正整数么，如果高频的离散随机信号的话，p肯定小于1，那怎么算？

octopussheng

如果p小于1，就取整数1了吧！

JulianChin

原帖由 octopussheng 于 2007-10-19 15:27 发表

                               
登录/注册后可看大图

如果p小于1，就取整数1了吧！

load sunspot.dat
year = sunspot(:,1);
wolfer = sunspot(:,2);
plot(year,wolfer)
title('Sunspot Data')
Y = fft(wolfer);
N = length(Y);
Y(1) = [];
power = abs(Y(1:N/2)).^2;
nyquist = 1/2;
freq = (1:N/2)/(N/2)*nyquist;
plot(freq,power), grid on
xlabel('cycles/year')
title('Periodogram')
period = 1./freq;
plot(period,power), axis([0 40 0 2e7]), grid on
ylabel('Power')
xlabel('Period(Years/Cycle)')


这种方法求出的平均频率只是真实频率的分数倍数吧，把真实频率看作1，求出的平均频率可能是零点几,

平均频率到底是真实的频率还是分频？？

waanm

用CC算法计算延迟时间和窗口时，是限定嵌入维数在2-5之间的吧，为什么会求出的tw大于tau的四倍呢？

shenyong3938021

这个真是一个好资料，我找了很久了，谢谢啊

cqupenghao

回复 winterdij 的帖子


    delt_s的第一个极小值对应的时间t就是延迟时间td，
s_cor的最小值对应的就是时间窗tw，
然后利用公式tw=（m-1）*td，这样就求出了嵌入维数m