Hello,
I'm trying to estimate the power spectrum of a time series using the
auto-correlation coefficients.
The auto-correlation function is symmetric about the origin, thus what
we need is the cosine transform of the auto-correlation values for lags
of 0 to "maxlag" time units. Firstly I "mirrored" the auto-correlation
as follows:
;-------------------------------
; 1.1
pc = File->pc ; the
time series
acr_srs_tmp = esacr( pc,maxlag ) ; auto-correlation
acr_srs = new( maxlag*2+1, float )
acr_srs(0:maxlag-1) = acr_srs_tmp(1::-1)
acr_srs(maxlag:) = acr_srs_tmp(:)
Then perform the spectral analysis:
;-------------------------------
; 1.2 Power spectrum
iopt = 0 ; remove series mean
jave = 0 ; no smoothing; return periodogram estimates
pct = 0 ; no tapering
return_value = specx_anal ( acr_srs(ipc,:),iopt,jave,pct )
raw_spec = return_value@spcx
frequency = return_value@frq
period = 1/frequency
pbar = 2*return_value@xvari ; the mean
spectral power
weights = (/0.25,0.5,0.25/)
smthed_spec = wgt_runave( raw_spec,weights,1 ) ; Hanning smoothing
The
resulting spectrum had the same shape as that obtained by applying
"specx_anal" directly to the original time series. (Please refer to the
attachment). Then I calculated the Markov spectrum and the lower and
upper confidence curves following WMO Technical Note No. 79 ("Climate
Change") as follows:
;----------------------------------
; 1.3 spectrum of the red noise
acr_pc = esacr( pc,1 )
acr1 = acr_pc(1) ; lag-1 auto-correlation
red_spec1 =
pbar*(1.-acr1*acr1)/(1-2*acr1*cos(2*pi*frequency)+acr1*acr1)
;----------------------------------
; 1.4 confidence bounds
b_lower =
0.05/(maxlag+1.0) ; 5% confidence level, WMO TN79, p41
b_upper = 1.- b_lower ; 95% confidence level
df =
(2.*ntime-0.5*maxlag)/maxlag ; WMO TN79, p40
xLower = chiinv
(b_lower,df)/df ; lower bound
xUpper = chiinv (b_upper,df)/df ; upper bound
data = new( (/4,nwave/),float )
data(0,:) = smthed_spec(:) ; smoothed
spectrum
data(1,:) = red_spec1(:)*xLower ; lower bound
data(2,:) = red_spec1(:) ; red noise
data(3,:) = red_spec1(:)*xUpper ; upper bound
;------------------------------------
; 1.5 plot settings ...
...
xsp = period
plot_spec = gsn_csm_xy(wks,xsp,data,ResL)
However it seemed that the
reference curves (red lines in the plots) were not consistent with the
spectrum of my data. Did I misunderstand anything or misuse any
function?
Thanks.
Hui
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 21 2005 - 01:08:36 MST