mean power spectrum computation

Hi everyone,

When I used Example 5 to in http://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml to calculate the mean power spectrum, there is still an error message.

In my code, I defined a matrix pw, which is the same as pc2dT:

-------------------------------------------------------------
ncl 217> printVarSummary(pc2dT)

Variable: pc2dT
Type: double
Total Size: 29520 bytes
            3690 values
Number of Dimensions: 2
Dimensions and sizes: [year | 30] x [day | 123]
Coordinates:
Number Of Attributes: 1
  _FillValue : -999000000

---------------------------------------------------------------

The difference between my code and Example is: I defined r1 and r1zsum before the loop, otherwise the loop cannot continue. Also I defined splt before using splt = specx_ci(df, 0.05, 0.95), because it always said splt is undefined.

Here is my code:
;===============================
; power spectra
;===============================
;x(nseg,ntim), nseg is nyrs here, ntim is nday here;
  ntim = nday
  nseg = nyrs
  pw = pc2dT

  d = 0 ; detrending opt: 0=>remove mean 1=>remove mean + detrend
  sm = 21 ; smooth: should be at least 3 and odd
  pct = 0.10 ; percent taper: (0.0 <= pct <= 1.0) 0.10 common.

;************************************************
  ; calculate mean spectrum spectrum and lag1 auto cor
  ;************************************************

  ; loop over each segment of length ntim
  r1=new(1,double)
  r1=0.0
  r1zsum=new(1,double)
  r1zsum=0.0

   spcavg = new ( ntim/2, typeof(pw))
   spcavg = 0.0

   do n=0,nseg-1
      dof = specx_anal(pw(n,:),d,sm,pct) ; current segment spc
      spcavg = spcavg + dof@spcx ; sum spc of each segment
      r1 = dof@xlag1 ; extract segment lag-1
      r1zsum = r1zsum + 0.5*(log((1.0+r1)/(1.0-r1))) ; sum the Fischer Z
   end do

   r1z = r1zsum/nseg ; average r1z
   r1 = (exp(2*r1z)-1)/(exp(2*r1z)+1) ; transform back
                                                ; this is the mean r1
   spcavg = spcavg/nseg ; average spectrum

 ;************************************************
  ; Assign mean spectrum to data object
  ;************************************************

   df = 2.0*nseg ; deg of freedom
                                                ; all segments
   df@spcx = spcavg ; assign the mean spc
   df@frq = dof@frq
   spcavg@xlag1 = r1 ; assign mean lag-1

  ;************************************************
  ; plotting
  ;************************************************
    wks = gsn_open_wks("ps","spec") ; Opens a ps file

    res = True
    res@tiMainString = "Mean Spectra: "+nseg+" segments, dof="+df ; title
    res@tiXAxisString = "Frequency (cycles/month)" ; xaxis
    res@tiYAxisString = "Variance" ; yaxis

    splt = new((/4,dimsizes(df@spcx)/),float)

    splt = specx_ci(df, 0.05, 0.95) ; confidence interval

    plot = gsn_csm_xy(wks, df@frq, splt ,res) ; create plot

;===================================

However, The code is still stuck in splt = specx_ci(df, 0.05, 0.95). The error message is

--------------------------

warning:Attempt to reference attribute (xlag1) which is undefined
fatal:Assignment type mismatch, right hand side can't be coerced to type of left hand side
fatal:["Execute.c":8128]:Execute: Error occurred at or near line 1348 in file $NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl

fatal:["Execute.c":8126]:Execute: Error occurred at or near line 216
----------------------------------------------------------

?Thank you very much.
Best regards
Yan

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Apr 20 15:38:58 2014