!------------------- 最大熵谱估计 ----------------------------------------!
!要求输入的有:                                                            !
!        1.X(n):需作最大熵谱分析的长度为n的序列;                          !
!        2.K:最高试验阶数,一般为N/2左右;                                  !
!        3.m:谱分析中的窗口长度;                                          !
!输出的有:(存放在filename文件中)                                          !
!        1.FPE(K):                                                        !
!        2.k0:最佳试验阶数;                                               !
!        3.各个周期对应的谱值;                                            !
!注意:                                                                    !
!        1.最大熵谱分析尤其适用于较短的时间序列;                          !
!        2.本程序调用了smaths.lib中的函数;                                !
!------------------------------------- 程正泉 2000.5 ---------------------!

!=========================================================================!
subroutine sub_mesa(n,K,m,X,filename)                                     !
implicit none                                                             !
integer,intent(in)::n,K,m                                                 !
real,dimension(n),intent(in)::X                                           !
real,allocatable,dimension(:)::A(:,:),Var,Rp,FPE,G                        !
integer::k0,i,j,ismin                                                     !
character*50 filename                                                     !
open(90,file=filename)                                                    !
                                                                          !
!----- 1.计算各阶试验模型 ----------------                                !
allocate(A(K,K),Var(0:K),Rp(0:K))                                         !
call sub_burg(n,K,X,A,Var,Rp)                                             !
                                                                          !
!----- 2.利用FPE准则确定最佳模型阶数k0 ---                                !
allocate(FPE(K))                                                          !
do i=1,K                                                                  !
   FPE(i)=Var(i)*(real(n)+real(i))/(real(n)-real(i))                      !
enddo                                                                     !
k0=ismin(K,FPE,K)                                                         !
                                                                          !
!----- 3.以k0阶的模型参数计算最大熵谱 ----                                !
allocate(G(0:m))                                                          !
call compute_spectrum(m,K,k0,A,Var,G)                                     !
                                                                          !
!----- 4.输出 ----------------------------                                !
write(90,*)'最大熵谱分析的输出结果'                                       !
write(90,*)                                                               !
do i=1,K                                                                  !
   write(90,'(1x,a4,i3,a3,f12.4)')'FPE(',i,')==',FPE(i)                   !
enddo                                                                     !
write(90,*)                                                               !
write(90,'(a16,i4)')'最佳模型参数为:',k0                                  !
write(90,*)                                                               !
do i=0,m                                                                  !
   write(90,'(a4,i4,4x,a4,f7.1,8x,f12.4)')'i==',i,'T==',2.0*m/i,G(i)      !
enddo                                                                     !
deallocate(A,Var,Rp,FPE,G)                                                !
end                                                                       !
!=========================================================================!


!=========================================================================!
subroutine sub_burg(n,K,X,A,Var,Rp)
! K:递推阶数
implicit none
integer,intent(in)::n,K
real,dimension(n),intent(in)::X
real,dimension(K,K),intent(out)::A,Var(0:K),Rp(0:K)
integer::i1,i2,i3
real::temp1,temp2,temp3,temp4

do i1=0,K
   if(i1 == 0)then
     Var(i1)=0.0
     do i2=1,n
        Var(i1)=Var(i1)+X(i2)**2/real(n)
     enddo
     Rp(i1)=Var(i1)
   elseif(i1 == 1)then
     temp1=0.0; temp2=0.0; temp3=0.0
	 do i2=2,n
	    temp1=temp1+X(i2)*X(i2-1)
		temp2=temp2+X(i2)*X(i2)
		temp3=temp3+X(i2-1)*X(i2-1)
	 enddo
	 A(i1,i1)=2*temp1/(temp2+temp3)
	 Rp(i1)=A(i1,i1)*Var(i1-1)
	 Var(i1)=(1-A(i1,i1)**2)*Var(i1-1)
   else
     temp3=0.0; temp4=0.0
     do i2=i1+1,n
        temp1=0.0; temp2=0.0
	    do i3=1,i1-1
	       temp1=temp1+A(i1-1,i3)*X(i2-i3)
		   temp2=temp2+A(i1-1,i3)*X(i2-i1+i3)
	    enddo
        temp3=temp3+(X(i2)-temp1)*(X(i2-i1)-temp2)
		temp4=temp4+(X(i2)-temp1)**2+(X(i2-i1)-temp2)**2
	 enddo
	 A(i1,i1)=2*temp3/temp4
	 do i2=1,i1-1
	    A(i1,i2)=A(i1-1,i2)-A(i1,i1)*A(i1-1,i1-i2)
	 enddo
	 Var(i1)=(1-A(i1,i1)**2)*Var(i1-1)
	 temp1=0.0
	 do i2=1,i1-1
	    temp1=temp1+A(i1-1,i2)*Rp(i1-i2)
	 enddo
	 Rp(i1)=temp1+A(i1,i1)*Var(i1-1)
   endif
enddo
end
!=========================================================================!


!=========================================================================!
subroutine compute_spectrum(m,K,k0,A,Var,G)
implicit none
integer,intent(in)::m,K,k0
real,parameter::pi=3.14149
real,dimension(K,K),intent(in)::A,Var(0:K)
real,dimension(0:m),intent(out)::G
real::tmp1,tmp2
integer::i,j
do i=0,m
   tmp1=0.0; tmp2=0.0
   do j=1,k0
      tmp1=tmp1+A(k0,j)*cos(pi*i*j/real(m))
	  tmp2=tmp2+A(k0,j)*sin(pi*i*j/real(m))
   enddo
   G(i)=Var(k0)/((1-tmp1)**2+tmp2**2)
enddo
end
!=========================================================================!