[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