Spectral Analysis

Background Information

The goal of spectral estimation is to describe the distribution (over frequency) of the power contained in a signal, based on a finite set of data. Estimation of power spectra is useful in a variety of applications, including the detection of signals buried in wideband noise.

The power spectral density (PSD) of a stationary random process x(n) is mathematically related to the autocorrelation sequence by the discrete-time Fourier transform. In terms of normalized frequency, this is given by

$$P_{xx}(\omega )=\frac{1}{2\pi }{\displaystyle \sum _{m=-\infty }^{\infty }{R}_{xx}(m)e^{-j\omega m}}.$$

This can be written as a function of physical frequency f (for example, in hertz) by using the relation ω = 2πf / f_s, where f_s is the sampling frequency:

$$P_{xx}(f)=\frac{1}{f_s}{\displaystyle \sum _{m=-\infty }^{\infty }{R}_{xx}(m)e^{-j2\pi mf/f_s}}.$$

The correlation sequence can be derived from the PSD by use of the inverse discrete-time Fourier transform:

$$R_{xx}(m)={\displaystyle \underset{-\pi }{\overset{\pi }{\int }}{P}_{xx}(\omega )e^{j\omega m}d\omega }={\displaystyle \underset{-f_s/2}{\overset{f_s/2}{\int }}{P}_{xx}(f)e^{j2\pi mf/f_s}df}.$$

The average power of the sequence x(n) over the entire Nyquist interval is represented by

$$R_{xx}(0)={\displaystyle \underset{-\pi }{\overset{\pi }{\int }}{P}_{xx}(\omega )d\omega ={\displaystyle \underset{-f_s/2}{\overset{f_s/2}{\int }}{P}_{xx}(f)df}}.$$

The average power of a signal over a particular frequency band [ω₁, ω₂], 0 ≤ ω₁ ≤ ω₂ ≤ π, can be found by integrating the PSD over that band:

$${\overline{P}}_{[{\omega }_1,{\omega }_2]}={\displaystyle {\int }_{{\omega }_1}^{{\omega }_2}{P}_{xx}}(\omega )d\omega ={\displaystyle {\int }_{-{\omega }_2}^{-{\omega }_1}{P}_{xx}}(\omega )d\omega .$$

You can see from the above expression that P_xx(ω) represents the power content of a signal in an infinitesimal frequency band, which is why it is called the power spectral density.

The units of the PSD are power (e.g., watts) per unit of frequency. In the case of P_xx(ω), this is watts/radian/sample or simply watts/radian. In the case of P_xx(f), the units are watts/hertz. Integration of the PSD with respect to frequency yields units of watts, as expected for the average power.

For real–valued signals, the PSD is symmetric about DC, and thus P_xx(ω) for 0 ≤ ω ≤ π is sufficient to completely characterize the PSD. However, to obtain the average power over the entire Nyquist interval, it is necessary to introduce the concept of the one-sided PSD.

The one-sided PSD is given by

$$P_{\text{one-sided}}(\omega )=\{\begin{array}{ll}0,\hfill & -\pi \le \omega <0,\hfill \\ 2P_{xx}(\omega ),\hfill & 0\le \omega \le \pi .\hfill \end{array}$$

The average power of a signal over the frequency band, [ω₁,ω₂] with 0 ≤ ω₁ ≤ ω₂ ≤ π, can be computed using the one-sided PSD as

$${\overline{P}}_{[{\omega }_1,{\omega }_2]}={\displaystyle {\int }_{{\omega }_1}^{{\omega }_2}{P}_{\text{one-sided}}}(\omega )d\omega .$$

Spectral Estimation Method

The various methods of spectrum estimation available in the toolbox are categorized as follows:

• Nonparametric methods

• Parametric methods

• Subspace methods

Nonparametric methods are those in which the PSD is estimated directly from the signal itself. The simplest such method is the periodogram. Other nonparametric techniques such as Welch's method [8], the multitaper method (MTM) reduce the variance of the periodogram.

Parametric methods are those in which the PSD is estimated from a signal that is assumed to be output of a linear system driven by white noise. Examples are the Yule-Walker autoregressive (AR) method and the Burg method. These methods estimate the PSD by first estimating the parameters (coefficients) of the linear system that hypothetically generates the signal. They tend to produce better results than classical nonparametric methods when the data length of the available signal is relatively short. Parametric methods also produce smoother estimates of the PSD than nonparametric methods, but are subject to error from model misspecification.

Subspace methods, also known as high-resolution methods or super-resolution methods, generate frequency component estimates for a signal based on an eigenanalysis or eigendecomposition of the autocorrelation matrix. Examples are the multiple signal classification (MUSIC) method or the eigenvector (EV) method. These methods are best suited for line spectra — that is, spectra of sinusoidal signals — and are effective in the detection of sinusoids buried in noise, especially when the signal to noise ratios are low. The subspace methods do not yield true PSD estimates: they do not preserve process power between the time and frequency domains, and the autocorrelation sequence cannot be recovered by taking the inverse Fourier transform of the frequency estimate.

All three categories of methods are listed in the table below with the corresponding toolbox function names. More information about each function is on the corresponding function reference page. See Parametric Modeling for details about lpc and other parametric estimation functions.

Spectral Estimation Methods/Functions

Method	Description	Functions
Periodogram	Power spectral density estimate	periodogram
Welch	Averaged periodograms of overlapped, windowed signal sections	pwelch, cpsd, tfestimate, mscohere
Multitaper	Spectral estimate from combination of multiple orthogonal windows (or “tapers”)	pmtm
Yule-Walker AR	Autoregressive (AR) spectral estimate of a time-series from its estimated autocorrelation function	pyulear
Burg	Autoregressive (AR) spectral estimation of a time-series by minimization of linear prediction errors	pburg
Covariance	Autoregressive (AR) spectral estimation of a time-series by minimization of the forward prediction errors	pcov
Modified Covariance	Autoregressive (AR) spectral estimation of a time-series by minimization of the forward and backward prediction errors	pmcov
MUSIC	Multiple signal classification	pmusic
Eigenvector	Pseudospectrum estimate	peig