Re: normalization in specx_anal

Hi Dennis,
It seems that multiplying the normalized spectrum by xvari/xvaro is not
the right way to obtain the raw non-normalized periodogram, since the
normalization factor is x_var/total_area. The built-in function
specx_anal gives me back x_var=xvaro but not the non-normalized area. So
the only way to obtain the raw non-normalized spectrum is to compute it
step by step.
Thank you for the script.

Franco

Il giorno ven, 13/03/2009 alle 11.32 -0600, Dennis Shea ha scritto:
> 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
> >
> Documento in testo semplice allegato (TST_spec_anal.ncl_v0)
> ;************************************************************
> ; 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
>
>

-- 
____________________________________________________
Eng. Franco Catalano
Ph.D. Student
D.I.T.S.
Department of Hydraulics, Transportation and Roads.
Via Eudossiana 18, 00184 Rome 
Sapienza University of Rome.
tel: +390644585218
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk