% 混沌時間序列的 Volterra 自適應預測(一步預測) -- 主函數% 參考文獻：朱昀.水聲信號非線性分析方法研究.西北工業大學博士學位論文.2002.06% P28-38clcclearclose all% dx/dt = sigma*(y(t)-x(t))% dy/dt = r*x(t) - y(t) - x(t)*z(t)% dz/dt = -b*z(t) + x(t)*y(t)% y = [-1;0;1];                   % 起始點 (3x1 的列向量)% h = 0.01;                       % 積分時間步長% 積分方程用四階 Runge-Kutta 法積分方程組得數值解% 數據除去前面 8000 個點，用后面 5000 個點load('LorenzData.mat')            % 加數 Lorenz 離散數據 x,y,z,共 5000 個點x = x(1:500);X = [x-mean(x)]/[max(x)-min(x)];  % 歸一化到均值為 0，振幅為 1Times = 5e+3       % 最小二乘迭代次數tau = 10           % 重構時延m = 3              % 重構維數p = 3              % Volterra 級數階數%--------------------------------------------------% 混沌序列的相空間重構 (phase space reconstruction)[xn,dn] = PhaSpaRecon(X,tau,m);% 輸入參數：    s          混沌序列%               tau        重構時延%               m          重構維數% 輸出參數：    xn         相空間中的點序列(每一列為相空間中一個點)%               dn         一步預測的目標%------------------------------------------------------interval = 4;[xn_train,dn_train,xn_test,dn_test] = TrainTestSample(interval,xn,dn);% 相空間中點的軌跡分解成：訓練樣本與測試樣本% 在所有樣本中，每隔 interval 個樣本為訓練樣本，其余為測試樣本%--------------------------------------------------% 訓練與測試[Hn,len_Hn] = Volterra_train(xn_train,dn_train,p,Times);% 輸入參數：    xn_train   訓練樣本(每一列為一個樣本)%               dn_train   訓練目標%               p          Volterra 級數階數%               Times      最小二乘估計迭代次數% 輸出參數：    Hn         最小二乘估計濾波器權矢量 Hn% 　　　　　    len_Hn　　 權矢量長度dn_pred = Volterra_test(xn_test,p,Hn);% 輸入參數：    xn_test    測試樣本%               p          Volterra 級數階數%               Hn         最小二乘估計濾波器權矢量 Hn% 輸出參數：    dn_pred　　一步預測值err = dn_test - dn_pred;err_mse = sqrt(sum(err.^2)/length(err))          % 預測的均方誤差n = 1:length(dn_test);plot(n,dn_test,'r+:',n,dn_pred,'bo:')title('+為真值，o為預測值')