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

[Sri nandini]

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)

-------------- next part --------------
A non-text attachment was scrubbed...
Name: PSD_wNA_Monthly_Hist.png
Type: image/png
Size: 60311 bytes
Desc: not available
URL: <https://mailman.ucar.edu/pipermail/ncl-talk/attachments/20210302/f8a45f78/attachment.png>