奇异谱分析MATLAB代码.zip

奇异谱分析MATLAB代码

function [drvr]=ssa（x1L）

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% -----------------------------------------------------------------                           
%    Author: Francisco Javier Alonso Sanchez    e-mail:fjas@unex.es
%    Departament of Electronics and Electromecanical Engineering
%    Industrial Engineering School
%    University of Extremadura
%    Badajoz
%    Spain
% -----------------------------------------------------------------
%
% SSA generates a trayectory matrix X from the original series x1
% by sliding a window of length L. The trayectory matrix is aproximated 
% using Singular Value Decomposition. The last step reconstructs
% the series from the aproximated trayectory matrix. The SSA applications
% include smoothing filtering and trend extraction.
% The algorithm used is described in detail in: Golyandina N. Nekrutkin 
% V. Zhigljavsky A. 2001. Analisys of Time Series Structure - SSA and 
% Related Techniques. Chapman & Hall/CR.

% x1 Original time series （column vector form）
% L  Window length
% y  Reconstructed time series 重构时间序列
% r  Residual time series r=x1-y 剩余时间序列
% vr Relative value of the norm of the approximated trajectory matrix with respect
%	  to the original trajectory matrix 相对于原始轨迹矩阵的近似轨道矩阵范数的相对值

% The program output is the Singular Spectrum of x1 （must be a column vector）
% using a window length L. You must choose the components be used to reconstruct 
%the series in the form [i1i2:ik...iL] based on the Singular Spectrum appearance.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%




% Step1 : Build trayectory matrix

   N=length（x1）; 
   if L>N/2;L=N-L;end
	K=N-L+1; 
   X=zeros（LK）;  
	for i=1:K
	  X（1:Li）=x1（i:L+i-1）; 
	end
    
% Step 2: SVD奇异值分解

   S=X*X‘; %
	[Uautoval]=eig（S）;%eig返回矩阵的特征值和特征向量，U是特征向量，autoval是特征值
	[di]=sort（-diag（autoval））;
    %x= diag（vk）以向量v的元素作为矩阵X的第k条对角线元素，当k=0时，v为X的主对角线；当k>0时，v为上方第k条对角线；当k<0时，v为下方第k条对角线。
    %sort（A，dim）表示对A中的各行元素升序排列，sort（A）默认的是升序，d是排序好的向量，i是向量d中每一项对应于A中项的索引
   d=-d;
   U=U（:i）;
   sev=sum（d）; %特征值的和
	plot（（d./sev）*100）hold onplot（（d./sev）*100‘rx‘）;
	title（‘Singular Spectrum‘）;xlabel（‘Eigenvalue Number‘）;ylabel（‘Eigenvalue （% Norm of trajectory matrix retained）‘）  %轨迹矩阵的范数
   V=（X‘）*U; 
   rc=U*V‘;
    %奇异值分解S=U*D*V
% Step 3: Grouping

   I=input（‘Choose the agrupation of components to reconstruct the series in the form I=[i1i2:ik...iL]  ‘）
   Vt=V‘;
   rca=U（:I）*Vt（I:）;

% Step 4: Reconstruction

   y=zeros（N1）;  
   Lp=min（LK）;
   Kp=max（LK）;

   for k=0:Lp-2
     for m=1:k+1;
      y（k+1）=y（k+1）+（1/（k+1））*rca（mk-m+2）;
     end
   end

   for k=Lp-1:Kp-1
     for m=1:Lp;
      y（k+1）=y（k+1）+（1/（Lp））*rca（mk-m+2）;
     end
   end

   for k=Kp:N
      for m=k-Kp+2:N-Kp+1;
       y（k+1）=y（k+1）+（1/（N-k））*rca（mk-m+2）;
     end
   end

   figure;subplot（21

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3427  2019-05-08 21:54  ssa.m