!------------多通道奇异谱分析(MSSA)方法子程序--------------------------!
!输入的有:                                                             !
!   1. L:通道数目;                                                     !
!   2. n:序列长度;                                                     !
!   3. m:滞后时间(嵌套维数);                                           !
!   4. itp:输出典型波形的波形(一般应取稍大);                           !
!   5. ki:对原始场的处理(ki= -1,原始场;ki= 0,距平场;ki= 1,标准化场)    !
!   6. X(L,n):由L个通道序列组成的矩阵;                                 !
!输出的有(存放于FileName中):                                           !
!   1. 特征值、累计特征值、方差贡献、累计方差贡献;                     !
!   2. 典型波形(特征向量);                                             !
!   3. 时间系数;(计算成对TPC之间的最大后延相关可求得显著耦合振荡信号)  !
!   4. 重建序列;(计算成对RCCS之间的最大后延相关可求得各耦合信号位相差) !
!注意:                                                                 !
!   1. 第26行的输出格式随itp而变;                                      !
!   2. 本程序用到了fortran 5.0的smaths.lib中的数学函数;                !
!   3. 本程序用到了sub_standard.f90文件;                               !
!-------------------------------------- 程正泉 2000.5 -----------------!
subroutine sub_mssa(L,n,m,itp,ki,X,FileName)
implicit none
integer,intent(in)::L,n,m,itp,ki
real,dimension(L,n),intent(in)::X
real,allocatable,dimension(:,:)::raX,raEgvt,raTime,Rccs(:,:,:)
integer::i,j,k
character*50 FileName
100 format(1x,i4,16f10.5)           ! 格式随itp的取值而变动

!----- 1.组成后延矩阵
allocate(raX(L*m,n-m+1))
do i=1,L*m; do j=1,n-m+1
  k=(i-1)/m+1
  raX(i,j)=X(k,j+(i-(k-1)*m)-1)
enddo;enddo

!----- 2.进行EOF分解
allocate(raEgvt(L*m,L*m),raTime(L*m,n-m+1))
call sub_mssa_eof(L*m,n-m+1,itp,ki,raX,raEgvt,raTime,FileName)

!----- 3.计算重建序列
allocate(Rccs(L*m,n,L))
call sub_mssa_Rc(n,m,l,raEgvt,raTime,Rccs)

!----- 4.输出
open(99,file=FileName,access='append')
do k=1,L
   write(99,'(/a2,i2,a35)')'第',k,'个通道的典型波形,每列为一个典型波形'
   do i=(k-1)*m+1,(k-1)*m+m
      write(99,100)i-(k-1)*m,(raEgvt(i,j),j=1,itp)
   enddo
enddo
do k=1,L
   write(99,'(/a2,i2,a45)')'第',k,'个通道的时间系数,每列为对应典型波形的时间系数'
   do j=1,n-m+1
      write(99,100)j,(raTime(i,j),i=1,itp)
   enddo
enddo
do k=1,L
   write(99,'(/a2,i2,a45/)')'第',k,'个通道的重建序列,每列为对应典型波形的重建序列'
   do j=1,n
      write(99,100)j,(Rccs(i,j,k),i=1,itp)
   enddo
enddo
close(99)
deallocate(raX,raEgvt,raTime,Rccs)
end

!※※※※※※※※※※※※※※※※※※※※※※※※※※※※※
subroutine sub_mssa_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


!----------------------------------
!MSSA的重建成份
subroutine sub_mssa_rc(n,m,L,raVector,raTime,Rccs)
implicit none
integer,intent(in)::n,m,L
real,dimension(L*m,L*m),intent(in)::raVector,raTime(L*m,n-m+1)
real,dimension(L*m,n,L),intent(out)::Rccs
integer::i,j,k,j2
real::con

do k=1,L; do i=1,L*m;;do j=1,n
  if((j >= 1) .and. (j <= m-1))then
    con=0.0
	do j2=1,j
       con=con+raVector((k-1)*m+j2,i)*raTime(i,j-j2+1)
    enddo
    Rccs(i,j,k)=con/float(j)
  elseif((j >= m) .and. (j <= n-m+1))then
    con=0.0
    do j2=1,m
       con=con+raVector((k-1)*m+j2,i)*raTime(i,j-j2+1)
    enddo
    Rccs(i,j,k)=con/float(m)
  elseif((j >= n-m+2) .and. (j <= n))then
    con=0.0
    do j2=j-n+m,m
       con=con+raVector((k-1)*m+j2,i)*raTime(i,j-j2+1)
    enddo
    Rccs(i,j,k)=con/float(n-j+1)
  endif
enddo;enddo;enddo

end


!=========================================================!
subroutine sub_standard2(m,n,x,ki)                        !
implicit none                                             !
integer,intent(in)::m,n,ki                                !
integer::i,j                                              !
real,dimension(m,n),intent(InOut)::x                      !
real::ave,xigma                                           !
real,intrinsic::sqrt                                      !
                                                          !
if(ki == -1)then                                          !
  do i=1,m; do j=1,n                                      !
     x(i,j)=x(i,j)                                        !
  enddo;enddo                                             !
elseif(ki == 0)then                                       !
  do i=1,m                                                !
     ave=0.0                                              !
     do j=1,n                                             !
        ave=ave+x(i,j)/float(n)                           !
     enddo                                                !
     do j=1,n                                             !
	    x(i,j)=x(i,j)-ave                                 !
	 enddo                                                !
  enddo                                                   !
elseif(ki == 1)then                                       !
  do i=1,m                                                !
     ave=0.0; xigma=0.0                                   !
     do j=1,n                                             !
        ave=ave+x(i,j)/float(n)                           !
     enddo                                                !
     do j=1,n                                             !
        xigma=xigma+(x(i,j)-ave)**2/float(n)              !
     enddo                                                !
     xigma=sqrt(xigma)                                    !
     do j=1,n                                             !
        x(i,j)=(x(i,j)-ave)/xigma                         !
     enddo                                                !
  enddo                                                   !
else                                                      !
  write(*,*)'ki必须取-1,0或1.'                            !
endif                                                     !
end                                                       !
!=========================================================!