Re: mean power spectrum computation

Hi Yan,

shouldn't xlag1 the attribute of df?

  df@xlag1 = r1

Maybe this helps,
Karin

Am 20.04.14 23:38, schrieb yjin7@masonlive.gmu.edu:
>
> 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

-- 
Dipl. Geophys. Karin Meier-Fleischer
Application Support, Visualization
Deutsches Klimarechenzentrum GmbH      E-Mail:   meier-fleischer@dkrz.de
Bundesstrasse 45a                      Internet: http://www.dkrz.de/
20146 Hamburg                          Phone:    +49 (0)40 460094 126
Germany                                Fax:      +49 (0)40 460094 270
Geschäftsführer: Prof. Dr. Thomas Ludwig
Sitz der Gesellschaft: Hamburg
Amtsgericht Hamburg HRB 39784

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Apr 21 04:47:42 2014