功率谱分析

大果子 · 发表于 2016-11-24 11:39:38

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

各位大神，我用下面这个fortran程序做了功率谱分析，但是做出来的图看着不太对，谁能帮我看看这个程序是否有问题？或者应该用哪些算出来的数据做图呢？

!------------------ 功率谱分析 -------------------------------!
!要求输入的有:                                              !
! X(n):原始序列;                                        !
! m:窗口长度;                                           !
! Alfa:检验的信度;                                        !
! filename:文件名,存放输出结果;                         !
!输出的内容有:                                              !
! 1.原序列为白噪声或红噪声;                               !
! 2.原序列的实际谱和给定信度alfa下的置信谱上界;          !
! 3.原序列的各个谱所能达到的信度水平;                   !
!注意: 该程序用到了sstats.lib中的统计函数;                   !
!-------------------------------------------------------------!
program power_pectrum
use imsl
use numerical_libraries
parameter(n=153,m=33,alfa=0.05)
character*50 FileName
real x(n)
FileName='150E-power.dat'
!-----  1.读入数据 -------------
open(1,file='k:\qj\hl\uave1981.grd',form='binary')
read(1)x
!read(1,*)x
!read(1,*)x

close(1)
111 format(f10.5)
!--------调用功率谱子程序--------------
call sub_spectrum(n,m,X,Alfa,filename)
end

!--------功率谱子程序-----------------------
subroutine sub_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 sub_one_backxg(N,m,x,R)
implicit none
integer::N,m,Tao
real,dimension(N)::x,R(0:m)
real::r_temp                                                                                                                                  !
do Tao=0,m
  call sub_correlation(N-Tao,x(1:N-Tao),x(1+Tao:N),r_temp)
  R(tao)=r_temp
enddo
end
!---------------------------------------------------------
subroutine sub_correlation(n,x,y,r)
implicit none
integer,intent(in)::n
real,dimension(n),intent(in)::x,y
real,intent(out)::r
real,intrinsic::sqrt,float
real::ave1,ave2,Var1,Var2,tmp
integer::i
ave1=0.0; ave2=0.0; Var1=0.0; Var2=0.0
do i=1,n
ave1=ave1+x(i)/float(n)
ave2=ave2+y(i)/float(n)
enddo
do i=1,n
Var1=Var1+(x(i)-ave1)**2
Var2=Var2+(y(i)-ave2)**2
enddo
tmp=0.0
do i=1,n
tmp=tmp+(x(i)-ave1)*(y(i)-ave2)
enddo
r=tmp/sqrt(Var1*Var2)
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

大果子 · 发表于 2016-11-24 11:41:33

补充，图中序列1是假设谱估计的置信上界，系列2是实际谱通过的信度检验，是用这两组数据画图吗？还是要用实际估计谱？

贫道敬孔 · 发表于 2016-11-26 18:44:51

同样不懂，关注一个。。

果糖子丫丫 · 发表于 2016-11-27 21:27:05

还缺一组数据吧

果糖子丫丫 · 发表于 2016-11-27 21:29:22

我的是这种

zhaojunri1020 · 发表于 2016-11-27 22:41:59

谢谢楼主！！！！

maiyatang_xihui · 发表于 2017-1-6 14:23:50

果糖子丫丫发表于 2016-11-27 21:29
我的是这种

请问你这个图95%检验和红噪声检验怎么做的？能分享下代码吗？

wuli茄茄 · 发表于 2017-3-21 12:52:06

同问也想知道红噪声的问题

yanana · 发表于 2018-3-11 18:04:12

谢谢分享!

ZMZ · 发表于 2019-4-8 23:51:27

楼主的sub_correlation(n,x,y,r)的y是什么

		自动登录	找回密码
密码			立即注册

[混合编程] 功率谱分析

登录后查看更多精彩内容~

本帖被以下淘专辑推荐: