!------------------ 奇异谱分析 ---------------------------------!
!要求输入的有:                                                  !
!    1.X(n):需作奇异谱分析的序列;                               !
!    2.m:嵌套维数;                                              !
!    3.itp:输出典型波形的个数(一般应取得稍大);                  !
!输出的有:(存放在filename文件中)                                !
!    1.特征值及其方差贡献;                                      !
!    2.典型波形(特征向量);(每列)                                !
!    3.时间系数;(每列)                                          !
!    4.重建序列;(每列)                                          !
!    5.相邻两个典型波形的滞后相关;                              !
!        第一列为TPC1和TPC2的相关;第二列为TPC2和TPC3的相关;     !
!        第三列为TPC3和TPC4的相关; ...                          !
!注意:                                                          !
!    1.第25,26行的输出格式与itp的取值有关;                      !
!    2.本程序用到了smaths.lib中的函数;                          !
!    3.本程序用到了sub_standard.f90和sub_correlation.f90文件;   !
!--------------------------------------- 程正泉 2000.5 ---------!
subroutine sub_ssa(n,m,itp,X,filename)
implicit none
integer,intent(in)::n,m,itp
real,dimension(n),intent(in)::X
real,allocatable,dimension(:,:)::raX,raEgvt,raTime,Rc,R
integer::i1,i2
character*50 filename
101 format(6i15)                    ! 格式随itp的取值而变动
102 format(i5,6f15.4)               ! 格式随itp的取值而变动

!----- 1.组成时间后延矩阵raX--------
allocate(raX(m,n-m+1))
do i1=1,m; do i2=1,n-m+1
   raX(i1,i2)=X(i2+i1-1)
enddo;enddo

!----- 2.EOF分解---------------------
allocate(raEgvt(m,m),raTime(m,n-m+1))
call sub_ssa_eof(m,n-m+1,itp,1,raX,raEgvt,raTime,filename)

!----- 3.计算重建序列----------------
allocate(Rc(m,n))
call sub_ssa_Rc(n,m,raEgvt,raTime,Rc)

!----- 4.计算相邻两个典型玻形之间的后延相关
allocate(R(itp-1,-50:50))
do i1=1,itp-1
   call sub_two_backxg(n-m+1,50,raTime(i1,1:n-m+1),raTime(i1+1,1:n-m+1),R(i1,-50:50))
enddo

!----- 5.输出
open(90,file=filename,access='append')
write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'
write(90,*)'典型波形:每一列为一个典型波形'
write(90,101)(i1,i1=1,itp)
do i1=1,m
   write(90,102)i1,(raEgvt(i1,i2),i2=1,itp)
enddo

write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'
write(90,*)'时间系数:每一列为对应典型波形的时间系数'
write(90,101)(i1,i1=1,itp)
do i2=1,n-m+1
   write(90,102)i2,(raTime(i1,i2),i1=1,itp)
enddo

write(90,'(/a28/)')'重建序列,每一列为一个重建序列'
do i2=1,n
   write(90,102)i2,(Rc(i1,i2),i1=1,itp)
enddo

write(90,'(/a48)')'相邻两个典型波形之间的后延相关,注意最后一列无效.'
do i2=-50,50
   write(90,102)i2,(R(i1,i2),i1=1,itp-1),0.0
enddo

deallocate(raX,raEgvt,raTime,Rc,R)
end

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
subroutine sub_ssa_eof(m,n,itp,ki,raX,raEgvt,raTime,filename)
implicit none
integer,intent(in)::m,n,itp,ki
real,dimension(m,n),intent(in)::raX
real,dimension(m,m),intent(out)::raEgvt,raTime(m,n)
real,allocatable,dimension(:,:)::raA,raEgvl
real::lamda,temp
integer::i,j
character*50 filename

open(90,file=filename)

!----- 0.资料处理 -----------------------!
call sub_standard2(m,n,raX,ki)

!----- 1.计算协方差矩阵raA=raX*raX' -----!
allocate(raA(m,m))
call mxytf(m,n,raX,m,              &
             m,n,raX,m,            &
		    m,m,raA,m)
do i=1,m; do j=1,m
   raA(i,j)=raA(i,j)/float(n)
enddo; enddo

!----- 2.计算raA的特征值和特征向量 ------!
allocate(raEgvl(m,4))
call evcsf(m,raA,m,raEgvl(1:m,1),raEgvt,m)
deallocate(raA)

!----- 3.计算时间系数raTime=raEgvt'*raX -!
call mxtyf(m,m,raEgvt,m,           &          
             m,n,raX,m,            &
			   m,n,raTime,m)

!----- 4.验证分解是否正确 ---------------!
write(90,*)'当下面两列值相等或非常近似时,SSA分解正确;否则,分解不正确'
write(90,'(/a8,6x,a12,7x,a22)')'序号','特征值','时间系数平方累加的平均'
do i=1,m
   lamda=0.0
   do j=1,n
      lamda=lamda+raTime(i,j)**2/float(n)
   enddo
   write(90,'(i7,2f20.5)')i,raEgvl(i,1),lamda
enddo
write(90,'(/a43)')'!-----!-----!-----!-----!-----!-----!-----!'

!----- 5.输出EOF分解的结果 --------------!
write(90,'(a4,8x,a6,6x,a10,4x,a8,4x,a12)')'序号','特征值','累计特征值','方差贡献','累计方差贡献'
temp=0.0
do i=1,m
   temp=temp+raEgvl(i,1)
   raEgvl(i,2)=temp
enddo
do i=1,m
   raEgvl(i,3)=raEgvl(i,1)/raEgvl(m,2)   
enddo
do i=1,m
   raEgvl(i,4)=0.0
   do j=1,i
      raEgvl(i,4)=raEgvl(i,4)+raEgvl(j,3)
   enddo
enddo
do i=1,m
   write(90,'(i3,f15.3,f15.3,f12.3,f14.3)')i,(raEgvl(i,j),j=1,4)
enddo
deallocate(raEgvl)
close(90)
end

!-----------------------------------------------


!----利用SSA重建m个Rc序列(每一列为一个Rc序列)--------
subroutine sub_ssa_Rc(n,m,raVector,raTime,Rccs)
implicit none
integer,intent(in)::n,m
real,dimension(m,m),intent(in)::raVector,raTime(m,n-m+1)
real,dimension(m,n),intent(out)::Rccs
integer::i,j,k
real::tmp
do k=1,m
do i=1,n
  if((i >= 1) .and. (i <= m-1))then
    tmp=0.0
	do j=1,i
      tmp=tmp+raVector(j,k)*raTime(k,i-j+1)
    enddo
    Rccs(k,i)=tmp/real(i)
  elseif((i >= m) .and. (i <= n-m+1))then
    tmp=0.0
    do j=1,m
      tmp=tmp+raVector(j,k)*raTime(k,i-j+1)
    enddo
    Rccs(k,i)=tmp/real(m)
  elseif((i >= n-m+2) .and. (i <= n))then
    tmp=0.0
    do j=i-n+m,m
      tmp=tmp+raVector(j,k)*raTime(k,i-j+1)
    enddo
    Rccs(k,i)=tmp/real(n-i+1)
  endif
enddo;enddo
end