Re: normalization in specx_anal

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Mon, 16 Mar 2009 08:05:02 -0600

Also, the script I sent allows you to use
different smoothing than the simple smoothing
used by specx_anal.

Good luck

Franco Catalano wrote:
> 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
>>
>>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Mar 16 2009 - 08:05:02 MDT

This archive was generated by hypermail 2.2.0 : Mon Mar 16 2009 - 10:31:01 MDT