Hello all,
I have two time series data and want to do cross spectrum analysis. Can you help me to check whether I have got the right coherence probility? My codes are as follow,
fn = sec1+"_daily.dat.bak"
data = asciiread(fn,(/nline,2/),"float")
xn = data(:,1)
fn = sec2+"_daily.dat.bak"
yn = data(:,1)
wgtt = filwgts_lanczos (301,2,1/400,1/80,1.0)
XN = wgt_runave (xn,wgtt,0)
YN = wgt_runave (yn,wgtt,0)
delete([/xn,yn/])
xn = XN((nwgt-1)/2:(dimsizes(XN)-(nwgt-1)/2-1))
yn = YN((nwgt-1)/2:(dimsizes(YN)-(nwgt-1)/2-1))
spec = specxy_anal(xn,yn,0,5,0.1)
dataL=new((/2,dimsizes(spec@coher)/),typeof(spec@coher))
dataL(0,:)=sqrt(spec@coher)
dataL(1,:)=spec@coher_probability(1)
phase=spec@phase
phase@_FillValue= 1e20
phase=where(sqrt(spec@coher).le.spec@coher_probability(1),phase@_FillValue,phase)
plot=gsn_csm_xy2(wks,spec@frq,dataL,phase,resL,resR); create plot
Am I rigth with the line
phase=where(sqrt(spec@coher).le.spec@coher_probability(1),phase@_FillValue,phase), in order to get the coherence satisfy the 95% confidence level?
Thanks.
2014-02-25
Xueming Zhu 朱学明
Key Laboratory of Research on Marine Hazards Forecasting (LoMF), SOA
National Marine Environmental Forecasting Center (NMEFC)
No.8, Dahuisi Road, Haidian District, Beijing, 100081
People's Republic of China
Tel:+86-10-82481923
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Feb 24 14:48:54 2014
This archive was generated by hypermail 2.1.8 : Mon Mar 03 2014 - 14:26:18 MST