Autoregressive power spectral density estimate — Burg’s method
the power spectral density (PSD) estimate,
pxx = pburg(
of a discrete-time signal,
x, found using Burg’s
x is a vector, it is treated as
a single channel. When
x is a matrix, the PSD
is computed independently for each column and stored in the corresponding
pxx is the
distribution of power per unit frequency. The frequency is expressed
in units of rad/sample.
order is the order of
the autoregressive (AR) model used to produce the PSD estimate.
pxx = pburg(
in the discrete Fourier transform (DFT). For real
even, and (
odd. For complex–valued
nfft. If you omit
or specify it as empty, then
pburg uses a default
DFT length of 256.
returns a frequency vector,
f] = pburg(___,
f, in cycles per unit time. The sampling
fs, is the number of samples per unit time. If the unit of time
is seconds, then
f is in cycles/second (Hz). For real-valued signals,
f spans the interval [0,
nfft is even and [0,
nfft is odd. For complex-valued signals,
f spans the
returns the two-sided AR PSD estimates at the frequencies specified in the vector,
f] = pburg(
f. The vector
f must contain at least two elements,
because otherwise the function interprets it as
nfft. The frequencies in
f are in cycles per unit time. The sampling frequency,
fs, is the number of samples per unit time. If the unit of time is
f is in cycles/second (Hz).
pxxc] = pburg(___,'ConfidenceLevel',
probability × 100%
confidence intervals for the PSD estimate in
pburg(___) with no output arguments
plots the AR PSD estimate in dB per unit frequency in the current
Burg PSD Estimate of an AR(4) Process
Create a realization of an AR(4) wide-sense stationary random process. Estimate the PSD using Burg's method. Compare the PSD estimate based on a single realization to the true PSD of the random process.
Create an AR(4) system function. Obtain the frequency response and plot the PSD of the system.
A = [1 -2.7607 3.8106 -2.6535 0.9238]; [H,F] = freqz(1,A,,1); plot(F,20*log10(abs(H))) xlabel('Frequency (Hz)') ylabel('PSD (dB/Hz)')
Create a realization of the AR(4) random process. Set the random number generator to the default settings for reproducible results. The realization is 1000 samples in length. Assume a sampling frequency of 1 Hz. Use
pburg to estimate the PSD for a 4th-order process. Compare the PSD estimate with the true PSD.
rng default x = randn(1000,1); y = filter(1,A,x); [Pxx,F] = pburg(y,4,1024,1); hold on plot(F,10*log10(Pxx)) legend('True Power Spectral Density','pburg PSD Estimate')
Reflection Coefficients for Model Order Determination
Create a realization of an AR(4) process. Use
arburg to determine the reflection coefficients. Use the reflection coefficients to determine an appropriate AR model order for the process. Obtain an estimate of the process PSD.
Create a realization of an AR(4) process 1000 samples in length. Use
arburg with the order set to 12 to return the reflection coefficients. Plot the reflection coefficients to determine an appropriate model order.
A = [1 -2.7607 3.8106 -2.6535 0.9238]; rng default x = filter(1,A,randn(1000,1)); [a,e,k] = arburg(x,12); stem(k,'filled') title('Reflection Coefficients') xlabel('Model Order')
The reflection coefficients decay to zero after order 4. This indicates an AR(4) model is most appropriate.
Obtain a PSD estimate of the random process using Burg's method. Use 1000 points in the DFT. Plot the PSD estimate.
Burg PSD Estimate of a Multichannel Signal
Create a multichannel signal consisting of three sinusoids in additive white Gaussian noise. The sinusoids' frequencies are 100 Hz, 200 Hz, and 300 Hz. The sampling frequency is 1 kHz, and the signal has a duration of 1 s.
Fs = 1000; t = 0:1/Fs:1-1/Fs; f = [100;200;300]; x = cos(2*pi*f*t)'+randn(length(t),3);
Estimate the PSD of the signal using Burg's method with a 12th-order autoregressive model. Use the default DFT length. Plot the estimate.
morder = 12; pburg(x,morder,,Fs)
C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.
Usage notes and limitations:
nfft argument is variable-size at compile time, then it must
not become a scalar or an empty array at runtime.
Introduced before R2006a