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