!------------------ 功率谱分析 -------------------------------!
!要求输入的有:                                                !
!    X(n):原始序列;                                           !
!    m:窗口长度;                                              !
!    Alfa:检验的信度;                                         !
!    filename:文件名,存放输出结果;                            !
!输出的内容有:                                                !
!    1.原序列为白噪声或红噪声;                                !
!    2.原序列的实际谱和给定信度alfa下的置信谱上界;            !
!    3.原序列的各个谱所能达到的信度水平;                      !
!注意:                                                        !
!    1.该程序用到了sstats.lib中的统计函数;                    !
!    2.本程序用到了sub_correlation.f90文件;                   !
!---------------------------------------- 程正泉 2000.5 ------!
subroutine sub_self_spectrum(n,m,X,Alfa,filename)
implicit none
real,intent(in)::Alfa
integer,intent(in)::m,n
real,dimension(n),intent(in)::X
integer::i,j
real,allocatable,dimension(:)::R,Gk1,Gk2
character*50 FileName

!-----  1.求样本自相关函数
allocate(R(0:m))
call sub_one_backxg(N,m,X,R)

!-----  2.求粗谱估计
allocate(Gk1(0:m))
call think_spectrum(m,R,Gk1)

!-----  3.计算加Tukey-Hanning窗函数的平滑谱
allocate(Gk2(0:m))
call smooth_spectrum(m,Gk1,Gk2)

!-----  4.功率谱的显著性检验并输出
call spectrum_test(Alfa,N,m,R,Gk2,FileName)
deallocate(R,Gk1,GK2)

end

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!粗谱估计
subroutine think_spectrum(m,R,Gk1)
!(R(0:m):样本自相关函数)//(Gk(0:m):粗谱)
implicit none
integer::m,Tao,k
real::temp,pi
real,dimension(0:m)::R,Gk1
pi=2.0*asin(1.0)
do k=0,m
   temp=0.0
   do Tao=1,m-1
      temp=temp+R(Tao)*cos(k*pi*Tao/float(m))
   enddo
   Gk1(k)=(R(0)+2*temp+R(m)*cos(k*pi))/float(m)
enddo
end

!-----------------------
!求加Tukey-Hanning窗函数的平滑谱
subroutine smooth_spectrum(m,Gk1,Gk2)
implicit none
integer::m,i
real,dimension(0:m)::Gk1,Gk2
Gk2(0)=(Gk1(0)  +Gk1(1))*0.5
Gk2(m)=(Gk1(m-1)+Gk1(m))*0.5
do i=1,m-1
   Gk2(i)=0.25*Gk1(i-1)+0.5*Gk1(i)+0.25*Gk1(i+1)
enddo
end

!-----------------------
!功率品的显著性检验
subroutine spectrum_test(Alfa,N,m,R,Gk2,FileName)
implicit none
integer::N,m,i,j
real::Alfa,DF,avGk2,tin,chiin,chidf
real::r_alfa,pi
real,dimension(0:m),intent(in)::R,Gk2
real,allocatable,dimension(:)::G0k,Gk_alfa
character*50 Name,FileName
pi=2.0*asin(1.0)
open(90,file=FileName)
!  1.确定假设总体谱的具体类型(白噪声序列或红噪声序列)
df=1.0*(N-2)
r_alfa=tin(1.0-Alfa/2.0,df)/sqrt(N-2+tin(1.0-Alfa/2.0,df)**2)
if(R(1) <= r_alfa)then
  Name='WHITE';write(*,*)'原序列为白噪声';write(90,'(a14)')'原序列为白噪声'
else
  write(*,*)'R(1)=',R(1),'R_alfa=',R_alfa
  Name='RED';write(*,*)'原序列为红噪声';write(90,'(a14)')'原序列为红噪声'
endif
!  2.估计原假设下的总体谱
allocate(G0k(0:m))
if(Name == 'WHITE')then
  do i=0,m
     G0k(i)=0.0
	 do j=0,m
	    G0k(i)=G0k(i)+Gk2(j)/float(m+1)
	 enddo
  enddo
elseif(Name == 'RED')then
  avGk2=0.0
  do i=0,m
     avGk2=avGk2+Gk2(i)/float(m+1)
  enddo
  do i=0,m
     G0k(i)=(1-R(1)**2)/(1+R(1)**2-2*R(1)*cos(i*pi/float(m)))
	 G0k(i)=G0k(i)*avGk2
  enddo
endif
!  3.给出检验依据
allocate(Gk_alfa(0:m))
df=(2*N-0.5*m)/float(m)
do i=0,m
   Gk_alfa(i)=G0k(i)*chiin(1.0-alfa,df)/df
enddo
write(90,'(a14,f4.2,a14)')'通过信度Alfa==',Alfa,'的显著周期为: '
do i=0,m
  if(Gk2(i) > Gk_alfa(i))then
    write(90,'(a4,i3,6x,a4,f7.1)')'k==',i,'T==',2.0*m/i
  else
  endif
enddo
write(90,'(/a4,a12,a24)')'m','实际估计谱','假设谱估计的置信上界'
do i=0,m
  write(90,'(i4,f10.4,f18.4)')i,Gk2(i),Gk_Alfa(i)
enddo
write(90,*)
write(90,*)'各个实际谱通过的信度检验'
do i=0,m
  write(90,'(a3,i3,a5,f6.1,a10,f6.4)')'m==',i,'T==',2.0*m/i,'alfa==',1-chidf(Gk2(i)*df/G0k(i),df)
enddo
deallocate(G0k,Gk_alfa)  
end