# Spectral density estimation

Last updated

In statistical signal processing, the goal of spectral density estimation (SDE) is to estimate the spectral density (also known as the power spectral density) of a random signal from a sequence of time samples of the signal. [1] Intuitively speaking, the spectral density characterizes the frequency content of the signal. One purpose of estimating the spectral density is to detect any periodicities in the data, by observing peaks at the frequencies corresponding to these periodicities.

## Contents

Some SDE techniques assume that a signal is composed of a limited (usually small) number of generating frequencies plus noise and seek to find the location and intensity of the generated frequencies. Others make no assumption on the number of components and seek to estimate the whole generating spectrum.

## Overview

Spectrum analysis, also referred to as frequency domain analysis or spectral density estimation, is the technical process of decomposing a complex signal into simpler parts. As described above, many physical processes are best described as a sum of many individual frequency components. Any process that quantifies the various amounts (e.g. amplitudes, powers, intensities) versus frequency (or phase) can be called spectrum analysis.

Spectrum analysis can be performed on the entire signal. Alternatively, a signal can be broken into short segments (sometimes called frames), and spectrum analysis may be applied to these individual segments. Periodic functions (such as ${\displaystyle \sin(t)}$) are particularly well-suited for this sub-division. General mathematical techniques for analyzing non-periodic functions fall into the category of Fourier analysis.

The Fourier transform of a function produces a frequency spectrum which contains all of the information about the original signal, but in a different form. This means that the original function can be completely reconstructed (synthesized) by an inverse Fourier transform. For perfect reconstruction, the spectrum analyzer must preserve both the amplitude and phase of each frequency component. These two pieces of information can be represented as a 2-dimensional vector, as a complex number, or as magnitude (amplitude) and phase in polar coordinates (i.e., as a phasor). A common technique in signal processing is to consider the squared amplitude, or power; in this case the resulting plot is referred to as a power spectrum.

Because of reversibility, the Fourier transform is called a representation of the function, in terms of frequency instead of time; thus, it is a frequency domain representation. Linear operations that could be performed in the time domain have counterparts that can often be performed more easily in the frequency domain. Frequency analysis also simplifies the understanding and interpretation of the effects of various time-domain operations, both linear and non-linear. For instance, only non-linear or time-variant operations can create new frequencies in the frequency spectrum.

In practice, nearly all software and electronic devices that generate frequency spectra utilize a discrete Fourier transform (DFT), which operates on samples of the signal, and which provides a mathematical approximation to the full integral solution. The DFT is almost invariably implemented by an efficient algorithm called fast Fourier transform (FFT). The array of squared-magnitude components of a DFT is a type of power spectrum called periodogram, which is widely used for examining the frequency characteristics of noise-free functions such as filter impulse responses and window functions. But the periodogram does not provide processing-gain when applied to noiselike signals or even sinusoids at low signal-to-noise ratios. In other words, the variance of its spectral estimate at a given frequency does not decrease as the number of samples used in the computation increases. This can be mitigated by averaging over time (Welch's method [2] )  or over frequency (smoothing). Welch's method is widely used for spectral density estimation (SDE). However, periodogram-based techniques introduce small biases that are unacceptable in some applications. So other alternatives are presented in the next section.

## Techniques

Many other techniques for spectral estimation have been developed to mitigate the disadvantages of the basic periodogram. These techniques can generally be divided into non-parametric, parametric, and more recently semi-parametric (also called sparse) methods. [3] The non-parametric approaches explicitly estimate the covariance or the spectrum of the process without assuming that the process has any particular structure. Some of the most common estimators in use for basic applications (e.g. Welch's method) are non-parametric estimators closely related to the periodogram. By contrast, the parametric approaches assume that the underlying stationary stochastic process has a certain structure that can be described using a small number of parameters (for example, using an auto-regressive or moving average model). In these approaches, the task is to estimate the parameters of the model that describes the stochastic process. When using the semi-parametric methods, the underlying process is modeled using a non-parametric framework, with the additional assumption that the number of non-zero components of the model is small (i.e., the model is sparse). Similar approaches may also be used for missing data recovery [4] as well as signal reconstruction.

Following is a partial list of non-parametric spectral density estimation techniques:

Below is a partial list of parametric techniques:

And finally some examples of semi-parametric techniques:

• SParse Iterative Covariance-based Estimation (SPICE) estimation, [3] and the more generalized ${\displaystyle (r,q)}$-SPICE. [5]
• Iterative Adaptive Approach (IAA) estimation. [6]
• Lasso, similar to Least-squares spectral analysis but with a sparsity enforcing penalty. [7]

### Parametric estimation

In parametric spectral estimation, one assumes that the signal is modeled by a stationary process which has a spectral density function (SDF) ${\displaystyle S(f;a_{1},\ldots ,a_{p})}$ that is a function of the frequency ${\displaystyle f}$ and ${\displaystyle p}$ parameters ${\displaystyle a_{1},\ldots ,a_{p}}$. [8] The estimation problem then becomes one of estimating these parameters.

The most common form of parametric SDF estimate uses as a model an autoregressive model ${\displaystyle {\text{AR}}(p)}$ of order ${\displaystyle p}$. [8] :392 A signal sequence ${\displaystyle \{Y_{t}\}}$ obeying a zero mean ${\displaystyle {\text{AR}}(p)}$ process satisfies the equation

${\displaystyle Y_{t}=\phi _{1}Y_{t-1}+\phi _{2}Y_{t-2}+\cdots +\phi _{p}Y_{t-p}+\epsilon _{t},}$

where the ${\displaystyle \phi _{1},\ldots ,\phi _{p}}$ are fixed coefficients and ${\displaystyle \epsilon _{t}}$ is a white noise process with zero mean and innovation variance${\displaystyle \sigma _{p}^{2}}$. The SDF for this process is

${\displaystyle S(f;\phi _{1},\ldots ,\phi _{p},\sigma _{p}^{2})={\frac {\sigma _{p}^{2}\Delta t}{\left|1-\sum _{k=1}^{p}\phi _{k}e^{-2i\pi fk\Delta t}\right|^{2}}}\qquad |f|

with ${\displaystyle \Delta t}$ the sampling time interval and ${\displaystyle f_{N}}$ the Nyquist frequency.

There are a number of approaches to estimating the parameters ${\displaystyle \phi _{1},\ldots ,\phi _{p},\sigma _{p}^{2}}$ of the ${\displaystyle {\text{AR}}(p)}$ process and thus the spectral density: [8] :452-453

• The Yule-Walker estimators are found by recursively solving the Yule-Walker equations for an ${\displaystyle {\text{AR}}(p)}$ process
• The Burg estimators are found by treating the Yule-Walker equations as a form of ordinary least squares problem. The Burg estimators are generally considered superior to the Yule-Walker estimators. [8] :452 Burg associated these with maximum entropy spectral estimation. [9]
• The forward-backward least-squares estimators treat the ${\displaystyle {\text{AR}}(p)}$ process as a regression problem and solves that problem using forward-backward method. They are competitive with the Burg estimators.
• The maximum likelihood estimators estimate the parameters using a maximum likelihood approach. This involves a nonlinear optimization and is more complex than the first three.

Alternative parametric methods include fitting to a moving average model (MA) and to a full autoregressive moving average model (ARMA).

## Frequency estimation

Frequency estimation is the process of estimating the complex frequency components of a signal in the presence of noise given assumptions about the number of the components. [10] This contrasts with the general methods above, which do not make prior assumptions about the components.

### Single tone

If one only wants to estimate the single loudest frequency, one can use a pitch detection algorithm. If the dominant frequency changes over time, then the problem becomes the estimation of the instantaneous frequency as defined in the time–frequency representation. Methods for instantaneous frequency estimation include those based on the Wigner-Ville distribution and higher order ambiguity functions. [11]

If one wants to know all the (possibly complex) frequency components of a received signal (including transmitted signal and noise), one uses a multiple-tone approach.

### Multiple tones

A typical model for a signal ${\displaystyle x(n)}$ consists of a sum of ${\displaystyle p}$ complex exponentials in the presence of white noise, ${\displaystyle w(n)}$

${\displaystyle x(n)=\sum _{i=1}^{p}A_{i}e^{jn\omega _{i}}+w(n)}$.

The power spectral density of ${\displaystyle x(n)}$ is composed of ${\displaystyle p}$ impulse functions in addition to the spectral density function due to noise.

The most common methods for frequency estimation involve identifying the noise subspace to extract these components. These methods are based on eigen decomposition of the autocorrelation matrix into a signal subspace and a noise subspace. After these subspaces are identified, a frequency estimation function is used to find the component frequencies from the noise subspace. The most popular methods of noise subspace based frequency estimation are Pisarenko's method, the multiple signal classification (MUSIC) method, the eigenvector method, and the minimum norm method.

Pisarenko's method
${\displaystyle {\hat {P}}_{\text{PHD}}\left(e^{j\omega }\right)={\frac {1}{\left|\mathbf {e} ^{H}\mathbf {v} _{\text{min}}\right|^{2}}}}$
MUSIC
${\displaystyle {\hat {P}}_{\text{MU}}\left(e^{j\omega }\right)={\frac {1}{\sum _{i=p+1}^{M}\left|\mathbf {e} ^{H}\mathbf {v} _{i}\right|^{2}}}}$,
Eigenvector method
${\displaystyle {\hat {P}}_{\text{EV}}\left(e^{j\omega }\right)={\frac {1}{\sum _{i=p+1}^{M}{\frac {1}{\lambda _{i}}}\left|\mathbf {e} ^{H}\mathbf {v} _{i}\right|^{2}}}}$
Minimum norm method
${\displaystyle {\hat {P}}_{\text{MN}}\left(e^{j\omega }\right)={\frac {1}{\left|\mathbf {e} ^{H}\mathbf {a} \right|^{2}}};\ \mathbf {a} =\lambda \mathbf {P} _{n}\mathbf {u} _{1}}$

## Example calculation

Suppose ${\displaystyle x_{n}}$, from ${\displaystyle n=0}$ to ${\displaystyle N-1}$ is a time series (discrete time) with zero mean. Suppose that it is a sum of a finite number of periodic components (all frequencies are positive):

{\displaystyle {\begin{aligned}x_{n}&=\sum _{k}A_{k}\sin(2\pi \nu _{k}n+\phi _{k})\\&=\sum _{k}A_{k}\left(\sin(\phi _{k})\cos(2\pi \nu _{k}n)+\cos(\phi _{k})\sin(2\pi \nu _{k}n)\right)\\&=\sum _{k}\left(\overbrace {a_{k}} ^{A_{k}\sin(\phi _{k})}\cos(2\pi \nu _{k}n)+\overbrace {b_{k}} ^{A_{k}\cos(\phi _{k})}\sin(2\pi \nu _{k}n)\right)\end{aligned}}}

The variance of ${\displaystyle x_{n}}$ is, for a zero-mean function as above, given by

${\displaystyle {\frac {1}{N}}\sum _{n=0}^{N-1}x_{n}^{2}.}$

If these data were samples taken from an electrical signal, this would be its average power (power is energy per unit time, so it is analogous to variance if energy is analogous to the amplitude squared).

Now, for simplicity, suppose the signal extends infinitely in time, so we pass to the limit as ${\displaystyle N\to \infty .}$ If the average power is bounded, which is almost always the case in reality, then the following limit exists and is the variance of the data.

${\displaystyle \lim _{N\to \infty }{\frac {1}{N}}\sum _{n=0}^{N-1}x_{n}^{2}.}$

Again, for simplicity, we will pass to continuous time, and assume that the signal extends infinitely in time in both directions. Then these two formulas become

${\displaystyle x(t)=\sum _{k}A_{k}\sin(2\pi \nu _{k}t+\phi _{k})}$

and

${\displaystyle \lim _{T\to \infty }{\frac {1}{2T}}\int _{-T}^{T}x(t)^{2}dt.}$

The root mean square of ${\displaystyle \sin }$ is ${\displaystyle 1/{\sqrt {2}}}$, so the variance of ${\displaystyle A_{k}\sin(2\pi \nu _{k}t+\phi _{k})}$ is ${\displaystyle {\tfrac {1}{2}}A_{k}^{2}.}$ Hence, the contribution to the average power of ${\displaystyle x(t)}$ coming from the component with frequency ${\displaystyle \nu _{k}}$ is ${\displaystyle {\tfrac {1}{2}}A_{k}^{2}.}$ All these contributions add up to the average power of ${\displaystyle x(t).}$

Then the power as a function of frequency is ${\displaystyle {\tfrac {1}{2}}A_{k}^{2},}$ and its statistical cumulative distribution function ${\displaystyle S(\nu )}$ will be

${\displaystyle S(\nu )=\sum _{k:\nu _{k}<\nu }{\frac {1}{2}}A_{k}^{2}.}$

${\displaystyle S}$ is a step function, monotonically non-decreasing. Its jumps occur at the frequencies of the periodic components of ${\displaystyle x}$, and the value of each jump is the power or variance of that component.

The variance is the covariance of the data with itself. If we now consider the same data but with a lag of ${\displaystyle \tau }$, we can take the covariance of ${\displaystyle x(t)}$ with ${\displaystyle x(t+\tau )}$, and define this to be the autocorrelation function ${\displaystyle c}$ of the signal (or data) ${\displaystyle x}$:

${\displaystyle c(\tau )=\lim _{T\to \infty }{\frac {1}{2T}}\int _{-T}^{T}x(t)x(t+\tau )dt.}$

If it exists, it is an even function of ${\displaystyle \tau .}$ If the average power is bounded, then ${\displaystyle c}$ exists everywhere, is finite, and is bounded by ${\displaystyle c(0),}$ which is the average power or variance of the data.

It can be shown that ${\displaystyle c}$ can be decomposed into periodic components with the same periods as ${\displaystyle x}$:

${\displaystyle c(\tau )=\sum _{k}{\frac {1}{2}}A_{k}^{2}\cos(2\pi \nu _{k}\tau ).}$

This is in fact the spectral decomposition of ${\displaystyle c}$ over the different frequencies, and is related to the distribution of power of ${\displaystyle x}$ over the frequencies: the amplitude of a frequency component of ${\displaystyle c}$ is its contribution to the average power of the signal.

The power spectrum of this example is not continuous, and therefore does not have a derivative, and therefore this signal does not have a power spectral density function. In general, the power spectrum will usually be the sum of two parts: a line spectrum such as in this example, which is not continuous and does not have a density function, and a residue, which is absolutely continuous and does have a density function.

## Related Research Articles

In probability theory, a normaldistribution is a type of continuous probability distribution for a real-valued random variable. The general form of its probability density function is

A chirp is a signal in which the frequency increases (up-chirp) or decreases (down-chirp) with time. In some sources, the term chirp is used interchangeably with sweep signal. It is commonly applied to sonar, radar, and laser systems, and to other applications, such as in spread-spectrum communications. This signal type is biologically inspired and occurs as a phenomenon due to dispersion. It is usually compensated for by using a matched filter, which can be part of the propagation channel. Depending on the specific performance measure, however, there are better techniques both for radar and communication. Since it was used in radar and space, it has been adopted also for communication standards. For automotive radar applications, it is usually called linear frequency modulated waveform (LFMW).

The power spectrum of a time series describes the distribution of power into frequency components composing that signal. According to Fourier analysis, any physical signal can be decomposed into a number of discrete frequencies, or a spectrum of frequencies over a continuous range. The statistical average of a certain signal or sort of signal as analyzed in terms of its frequency content, is called its spectrum.

In signal processing, a periodogram is an estimate of the spectral density of a signal. The term was coined by Arthur Schuster in 1898. Today, the periodogram is a component of more sophisticated methods. It is the most common tool for examining the amplitude vs frequency characteristics of FIR filters and window functions. FFT spectrum analyzers are also implemented as a time-sequence of periodograms.

In statistics, a mixture model is a probabilistic model for representing the presence of subpopulations within an overall population, without requiring that an observed data set should identify the sub-population to which an individual observation belongs. Formally a mixture model corresponds to the mixture distribution that represents the probability distribution of observations in the overall population. However, while problems associated with "mixture distributions" relate to deriving the properties of the overall population from those of the sub-populations, "mixture models" are used to make statistical inferences about the properties of the sub-populations given only observations on the pooled population, without sub-population identity information.

Array processing is a wide area of research in the field of signal processing that extends from the simplest form of 1 dimensional line arrays to 2 and 3 dimensional array geometries. Array structure can be defined as a set of sensors that are spatially separated, e.g. radio antenna and seismic arrays. The sensors used for a specific problem may vary widely, for example microphones, accelerometers and telescopes. However, many similarities exist, the most fundamental of which may be an assumption of wave propagation. Wave propagation means there is a systemic relationship between the signal received on spatially separated sensors. By creating a physical model of the wave propagation, or in machine learning applications a training data set, the relationships between the signals received on spatially separated sensors can be leveraged for many applications.

Estimation theory is a branch of statistics that deals with estimating the values of parameters based on measured empirical data that has a random component. The parameters describe an underlying physical setting in such a way that their value affects the distribution of the measured data. An estimator attempts to approximate the unknown parameters using the measurements. Read In estimation theory, two approaches are generally considered.

In statistics, a probit model is a type of regression where the dependent variable can take only two values, for example married or not married. The word is a portmanteau, coming from probability + unit. The purpose of the model is to estimate the probability that an observation with particular characteristics will fall into a specific one of the categories; moreover, classifying observations based on their predicted probabilities is a type of binary classification model.

In radiometry, radiant flux or radiant power is the radiant energy emitted, reflected, transmitted or received, per unit time, and spectral flux or spectral power is the radiant flux per unit frequency or wavelength, depending on whether the spectrum is taken as a function of frequency or of wavelength. The SI unit of radiant flux is the watt (W), that is the joule per second in SI base units, while that of spectral flux in frequency is the watt per hertz and that of spectral flux in wavelength is the watt per metre —commonly the watt per nanometre.

MUSIC is an algorithm used for frequency estimation and radio direction finding.

Geophysical survey is the systematic collection of geophysical data for spatial studies. Detection and analysis of the geophysical signals forms the core of Geophysical signal processing. The magnetic and gravitational fields emanating from the Earth's interior hold essential information concerning seismic activities and the internal structure. Hence, detection and analysis of the electric and Magnetic fields is very crucial. As the Electromagnetic and gravitational waves are multi-dimensional signals, all the 1-D transformation techniques can be extended for the analysis of these signals as well. Hence this article also discusses multi-dimensional signal processing techniques.

In signal processing, the multitaper method is a technique developed by David J. Thomson to estimate the power spectrum SX of a stationary ergodic finite-variance random process X, given a finite contiguous realization of X as data. It is one of a number of approaches to spectral density estimation.

The method of reassignment is a technique for sharpening a time-frequency representation by mapping the data to time-frequency coordinates that are nearer to the true region of support of the analyzed signal. The method has been independently introduced by several parties under various names, including method of reassignment, remapping, time-frequency reassignment, and modified moving-window method. In the case of the spectrogram or the short-time Fourier transform, the method of reassignment sharpens blurry time-frequency data by relocating the data according to local estimates of instantaneous frequency and group delay. This mapping to reassigned time-frequency coordinates is very precise for signals that are separable in time and frequency with respect to the analysis window.

Least-squares spectral analysis (LSSA) is a method of estimating a frequency spectrum, based on a least squares fit of sinusoids to data samples, similar to Fourier analysis. Fourier analysis, the most used spectral method in science, generally boosts long-periodic noise in long gapped records; LSSA mitigates such problems. An intelligible derivation of the relationship between the Discrete Fourier Transform and signal modeling using the least squares method is given in.

Bilinear time–frequency distributions, or quadratic time–frequency distributions, arise in a sub-field of signal analysis and signal processing called time–frequency signal processing, and, in the statistical analysis of time series data. Such methods are used where one needs to deal with a situation where the frequency composition of a signal may be changing over time; this sub-field used to be called time–frequency signal analysis, and is now more often called time–frequency signal processing due to the progress in using these methods to a wide range of signal-processing problems.

In statistics, signal processing, and time series analysis, a sinusoidal model to approximate a sequence Yi is:

Kernel density estimation is a nonparametric technique for density estimation i.e., estimation of probability density functions, which is one of the fundamental questions in statistics. It can be viewed as a generalisation of histogram density estimation with improved statistical properties. Apart from histograms, other types of density estimators include parametric, spline, wavelet and Fourier series. Kernel density estimators were first introduced in the scientific literature for univariate data in the 1950s and 1960s and subsequently have been widely adopted. It was soon recognised that analogous estimators for multivariate data would be an important addition to multivariate statistics. Based on research carried out in the 1990s and 2000s, multivariate kernel density estimation has reached a level of maturity comparable to its univariate counterparts.

The Kolmogorov–Zurbenko (KZ) filter was first proposed by A. N. Kolmogorov and formally defined by Zurbenko. It is a series of iterations of a moving average filter of length m, where m is a positive, odd integer. The KZ filter belongs to the class of low-pass filters. The KZ filter has two parameters, the length m of the moving average window and the number of iterations k of the moving average itself. It also can be considered as a special window function designed to eliminate spectral leakage.

In statistics, the variance function is a smooth function which depicts the variance of a random quantity as a function of its mean. The variance function is a measure of heteroscedasticity and plays a large role in many settings of statistical modelling. It is a main ingredient in the generalized linear model framework and a tool used in non-parametric regression, semiparametric regression and functional data analysis. In parametric modeling, variance functions take on a parametric form and explicitly describe the relationship between the variance and the mean of a random quantity. In a non-parametric setting, the variance function is assumed to be a smooth function.

Power spectral estimation forms the basis for distinguishing and tracking signals in the presence of noise and extracting information from available data. One dimensional signals are expressed in terms of a single domain while multidimensional signals are represented in wave vector and frequency spectrum. Therefore, spectral estimation in the case of multidimensional signals gets a bit tricky.

## References

1. P Stoica and R Moses, Spectral Analysis of Signals, Prentice Hall, 2005.
2. Welch, P. D. (1967), "The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms", IEEE Transactions on Audio and Electroacoustics, AU-15 (2): 70–73, Bibcode:1967ITAE...15...70W, doi:10.1109/TAU.1967.1161901
3. Stoica, Petre; Babu, Prabhu; Li, Jian (January 2011). "New Method of Sparse Parameter Estimation in Separable Models and Its Use for Spectral Analysis of Irregularly Sampled Data". IEEE Transactions on Signal Processing. 59 (1): 35–47. doi:10.1109/TSP.2010.2086452. ISSN   1053-587X.
4. Stoica, Petre; Li, Jian; Ling, Jun; Cheng, Yubo (April 2009). "Missing data recovery via a nonparametric iterative adaptive approach". 2009 IEEE International Conference on Acoustics, Speech and Signal Processing. IEEE. doi:10.1109/icassp.2009.4960347.
5. Sward, Johan; Adalbjornsson, Stefan Ingi; Jakobsson, Andreas (March 2017). "A generalization of the sparse iterative covariance-based estimator". 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE. doi:10.1109/icassp.2017.7952898.
6. Yardibi, Tarik; Li, Jian; Stoica, Petre; Xue, Ming; Baggeroer, Arthur B. (January 2010). "Source Localization and Sensing: A Nonparametric Iterative Adaptive Approach Based on Weighted Least Squares". IEEE Transactions on Aerospace and Electronic Systems. 46 (1): 425–443. doi:10.1109/TAES.2010.5417172. ISSN   0018-9251.
7. Panahi, Ashkan; Viberg, Mats (February 2011). "On the resolution of The LASSO-based DOA estimation method". 2011 International ITG Workshop on Smart Antennas. IEEE. doi:10.1109/wsa.2011.5741938.
8. Percival, Donald B.; Walden, Andrew T. (1992). Spectral Analysis for Physical Applications. Cambridge University Press. ISBN   9780521435413.
9. Burg, J.P. (1967) "Maximum Entropy Spectral Analysis", Proceedings of the 37th Meeting of the Society of Exploration Geophysicists, Oklahoma City, Oklahoma.
10. Hayes, Monson H., Statistical Digital Signal Processing and Modeling, John Wiley & Sons, Inc., 1996. ISBN   0-471-59431-8.
11. Lerga, Jonatan. "Overview of Signal Instantaneous Frequency Estimation Methods" (PDF). University of Rijeka. Retrieved 22 March 2014.
• Porat, B. (1994). Digital Processing of Random Signals: Theory & Methods. Prentice Hall. ISBN   978-0-13-063751-2.
• Priestley, M.B. (1991). Spectral Analysis and Time Series. Academic Press. ISBN   978-0-12-564922-3.