spectral analysis

From: Hui Wan <wan_at_nyahnyahspammersnyahnyah>
Date: Mon, 21 Nov 2005 09:08:36 +0100
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

This archive was generated by hypermail 2.2.0 : Mon Nov 21 2005 - 10:28:39 MST