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 是否可以验证以上程序正确呢？ 其中AutoCorrelation(X,maxLags,IsPlot)函数如下：
function = 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
回复 #30 winterdij 的帖子
大家都是同命中人啊，俺也得在十月分左右把先前做的东西总结分析，想把这部分东西加进去，就看出来出不来了 革命尚未成功同志仍需努力啊^_^不管用什么方法 时间延迟是最终要的也是必须的啊 G-P算法计算时就用到了时间延迟 这样肯定行不通了
回复 #31 winterdij 的帖子
你的程序我算得Lorenz的时间延迟是12单纯地从像空间重构来说，自相关法也是值得肯定的，也比较成熟，它认为时间延迟和嵌入维数的选取是互相独立的（相空间重构奠基人Takens认为是独立的）。吕金虎论文及书、相关期刊都说此方法使用于线性线性相关性，对于高维系统不太适合，对此我也没有理解太深
回复 #35 sssssxxxxx921 的帖子
我今天这么做了一下，你看看是否可以。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;
= PhaSpaRecon(X,tau,m);
xn = xn'; % 重构相空间点集，每一行为一个点
r = 0.04; % r 相对于吸引子直径的大小,(0,1)之间
%---------------------------------------------------
% 调用 mex 函数
Dim = CorrelationDimension(xn,r)
————————————————————————————————————————
其中normalize_1函数如下：
function = normalize_1(sig_input)
% 信号归一化到均值为 0,振幅为 1
% = normalize_sig(sig_input)
% 输入参数：sig_input输入信号(可以批处理)
% 输出参数：sig_output 标准化的信号
= 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 = PhaSpaRecon(s,tau,m)
% 混沌序列的相空间重构 (phase space reconstruction)
% = PhaSpaRecon(s, tau, m)
% 输入参数： s 混沌序列(列向量)
% tau 重构时延
% m 重构维数
% 输出参数： xn 相空间中的点序列(每一列为一个点)
% dn 一步预测的目标(行向量)
= 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
回复 #35 sssssxxxxx921 的帖子
我给你个我的qq号26156847回复 #37 winterdij 的帖子
这几天有事，你的程序验证没，最大李雅普诺夫指数出来和权威结果对不回复 #1 winterdij 的帖子
时间延迟是指delt_s的第一个极小值所对应的t再乘以你的采样时间嵌入维数是由s_cor的最小值所对应的t乘以采样时间得到的延时时间窗口再除以时间延迟再加1所得 tau应该是一个正整数么？
平均周期p应该是一个正整数么，如果高频的离散随机信号的话，p肯定小于1，那怎么算？
回复 #40 JulianChin 的帖子
如果p小于1，就取整数1了吧！ 原帖由 octopussheng 于 2007-10-19 15:27 发表 http://www.chinavib.com/forum/images/common/back.gif如果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(), grid on
ylabel('Power')
xlabel('Period(Years/Cycle)')
这种方法求出的平均频率只是真实频率的分数倍数吧，把真实频率看作1，求出的平均频率可能是零点几,
平均频率到底是真实的频率还是分频？？ 用CC算法计算延迟时间和窗口时，是限定嵌入维数在2-5之间的吧，为什么会求出的tw大于tau的四倍呢？
好资料
这个真是一个好资料，我找了很久了，谢谢啊 回复 winterdij 的帖子delt_s的第一个极小值对应的时间t就是延迟时间td，
s_cor的最小值对应的就是时间窗tw，
然后利用公式tw=（m-1）*td，这样就求出了嵌入维数m