The short answer is "yes"
------
Attached is a very simple script that
calculates the power spectrum using the builtin
function and then a step-by-step example.
The raw non-normalized periodogram is line 69.
Then a weighted running average is passed over the
periodogram estimates [ line 79 ]. This gives the power spectrum.
The non-normalized area is at lines 85-86.
The normalization is at line 91,
--- Franco Catalano wrote: > Dear ncl developers and users, > I want to obtain the non-normalized spectrum from the function > specx_anal. > I guess I just have to multiply the normalized spectrum spcx for this > normalization factor: xvari/xvaro > Is it correct? Otherwise, what is the right way to obtain the > non-normalized spectrum from spcx? > Thank you for your suggestions. > > Franco >
;************************************************************
; purpose : Compare NCL function specx_anal explicit steps
;************************************************************
begin
; define a time series
x = (/ 1002, 1017, 1018, 1020, 1018, 1027, \
1028, 1030, 1012, 1012, 982, 1012, \
1001, 996, 995, 1011, 1027, 1025, \
1030, 1016, 996, 1006, 1002, 982 /)
N = dimsizes(x)
;************************************************
; set "specx_anal" function arguments
;************************************************
d = 1 ; detrending opt: 0=>remove mean 1=>remove mean and detrend
sm = 3 ; smoothing periodogram: should be at least 3 and odd
pct = 0.10 ; percent tapered: (0.0 <= pct <= 1.0) 0.10 common.
;************************************************
; calculate spectrum using the built-in function
;************************************************
spec = specx_anal(x,d,sm,pct)
print (spec)
;************************************************
; use spectral estimates to calculate the variance
;************************************************
df = 1.0/N
area_under = (spec_at_spcx(0) + spec_at_spcx(N/2 - 1))*(df/2)\ ; half width
+ sum(spec_at_spcx(1:(N/2-2)))*df ; full width
;***********************************************
; step by step to calculation according to
; the introdcution on the web site:
; http://www.ncl.ucar.edu/Document/Functions/Built-in/specx_anal.shtml
;***********************************************
;step 1: detrend the series (The mean is also removed)
;***********************************************
x_dtrend = dtrend(x,False)
;************************************************
;step 2: taper the series
;************************************************
x_taper = taper(x_dtrend,0.1,0)
;************************************************
;step 3: Calculate the variance of the series
; Modify variance to match specx_anal
;************************************************
x_var1 = variance(x)*(N-1.)/N ; raw series
x_var = variance(x_dtrend)*(N-1.)/N ; detrended series
;************************************************
;step 4: perform a forward FFT the series
;************************************************
cf = ezfftf(x_taper)
;************************************************
;step 5: square the real (rc) and imaginary (ic) coef from [d].
; This gives the periodogram
;************************************************
perio_dog = cf(0,:)^2 + cf(1,:)^2
print("======================")
print("Real coef Imag coef")
print(cf(0,:) +" "+cf(1,:)) ; real/imag
;************************************************
;step 6: perform a running average over the periodogram values using
; wrunave. This creates the spectrum.
;************************************************
spectrum = wgt_runave(perio_dog,(/0.5,1,0.5/),1)
;************************************************
;step 7: Calculate the unnormalized variance
; Remember the first and last bandwidths are only 0.5*df
;************************************************
total_area = (spectrum(0) + spectrum(N/2-1))*(df/2) + \
sum(spectrum(1:N/2-2))*df
;************************************************
;step 8: normalize the area under the curve to the input variance.
;************************************************
spec_est = spectrum*(x_var/total_area)
; normalized area
nml_area = (spec_est(0) + spec_est(N/2-1))*(df/2) + \
sum(spec_est(1:N/2-2))*df
print("**********************")
print("NCL specx_anal function results")
print("**********************")
print("the spectral band width = " + spec_at_bw)
print("the variance of the x series on input = "+ spec_at_xvari)
print("the variance of the x series after detrending = " + spec_at_xvaro)
print("df=(1/N)=frequency spacing = "+ df)
print("the area under the curve = " + area_under)
print("**********************")
print("spectral estimate at "+spec_at_frq + " is "+spec_at_spcx)
print(" ")
print("&&&&&&&&&&&&&&&&&&&&&&")
print("step by step results")
print("&&&&&&&&&&&&&&&&&&&&&&")
print(" ")
frq = fspan(df,0.5,N/2)
print("the variance of the x series on input = "+ x_var1)
print("the variance of the x series after detrending = " + x_var)
print("the area under the curve (before normalized) = " + total_area)
print("the area under the curve (normalized) = " + nml_area)
print("**********************")
print("spectral estimate at "+frq+" is "+spec_est)
end
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Mar 13 2009 - 11:32:21 MDT
This archive was generated by hypermail 2.2.0 : Mon Mar 16 2009 - 10:31:01 MDT