Re: spectral analysis

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