- 积分
- 7263
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-12-20
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 nickbsb 于 2012-4-19 13:00 编辑
c-----*----------------------------------------------------6---------7--
c Subroutine for continuous spectrum analysis of an one-dimensional series
c x(i) (i=1,...,n).
c Input parameters and arrays: n,m, x(n),
c n: number of data
c m: biggest lag time length. Generally, m is between n/3 and n/10.
c x(n): raw series
c Output variables:
c ol(0:m): frequency array
c tl(0:m): periodic array
c sl(0:m): power spectrum. Its sum from 0 to m is 1.
c st95(0:m): 95% confidence upper limit of red or white noise spectrum
c strw(0:m): spectrum density of red or white noise
c By Dr. Li Jianping, May 12,1998.
PROGRAM Spectrum
c 资料序列长度，时长间隔 --m is between n/3 and n/10.
parameter(n=3630,m=1210)
real x(n),sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
open(10,file='e:\aam\zonalaam.grd',form='binary')
do i=1,n
read(10) x(i)
enddo
close(10)
call cspectrum(n,m,x,ol,tl,sl,st95,strw)
open(20,file='e:\aamps\ps.grd',form='binary')
! open(20,file='e:\aamps\ps.txt')
ccccc
ccccc 注意此处为 do i=0,m
ccccc
do i=0,m
write(20) tl(i),sl(i),st95(i),strw(i)
c write(*) tl(i),sl(i),st95(i),strw(i)
enddo
!close(20)
stop '程序计算结束，感谢气象家园言深深版主的指导！'
end
subroutine cspectrum(n,m,x,ol,tl,sl,st95,strw)
dimension x(n)
dimension sl(0:m),ol(0:m),tl(0:m),st95(0:m),strw(0:m)
dimension r(0:m) ! Work array
pi=3.1415926
call autocorrelation(n,m,x,r)
do 10 l=0,m
sl(l)=r(0)
do i=1,m-1
cc=pi*i/float(m)
sl(l)=sl(l)+r(i)*(1.+cos(cc))*cos(l*cc)
enddo
10 continue
sl(0)=sl(0)/2.
sl(m)=sl(m)/2.
do 20 l=0,m
sl(l)=sl(l)/m
20 continue
do 30 l=1,m
ol(l)=pi*l/m
tl(l)=2.*m/l
30 continue
ol(0)=0
tl(0)=100.*tl(1)
call cnoisetest(n,m,r(1),sl,st95,strw)
return
end
c-----*----------------------------------------------------6---------7--
c This is a subroutine for calculating the lag autocorrelations
c r(m) of a series x(i).
c input: n,m,and x(n)
c n: number of the data x(n)
c x: the raw series
c m: maximum lagged length; m: maximum lag time
c output: r(0:m)
c r: lag autocorrelations from 0 to m.
c By Dr. LI Jianping, April 3, 1998.
subroutine autocorrelation(n,m,x,r)
dimension x(n),r(0:m)
call meanvar(n,x,ax,sx,vx)
r(0)=1.
do i=1,n
x(i)=x(i)-ax
enddo
do 10 i=1,m
sxy=0.
do j=1,n-i
sxy=sxy+x(j)*x(i+j)
enddo
sxy=sxy/float(n)
r(i)=sxy/vx
10 continue
return
end
c-----*----------------------------------------------------6---------7--
c Computing the mean ax, standard deviation sx
c and variance vx of a series x(i) (i=1,...,n).
c input: n and x(n)
c n: number of raw series
c x(n): raw series
c output: ax, sx and vx
c ax: the mean value of x(n)
c sx: the standard deviation of x(n)
c vx: the variance of x(n)
c By Dr. LI Jianping, May 6, 1998.
subroutine meanvar(n,x,ax,sx,vx)
dimension x(n)
ax=0.
vx=0.
sx=0.
do 10 i=1,n
ax=ax+x(i)
10 continue
ax=ax/float(n)
do 20 i=1,n
vx=vx+(x(i)-ax)**2
20 continue
vx=vx/float(n)
sx=sqrt(vx)
return
end
c-----*----------------------------------------------------6---------7--
c Noise test for continuous spectrum (95% confidence level).
c strw: spectrum density of red or white noise
c chi**2 distribution with degrees of freedom v=(2*n-m/2)/m is used.
c By Dr. Li Jianping, May 12,1998.
subroutine cnoisetest(n,m,r1,sl,st95,strw)
dimension sl(0:m),st95(0:m),strw(0:m)
dimension chi95(30),chi99(30)
call chi_table(chi95,chi98,chi99)
pi=3.1415926
as=0.
do 10 l=0,m
as=as+sl(l)
10 continue
as=as/float(m+1)
r2=r1*r1
v=(2.*n-m/2.)/float(m)
if(v.gt.30.)then !! invaliable from mathmatic table for v>30
write(*,*)'Unable continue... Please reinput m !'
stop
endif
iv=int(v)
c95=chi95(iv)+(v-iv)*(chi95(iv+1)-chi95(iv))
if(r1.gt.0.)then
c tname='red noise test'
do l=0,m
r3=1+r2-2*r1*cos(l*pi/m)
strw(l)=as*(1.-r2)/r3
st95(l)=strw(l)*c95/v
enddo
else
c tname='white noise test'
do l=0,m
strw(l)=as
st95(l)=as*c95/v
enddo
endif
return
end
c-----*----------------------------------------------------6---------7--
c chi-square distribution
c Chi-table with right tail probabilities
c P(X**2>=Xa**2)=a
c where a (alpha) significance level and (1-a)*100% is confindence level.
c chi95: chi-square distribution at 95% confidence level.
c chi98: chi-square distribution at 98% confidence level.
c chi99: chi-square distribution at 99% confidence level.
c By Dr. Li Jianping, May 12,1998.
subroutine chi_table(chi95,chi98,chi99)
dimension chi95(30),chi98(30),chi99(30)
dimension c95(30),c98(30),c99(30) !work array
data c95/ 3.841, 5.991, 7.815, 9.488,11.070,
* 12.592,14.067,15.507,16.919,18.307,
* 19.675,21.026,22.362,23.685,24.996,
* 26.296,27.587,28.869,30.144,31.410,
* 32.671,33.924,35.172,36.415,37.652,
* 38.885,40.113,41.337,42.557,43.773/
data c98/ 5.412, 7.824, 9.837,11.668,13.388,
* 15.033,16.622,18.168,19.679,21.161,
* 22.618,24.054,25.472,26.873,28.259,
* 29.633,30.995,32.346,33.687,35.020,
* 36.343,37.659,38.968,40.270,41.566,
* 42.856,44.140,45.419,46.693,47.662/
data c99/ 6.635, 9.210,11.345,13.277,15.086,
* 16.812,18.475,20.090,21.666,23.209,
* 24.725,26.217,27.688,29.141,30.578,
* 32.000,33.409,34.805,36.191,37.566,
* 38.932,40.289,41.638,42.980,44.314,
* 45.642,46.963,48.278,49.588,50.892/
do i=1,30
chi95(i)=c95(i)
chi98(i)=c98(i)
chi99(i)=c99(i)
enddo
return
end
要输出这些物理量对吧
ol(0:m): frequency array
tl(0:m): periodic array
sl(0:m): power spectrum. Its sum from 0 to m is 1.
st95(0:m): 95% confidence upper limit of red or white noise spectrum
strw(0:m): spectrum density of red or white noise
但画图时怎么很多是直线呢
|
-
-
30-60.grd
14.18 KB, 下载次数: 18, 下载积分: 金钱 -5
输入资料
|
窥视卡
雷达卡
发表于 2012-4-19 12:58:37
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
