Signal Processing Toolbox™

pburg - PSD using Burg method

Syntax

Pxx = pburg(x,p)
Pxx = pburg(x,p,nfft)
[Pxx,w] = pburg(...)
[Pxx,w] = pburg(x,p,w)
Pxx = pburg(x,p,nfft,fs)
Pxx = pburg(x,p,f,fs)
[Pxx,f] = pburg(x,p,nfft,fs)
[Pxx,f] = pburg(x,p,f,fs)
[Pxx,f] = pburg(x,p,nfft,fs,'range')
[Pxx,w] = pburg(x,p,nfft,'range')
pburg(...)

Description

Pxx = pburg(x,p) implements the Burg algorithm, a parametric spectral estimation method, and returns Pxx, an estimate of the power spectral density (PSD) of the vector x. The entries of x represent samples of a discrete-time signal, and p is the integer specifying the order of an autoregressive (AR) prediction model for the signal, used in estimating the PSD.

The power spectral density is calculated in units of power per radians per sample. Real-valued inputs produce full power one-sided (in frequency) PSDs (by default), while complex-valued inputs produce two-sided PSDs.

In general, the length of the FFT and the values of the input x determine the length of Pxx and the range of the corresponding normalized frequencies. For this syntax, the (default) FFT length is 256. The following table indicates the length of Pxx and the range of the corresponding normalized frequencies for this syntax.

PSD Vector Characteristics for an FFT Length of 256 (Default)

Pxx = pburg(x,p,nfft) uses the integer FFT length nfft to calculate the PSD vector Pxx.

[Pxx,w] = pburg(...) also returns w, a vector of normalized angular frequencies at which the two-sided PSD is estimated. Pxx and w have the same length. The units for w are rad/sample.

The length of Pxx and the frequency range for w depend on nfft and the values of the input x. The following table indicates the length of Pxx and the frequency range for w in this syntax.

PSD and Frequency Vector Characteristics

[Pxx,w] = pburg(x,p,w) uses a vector of normalized frequencies w with two or more elements to compute the PSD at those frequencies and returns a two-sided PSD.

Pxx = pburg(x,p,nfft,fs)

or

Pxx = pburg(x,p,f,fs) uses the integer FFT length nfft to calculate the PSD vector Pxx or uses the vector of frequencies f in Hz and the sampling frequency fs to compute the two-sided PSD vector Pxx at those frequencies. If you specify nfft as the empty vector [], it uses the default value of 256. If you specify fs as the empty vector [], the sampling frequency fs defaults to 1 Hz. The spectral density produced is calculated in units of power per Hz.

[Pxx,f] = pburg(x,p,nfft,fs)

or

[Pxx,f] = pburg(x,p,f,fs) returns the frequency vector f. In this case, the units for the frequency vector are in Hz. The frequency range for f depends on nfft, fs, and the values of the input x. The length of Pxx is the same as in the table above. The following table indicates the frequency range for f for this syntax.

PSD and Frequency Vector Characteristics with fs Specified

[Pxx,f] = pburg(x,p,nfft,fs,'range') or

[Pxx,w] = pburg(x,p,nfft,'range') specifies the range of frequency values to include in f or w. This syntax is useful when x is real. 'range' can be either:

• 'twosided': Compute the two-sided PSD over the frequency range [0,fs). This is the default for determining the frequency range for complex-valued x.

  • If you specify fs as the empty vector, [], the frequency range is [0,1).

  • If you don't specify fs, the frequency range is [0, 2π).

• 'onesided': Compute the one-sided PSD over the frequency ranges specified for real x. This is the default for determining the frequency range for real-valued x. Note that 'onesided' is not valid if you pass in a vector of frequencies (f or w).

Note   You can put the string argument 'range' anywhere in the input argument list after p.

pburg(...) with no outputs plots the PSD in the current figure window. The frequency range on the plot is the same as the range of output w (or f) for a given set of parameters.

Remarks

The power spectral density is computed as the distribution of power per unit frequency. This algorithm depends on your selecting an appropriate model order for your signal.

Examples

The Burg method estimates the spectral density by fitting an AR prediction model of a given order to the signal, so first generate a signal from an AR (all-pole) model of a given order. Use freqz to check the magnitude of the frequency response of your AR filter. Then, generate the input signal x by filtering white noise through the AR filter. Estimate the PSD of x based on a fourth-order AR prediction model because in this case we know that the original AR system model a has order 4:

% Define AR filter coefficients
a = [1 -2.2137 2.9403 -2.1697 0.9606];

[H,w] = freqz(1,a,256);     % AR filter freq response

% Scale to make one-sided PSD
Hp = plot(w/pi,20*log10(2*abs(H)/(2*pi)),'r'); 
hold on;
randn('state',1);
x = filter(1,a,randn(256,1));           % AR system output
pburg(x,4,511); 
xlabel('Normalized frequency (\times \pi rad/sample)')
ylabel('One-sided PSD (dB/rad/sample)')
legend('PSD of model output','PSD estimate of x')

Algorithm

You can use linear prediction filters to model the second-order statistical characteristics of a signal. The prediction filter output can be used to model the signal when the input is white noise.

The Burg method fits an AR linear prediction filter model of the specified order to the input signal by minimizing (using least squares) the arithmetic mean of the forward and backward prediction errors. The spectral density then is computed from the frequency response of the prediction filter. The AR filter parameters are constrained to satisfy the Levinson-Durbin recursion.

References

[1] Marple, S.L. Digital Spectral Analysis, Englewood Cliffs, NJ, Prentice-Hall, 1987, Chapter7.

[2] Stoica, P., and R.L. Moses, Introduction to Spectral Analysis, Prentice-Hall, 1997.

See Also

arburg, lpc, pcov, peig, periodogram, pmcov, pmtm, pmusic, pwelch, pyulear

parzenwin

pcov