[ncl-talk] compute frequency power density power spectra for selected regions

Dennis Shea shea at ucar.edu
Wed Mar 3 07:03:53 MST 2021


In my opinion, the spectrum reveals a classic red-noise spectrum. Only low
frequency variance is present.

On Tue, Mar 2, 2021 at 5:53 AM Sri nandini via ncl-talk <
ncl-talk at mailman.ucar.edu> wrote:

> Hello dear ncl-users,
>
> I have a question regarding computing frequency power spectra density
> for my selected regions, e.g. Northwestern Atlantic SSH anomalies.
>
> 1) I calculated the frequency spectra of SSH anomalies by (1) extracting
> region and subtracting the ensemble mean, (2) temporal removal of mean
> and detrending was applied, as well as smoothing (periodogram estimates
> using modified Daniel) and percent tapering of the anomalies, (3) the
> function for spectral analysis returns the degrees of freedom -spectra
> for each ensemble, so i used a loop over all ensemble members to get the
> ensemble mean spectra (4) the variables returned are variance/(unit
> frequency interval) (y-axis) and frequency (cycles/months) for
> x-axis==see attached plot for power Spectra for historical Northwestern
> Atlantic SSH anomalies. I don't have a reference for this figure to
> compare with but it doesn't look meaningful, but do you think this is
> correct?
>
> Would be grateful, below is my code
>
> Sri
>
>
>   f     = addfile ("Hist_monthly_corrected_pi_trend.nc", "r")
>    hist_anom    = f->hist_trend(1632:1871,:,{-30:60},{0:-40})
>    hist_anom_NY1=dim_avg_n_Wrap(hist_anom,(/2,3/))
>    printVarSummary(hist_anom_NY1)
>    hist_anom_NY1=hist_anom_NY1*100
>     dimx = dimsizes(hist_anom_NY1)
>     ntim = dimx(0)          ; 240
>     nens = dimx(1)         ;100
>    hist_ensavg=dim_avg_n(hist_anom_NY1,1)
>    printVarSummary(hist_ensavg)
>    hist_anom1 = hist_anom_NY1
>    do n=0,nens-1
>    hist_anom1(:,n) = hist_anom_NY1(:,n) - hist_ensavg
>    end do
>
>    nLen = ntim/2
>    spec  = new ( nLen, typeof(hist_anom1))
>    spec  = 0.0
>    z1    = 0.0
>    xtsVar= 0.0
>    printVarSummary(nLen)
>    printVarSummary(spec)
>    printVarSummary(hist_anom1)
>   ;-------------------------------------;
>   ;set function agurments
>   ;-------------------------------------;
>
>    iopt  = 1                           ; detrending opt: 0=>remove mean
> 1=>remove mean and detrend
>    jave  = 7                           ; Average 7 periodogram estimates
> using modified Daniel
>    pct   = 0.1                         ; percent tapered: (0.0 <= pct <=
> 1.0) 0.10 common.
>   ;-------------------------------------;
>   ; calculate mean spectrum:specx_anal returns the degrees of freedom as
> a scalar.
>   ; calculate confidence interval [here 5 and 95%] and lag1 auto cor.
>   ; return 4 curves to be plotted
>   ; 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
>     printVarSummary(z1)
>     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
>    printVarSummary(spec)
>    printVarSummary(z1)
>    printVarSummary(xvari)
>    printVarSummary(sdof)
>
>   ;-----------------------------------------------------------------------------;
>   ; create a new variable
>
>   ;-----------------------------------------------------------------------------;
>
>    ss       = (/ sdof /)
>    ss at dof   = 2*nens                       ; # degrees freedom
>    ss at bw    = sdof at bw                      ; representing the spectral
> band width.
>    ss at spcx  = spec                         ;One-dimensional array of
> length N/2.units are variance/(unit frequency interval).
>    ss at xvari = xvari                        ; representing the variance
> of the x series on input.
>    ss at frq     = sdof at frq                   ; A one-dimensional array of
> length N/2 representing frequency (cycles/time).
>    ss at xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1)  ; match specx_anal
>    printVarSummary(ss)
>    splt = specx_ci(ss, 0.05, 0.95)               ; calc confidence interval
>    printVarSummary(splt)
>
> ;-----------------------------------------------------------------------------;
> ; plotting parameters
>
> ;-----------------------------------------------------------------------------;
>    wks  = gsn_open_wks("png","PSD_wNA_Monthly_Hist")
>
>    res  = True                                                  ; plot
> mods desired
>    res at gsnMaximize   = True                        ; blow up plot
>    res at gsnDraw = False
>    res at gsnFrame = False
>
>    res at vpXF = 0.1
>    res at vpYF = 0.7
>    res at vpWidthF = 0.7
>    res at vpHeightF = 0.25
>
>    res at trXMinF = 0
>    res at trXMaxF = .50
>    res at trYMinF = 20
>    res at trYMaxF = 0
>
>    res at gsnStringFontHeightF = 0.018
>    res at gsnRightString = " "
>    res at tiXAxisString = "Frequency (cycles/month)"      ;the fraction (or
> use Period(months))
>    res at tiXAxisFontHeightF = 0.018
>
>    res at tiYAxisString = "PSD (cm^2/month^-1)" ;"Variance/frq_interval
>    res at tiYAxisFontHeightF = 0.018
>
>    res at tmXBLabelFontHeightF = 0.016
>    res at tmYLLabelFontHeightF = 0.016
>
>   ;-------------------------------------;
>   ; 4-2. Timeseries plots
>   ;-------------------------------------;
>
>    colors = (/2, 71, 48/)
>    colors_em = (/15, 84, 60/)
>
>    res at gsnLeftString  = "a) Historical: Northwest Atlantic, SSHa
> Frequency Spectra"
>    res at xyLineThicknesses   = (/2.,1.,1.,1./)       ; Define line
> thicknesses
>    res at xyDashPatterns      = (/0,0,1,1/)           ; Dash patterns
>
> res at xyLineColor = colors(0)
> plot = gsn_csm_xy(wks,ss at frq, ss at spcx,res)
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at mailman.ucar.edu
> List instructions, subscriber options, unsubscribe:
> https://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210303/2f282a49/attachment.html>


More information about the ncl-talk mailing list