Hello
You left out some information.
[1] size of pc?
[2] what is maxlag? According to Chatfield [see esacr documentation
for reference] maxlag should not be larger than dimsizes(pc)/4.
Using this length means that the autocorrelation approach
will not yield the same information for higher freq.
[3] What is "ipc" subscript?
[4] by setting iopt/jave/pct=0, then spec_anal
[a] doing a straight fft on the series
[b] squaring the coef
[c] normalizing the area under the periodogram
so that the variance under the curve matches
that of the input series.
The shape/details of spectrum are dependent upon
maxlag. Also, if you used maxlag> dimsizes(pc)/4,
remember that the series is normalized by the overall
series variance to minimize bias.
see comments below
>Hello,
>I'm trying to estimate the power spectrum of a time series using the
>auto-correlation coefficients.
why are you using this approach?
Nothing wrong with it ... just curious.
>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(:)
acr_srs(maxlag:) = acr_srs_tmp ; more efficient to not use (:)
> 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_at_spcx
> frequency = return_value_at_frq
> period = 1/frequency
> pbar = 2*return_value_at_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:
This is basically the same approach used in "specx_ci"
located in "shea_util.ncl"
;----------------------------------
; 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?
Well, when I look at the black dots ... I'd
say your spectrum is not consistent with a
first order Markov process.
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Nov 21 2005 - 11:22:54 MST
This archive was generated by hypermail 2.2.0 : Mon Nov 21 2005 - 17:26:39 MST