# [ncl-talk] calculate the ensemble mean power spectrum

Dennis Shea shea at ucar.edu
Tue Mar 3 20:29:02 MST 2020

```This type of question was recently answered on ncl-talk.

To avoid repeated emails, I extracted the pertinent section.

This has not been tested on your data. There may be errors.
It is your responsibility to fix them.
+++++++++++++++++++++++++++++++++++++++++++++++++
hist_anom1=hist_anom(:,:,{-7},{160})     ;(ntim, ens) ; (240,100)
=====
nLen = ntim/2
spec  = new ( nLen, typeof(hist_anom1))
spec  = 0.0
z1    = 0.0
xtsVar= 0.0

iopt  = 0
jave  = 1                           ; periodigram for each segment
pct   = 0.1
; loop over  every ensemble
member
do ns=0,nens-1                 ; spectra for each ensemble member
sdof   = specx_anal(hist_anom1(:.ns),iopt,jave,pct)
spec   = spec + sdof at spcx      ; sum spectra
xtsVar = xtsVar + sdof at xvari   ; sum ensemble member variances
; sum z-traansformed
lag-1 correlations
z1     = z1   + 0.5*log((1+sdof at xlag1)/(1-sdof at xlag1)) ; Fischer
z-transform
end do

spec  = spec/nens                      ; average spectra ==> ensemble
mean spectra
z1    = z1/nens                           ; ensemble mean z=transformed
lag1
xvari = xtsVar/nens                    ; mean variance of ensemble members

; create a new variable
ss           = (/ sdof /)
ss at dof   = 2*nens                       ; # degrees freedom
ss at bw    = sdof at bw
ss at spcx  = spec
ss at xvari = xvari
ss at frq     = sdof at frq
ss at xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1)  ; match specx_anal

printVarSummary(ss)

;************************************************
; plotting parameters
;************************************************
pltDir  = "./"
pltName = "ens_mean"
pltPath = pltDir+pltName
pltType = "png"

wks  = gsn_open_wks(pltType,pltPath)

r  = True                                                  ; plot mods
desired
r at gsnMaximize   = True                        ; blow up plot
r at tiMainString  = "???"                          ; title
r at tiXAxisString = "Frequency (cycles/day)"    ; xaxis
r at tiYAxisString = "Variance"                  ; yaxis
splt = specx_ci(ss, 0.05, 0.95)                ; calc confidence interval
r at xyLineThicknesses   = (/2.,1.,1.,1./)   ; Define line thicknesses
r at xyDashPatterns      = (/0,0,1,1/)          ; Dash patterns
plot = gsn_csm_xy(wks,ss at frq, splt,r)

On Tue, Mar 3, 2020 at 7:08 AM Sri.durgesh Nandini-Weiss via ncl-talk <
ncl-talk at ucar.edu> wrote:

> Dear ncl users,
>
> I am trying to perform power spectral analysis of a time series of
> ensemble data (time*ens==>20years*100ensemble to compute over the year
> and each ensemble member and than calculate the ensemble mean of the
> spectra.
>
> The best way would be to use specx_anal() and write it as a loop but
> this is where i am failing.
>
> I would later on like to know how to add the Markov red noise and the
> upper and lower confidence levels?
>
> would someone advice me on my script below?
>
>
> begin
>
>     f     = addfile ("anom_hist_zo_slr.nc", "r")
>     hist_anom    = f->hist_anom
>     printVarSummary(hist_anom)
>
>     f1     = addfile ("anom_rcp45_zo_slr.nc", "r")
>     rcp45_anom   = f1->rcp45_anom
>     printVarSummary(rcp45_anom)
>
>     dimx = dimsizes(hist_anom)
>     ntim = dimx(0)           ; 240
>     nens = dimx(1)           ; 100
>     nlat = dimx(2)           ; 45
>     nlon = dimx(3)           ; 90
>     lat1=nlat
>     lon1=nlon
>
>     nmos = 12
>     nyrs   = ntim/nmos       ; 20
>     printVarSummary(nyrs)
>
>
> ;==================================================================
> ; print ("Data read in - start calculations")
> ; Region domain lat and lon, ensemble average
> ;==================================================================
>
>    hist_anom1=hist_anom(:,:,{-7},{160}) ;(ntim, ens)
>    printVarSummary(hist_anom1)
>    aveX=dim_avg_n_Wrap(hist_anom1,(/1/))
>    printVarSummary(aveX)     ; [time]
>
>    rcp45_anom1=rcp45_anom(:,:,{-7},{160})                ;(ntim, ens)
>    printVarSummary(rcp45_anom1)
>    aveY=dim_avg_n_Wrap(rcp45_anom1,(/1/))
>    printVarSummary(aveY)      ; [time]
>
>
>    d   = 0                      ; detrending opt: 0=>remove mean
> 1=>remove mean + detrend
>    sm  = 7                    ; smooth periodogram:: should be at least
> 3 and odd
>                                     ; Average 7 periodogram estimates
> using modified Daniell
>    pct = 0.10                 ; percent taper: (0.0 <= pct <= 1.0) 0.10
> common.
>
> ;==================================================================
>    ;i can either use ndtooned and carry out specx_anal on this or make a
> loop.
>
>    aveX = ndtooned(hist_anom1)
>    printVarSummary(aveX)
>
>    aveY = ndtooned(rcp45_anom1)
>    printVarSummary(aveY)
>
>    spec = specx_anal(aveX,d,sm,pct)
>    printVarSummary(spec)
>
>    spec1 = specx_anal(aveY,d,sm,pct)
>    printVarSummary(spec1)
>
>
> Thanx
>
> Sri
>
>
> --
> Dr. Sri Nandini-Weiß
> Research Scientist
>
> Universität Hamburg
> Center for Earth System Research and Sustainability (CEN)
> Cluster of Excellence 'Climate, Climatic Change, and Society' (CLICCS)
>
> Bundesstrasse 53, 20146 Hamburg
> Tel: +49 (0) 40 42838 7472
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20200303/4124ef42/attachment.html>
```