main_volterra.m

时间序列分析G-P method to calculate the correlation dimension of Matlab (Mex version) not use "temporary s

% 混沌时间序列的 volterra 预测(一步预测) -- 主函数
% 使用平台 - Matlab7.0
% 作者：陆振波，海军工程大学
% 欢迎同行来信交流与合作，更多文章与程序下载请访问我的个人主页
% 电子邮件：[e-mail]
% 个人主页：http://luzhenbo.88uu.com.cn

clc
clear
close all

%--------------------------------------------------
% 公共参数
k1 = 6000;              % 前面的迭代点数
k2 = 2000;              % 后面的迭代点数 (总样本数)

%--------------------------------------------------
% 产生 Lorenz 时间序列
% dx/dt = sigma*(y-x)
% dy/dt = r*x - y - x*z
% dz/dt = -b*z + x*y

sigma = 16;          % Lorenz方程参数
r = 45.92;               
b = 4;          

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

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

Z = LorenzData(y,h,k1+k2,sigma,r,b);
X = Z(k1+1:end,1);   % 时间序列(列向量)
        
%--------------------------------------------------

X = normalize_a(X,1);   % 信号归一化到均值为0,振幅为1
tau = 1;               % 时延
m = 3;                  % 嵌入维数
p = 3;                  % Volterra阶数

%--------------------------------------------------

train_num = 500        % 训练样本数
test_num = 1500         % 测试样本数

%--------------------------------------------------
% 混沌序列的相空间重构 (phase space reconstruction)

x_train = X(1:train_num);
x_test =  X(train_num+1:train_num+test_num);

[xn_train,dn_train] = PhaSpaRecon(x_train,tau,m);
[xn_test,dn_test] = PhaSpaRecon(x_test,tau,m);

%------------------------------------------------------
% 训练与预测

% 训练算法求权矢量
[Wn,err_mse1] = volterra_train_lu(xn_train,dn_train,p);
err_mse1_std = err_mse1/var(X)

% 测试样本的一步预测
dn_pred = volterra_test(xn_test,p,Wn);
err2 = dn_test - dn_pred;
err_mse2 = sum(err2.^2)/length(err2);
err_mse2_std =  err_mse2/var(X)

%----------------------------------------------------
% 结果作图

figure;
subplot(211);
plot(1:length(err2),dn_test,'r+:',1:length(err2),dn_pred,'bo-');
title('真实值(+)与预测值(o)')
subplot(212);
plot(err2,'k');
title('预测绝对误差')