Log-scale power spectral density estimates


Univariate power spectral density on a logarithmic scale can be performed by the LPSD algorithm, which is an application of Welch's averaged, modified periodogram method, spectral density estimates are not evaluated at freqencies which are linear multiples of the minimum frequency resolution 1/T where T is the window lenght, but on a logarithmic scale. The algorithm takes care of calculating the frequencies at which to evaluate the spectral estimate, aiming at minimizing the uncertainty in the estimate itself, and to recalculate a suitable window length for each frequency bin.

ltpda_lpsd estimates the power spectral density of time-series signals, included in the input AOs. Data are windowed prior to the estimation of the spectrum, by multiplying it with a spectral window object, and can be detrended by polinomial of time in order to reduce the impact of the border discontinuities. Detrending is performed on each individual window. The user can choose the quantity being given in output among ASD (amplitude spectral density), PSD (power spectral density), AS (amplitude spectrum), and PS (power spectrum).

Syntaxis

	b = ltpda_lpsd(a1,a2,a3,...,pl)
  

a1, a2, a3, ... are AO(s) containing the input time series to be evaluated. b includes the output object(s). The parameter list pl includes the following parameters:

The length of the window is recalculated for each freqency bin, so the 'Win' parameter is used only to provide the key features of the window, i.e. the name and, for Keiser windows, the PSLL.
  If the user doesn't specify the value of a given parameter, the default value is used.
  

Examples

1. Evaluation of the ASD of a time-series represented by a low frequency sinewave signal, superimposed to white noise. Comparison of the effect of using standard Pwelch and LPSD on the estimate of the white noise level and on resolving the signal.

	x1 = ao(plist('waveform','sine wave','f',0.1,'A',1,'nsecs',1000,'fs',10)); 
	x2 = ao(plist('waveform','noise','type','normal','nsecs',1000,'fs',10));
    x = x1 + x2;
    pl = plist('scale','ASD','order',-1,'win',specwin(plist('name','BH92')));
	y1 = ltpda_pwelch(x, pl);
    y2 = ltpda_lpsd(x,pl);
	iplot(y1, y2) 
  




©LTP Team