From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>

Date: Mon, 21 Nov 2005 11:22:54 -0700 (MST)

Date: Mon, 21 Nov 2005 11:22:54 -0700 (MST)

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:33 MST
*