<div dir="ltr">In my opinion, the spectrum reveals a classic red-noise spectrum. Only low frequency variance is present. <br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Mar 2, 2021 at 5:53 AM Sri nandini via ncl-talk <<a href="mailto:ncl-talk@mailman.ucar.edu">ncl-talk@mailman.ucar.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hello dear ncl-users,<br>
<br>
I have a question regarding computing frequency power spectra density <br>
for my selected regions, e.g. Northwestern Atlantic SSH anomalies.<br>
<br>
1) I calculated the frequency spectra of SSH anomalies by (1) extracting <br>
region and subtracting the ensemble mean, (2) temporal removal of mean <br>
and detrending was applied, as well as smoothing (periodogram estimates <br>
using modified Daniel) and percent tapering of the anomalies, (3) the <br>
function for spectral analysis returns the degrees of freedom -spectra <br>
for each ensemble, so i used a loop over all ensemble members to get the <br>
ensemble mean spectra (4) the variables returned are variance/(unit <br>
frequency interval) (y-axis) and frequency (cycles/months) for <br>
x-axis==see attached plot for power Spectra for historical Northwestern <br>
Atlantic SSH anomalies. I don't have a reference for this figure to <br>
compare with but it doesn't look meaningful, but do you think this is <br>
correct?<br>
<br>
Would be grateful, below is my code<br>
<br>
Sri<br>
<br>
<br>
f = addfile ("Hist_monthly_corrected_pi_<a href="http://trend.nc" target="_blank">trend.nc</a>", "r")<br>
hist_anom = f->hist_trend(1632:1871,:,{-30:60},{0:-40})<br>
hist_anom_NY1=dim_avg_n_Wrap(hist_anom,(/2,3/))<br>
printVarSummary(hist_anom_NY1)<br>
hist_anom_NY1=hist_anom_NY1*100<br>
dimx = dimsizes(hist_anom_NY1)<br>
ntim = dimx(0) ; 240<br>
nens = dimx(1) ;100<br>
hist_ensavg=dim_avg_n(hist_anom_NY1,1)<br>
printVarSummary(hist_ensavg)<br>
hist_anom1 = hist_anom_NY1<br>
do n=0,nens-1<br>
hist_anom1(:,n) = hist_anom_NY1(:,n) - hist_ensavg<br>
end do<br>
<br>
nLen = ntim/2<br>
spec = new ( nLen, typeof(hist_anom1))<br>
spec = 0.0<br>
z1 = 0.0<br>
xtsVar= 0.0<br>
printVarSummary(nLen)<br>
printVarSummary(spec)<br>
printVarSummary(hist_anom1)<br>
;-------------------------------------;<br>
;set function agurments<br>
;-------------------------------------;<br>
<br>
iopt = 1 ; detrending opt: 0=>remove mean <br>
1=>remove mean and detrend<br>
jave = 7 ; Average 7 periodogram estimates <br>
using modified Daniel<br>
pct = 0.1 ; percent tapered: (0.0 <= pct <= <br>
1.0) 0.10 common.<br>
;-------------------------------------;<br>
; calculate mean spectrum:specx_anal returns the degrees of freedom as <br>
a scalar.<br>
; calculate confidence interval [here 5 and 95%] and lag1 auto cor.<br>
; return 4 curves to be plotted<br>
; loop over every ensemble member<br>
;-------------------------------------;<br>
<br>
do ns=0,nens-1 ; spectra for each ensemble member<br>
sdof = specx_anal(hist_anom1(:,ns),iopt,jave,pct)<br>
spec = spec + sdof@spcx ; sum spectra<br>
xtsVar = xtsVar + sdof@xvari ; sum ensemble member variances<br>
; sum z-traansformed lag-1 correlations<br>
z1 = z1 + 0.5*log((1+sdof@xlag1)/(1-sdof@xlag1)) ; Fischer <br>
z-transform<br>
printVarSummary(z1)<br>
end do<br>
<br>
spec = spec/nens ; average spectra ==> ensemble <br>
mean spectra<br>
z1 = z1/nens ; ensemble mean z=transformed lag1<br>
xvari = xtsVar/nens ; mean variance of ensemble members<br>
printVarSummary(spec)<br>
printVarSummary(z1)<br>
printVarSummary(xvari)<br>
printVarSummary(sdof)<br>
;-----------------------------------------------------------------------------;<br>
; create a new variable<br>
;-----------------------------------------------------------------------------;<br>
<br>
ss = (/ sdof /)<br>
ss@dof = 2*nens ; # degrees freedom<br>
ss@bw = sdof@bw ; representing the spectral <br>
band width.<br>
ss@spcx = spec ;One-dimensional array of <br>
length N/2.units are variance/(unit frequency interval).<br>
ss@xvari = xvari ; representing the variance <br>
of the x series on input.<br>
ss@frq = sdof@frq ; A one-dimensional array of <br>
length N/2 representing frequency (cycles/time).<br>
ss@xlag1 = (exp(2*z1)-1)/(exp(2*z1)+1) ; match specx_anal<br>
printVarSummary(ss)<br>
splt = specx_ci(ss, 0.05, 0.95) ; calc confidence interval<br>
printVarSummary(splt)<br>
;-----------------------------------------------------------------------------;<br>
; plotting parameters<br>
;-----------------------------------------------------------------------------;<br>
wks = gsn_open_wks("png","PSD_wNA_Monthly_Hist")<br>
<br>
res = True ; plot <br>
mods desired<br>
res@gsnMaximize = True ; blow up plot<br>
res@gsnDraw = False<br>
res@gsnFrame = False<br>
<br>
res@vpXF = 0.1<br>
res@vpYF = 0.7<br>
res@vpWidthF = 0.7<br>
res@vpHeightF = 0.25<br>
<br>
res@trXMinF = 0<br>
res@trXMaxF = .50<br>
res@trYMinF = 20<br>
res@trYMaxF = 0<br>
<br>
res@gsnStringFontHeightF = 0.018<br>
res@gsnRightString = " "<br>
res@tiXAxisString = "Frequency (cycles/month)" ;the fraction (or <br>
use Period(months))<br>
res@tiXAxisFontHeightF = 0.018<br>
<br>
res@tiYAxisString = "PSD (cm^2/month^-1)" ;"Variance/frq_interval<br>
res@tiYAxisFontHeightF = 0.018<br>
<br>
res@tmXBLabelFontHeightF = 0.016<br>
res@tmYLLabelFontHeightF = 0.016<br>
<br>
;-------------------------------------;<br>
; 4-2. Timeseries plots<br>
;-------------------------------------;<br>
<br>
colors = (/2, 71, 48/)<br>
colors_em = (/15, 84, 60/)<br>
<br>
res@gsnLeftString = "a) Historical: Northwest Atlantic, SSHa <br>
Frequency Spectra"<br>
res@xyLineThicknesses = (/2.,1.,1.,1./) ; Define line thicknesses<br>
res@xyDashPatterns = (/0,0,1,1/) ; Dash patterns<br>
<br>
res@xyLineColor = colors(0)<br>
plot = gsn_csm_xy(wks,ss@frq, ss@spcx,res)<br>
<br>
_______________________________________________<br>
ncl-talk mailing list<br>
<a href="mailto:ncl-talk@mailman.ucar.edu" target="_blank">ncl-talk@mailman.ucar.edu</a><br>
List instructions, subscriber options, unsubscribe:<br>
<a href="https://mailman.ucar.edu/mailman/listinfo/ncl-talk" rel="noreferrer" target="_blank">https://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote></div>