Re: normalization in specx_anal

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