<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body>
<p>Hello dear ncl-users,</p>
<p>Would anyone know how to remove this red-noise spectrum so i can
show the peaks of high variance for my plot?</p>
<p>Best</p>
<p>sri<br>
</p>
<div class="moz-cite-prefix">On 03.03.21 15:03, Dennis Shea wrote:<br>
</div>
<blockquote type="cite"
cite="mid:CAOF1d_6-bkQbJizxrfFRsLsfV+WpCah0=GRaH7QBCpy3LB4VEw@mail.gmail.com">
<meta http-equiv="content-type" content="text/html; charset=UTF-8">
<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"
moz-do-not-send="true">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"
moz-do-not-send="true">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"
moz-do-not-send="true">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" moz-do-not-send="true">https://mailman.ucar.edu/mailman/listinfo/ncl-talk</a></blockquote>
</div>
</blockquote>
</body>
</html>