View source for Spectral density estimation

{{Short description|Signal processing technique}}
{{For|the statistical method|Probability density estimation}}
{{Broader|Spectral density}}
In [[statistical signal processing]], the goal of '''spectral density estimation''' ('''SDE''') or simply '''spectral estimation''' is to [[estimation theory|estimate]] the [[spectral density]] (also known as the [[power spectrum|power spectral density]]) of a signal from a sequence of time samples of the signal.<ref>[[Peter Stoica|P Stoica]] and R Moses, Spectral Analysis of Signals, Prentice Hall, 2005.</ref> Intuitively speaking, the spectral density characterizes the [[frequency]] content of the signal. One purpose of estimating the spectral density is to detect any [[Periodic function|periodicities]] in the data, by observing peaks at the frequencies corresponding to these periodicities.

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 ==
{{merge cleanup|Frequency domain|date=June 2015}}

[[Image:Voice waveform and spectrum.png|thumb|Example of voice waveform and its frequency spectrum]]
[[Image:triangle-td and fd.png|thumb|A periodic waveform ([[triangle wave]]) and its frequency spectrum, showing a "fundamental" frequency at 220 Hz followed by multiples (harmonics) of 220 Hz]]
[[File:Comparison of periodogram and Welch methods of spectral density estimation.png|thumb|The power spectral density of a segment of music is estimated by two different methods, for comparison]]

'''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 (waves)|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 function]]s (such as <math>\sin (t)</math>) 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 (waves)|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 (physics)|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 [[nonlinearity|non-linear]] or [[Time-variant system|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 [[sampling (signal processing)|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 [[Impulse response|filter impulse responses]] and [[Window function|window functions]].  But the periodogram does not provide processing-gain when applied to noiselike signals or even sinusoids at low signal-to-noise ratios{{why|reason=FFT does provide higher process gain as the bin count increases. Why do we suggest that it doesn't at low SNR?|date=November 2024}}. 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]]<ref>{{Citation|last=Welch|first=P. D.|title=The use of Fast Fourier Transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms|journal=IEEE Transactions on Audio and Electroacoustics|volume=AU-15|issue=2|year=1967|pages=70–73|doi=10.1109/TAU.1967.1161901|bibcode=1967ITAE...15...70W|s2cid=13900622 }}</ref>)&nbsp; 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 statistics|non-parametric]],'' ''[[parametric estimation|parametric]],'' and more recently [[semiparametric model|semi-parametric]] (also called sparse) methods.<ref name=":0">{{Cite journal|last1=Stoica|first1=Petre|last2=Babu|first2=Prabhu|last3=Li|first3=Jian|date=January 2011|title=New Method of Sparse Parameter Estimation in Separable Models and Its Use for Spectral Analysis of Irregularly Sampled Data|url=https://ieeexplore.ieee.org/document/5599897|journal=IEEE Transactions on Signal Processing|volume=59|issue=1|pages=35–47|doi=10.1109/TSP.2010.2086452|bibcode=2011ITSP...59...35S |s2cid=15936187 |issn=1053-587X}}</ref> 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 process|stationary stochastic process]] has a certain structure that can be described using a small number of parameters (for example, using an [[autoregressive moving-average model|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<ref>{{Cite book|last1=Stoica|first1=Petre|last2=Li|first2=Jian|last3=Ling|first3=Jun|last4=Cheng|first4=Yubo|title=2009 IEEE International Conference on Acoustics, Speech and Signal Processing |chapter=Missing data recovery via a nonparametric iterative adaptive approach |date=April 2009|chapter-url=http://dx.doi.org/10.1109/icassp.2009.4960347|pages=3369–3372 |publisher=IEEE|doi=10.1109/icassp.2009.4960347|isbn=978-1-4244-2353-8 }}</ref> as well as [[Compressed sensing|signal reconstruction]].

Following is a partial list of spectral density estimation techniques:
* ''Non-parametric methods'' for which the signal samples can be [[unevenly spaced time series|unevenly spaced in time]] ('''records can be incomplete''')
** [[Least-squares spectral analysis]], based on [[least squares]] fitting to known frequencies
** [[Lomb–Scargle periodogram]], an approximation of the [[Least-squares spectral analysis]]
** [[Non-uniform discrete Fourier transform]]

* ''Non-parametric methods'' for which the signal samples must be evenly spaced in time ('''records must be complete'''):
** [[Periodogram]], the [[modulus squared]] of the discrete Fourier transform
** [[Bartlett's method]] is the average of the periodograms taken of multiple segments of the signal to reduce variance of the spectral density estimate
** [[Welch's method]] a windowed version of Bartlett's method that uses overlapping segments
** [[Multitaper]] is a periodogram-based method that uses multiple tapers, or windows, to form independent estimates of the spectral density to reduce variance of the spectral density estimate
** [[Singular spectrum analysis]] is a nonparametric method that uses a [[singular value decomposition]] of the [[covariance matrix]] to estimate the spectral density
** [[Short-time Fourier transform]]
** [[Information field theory#Critical filter|Critical filter]] is a nonparametric method based on [[information field theory]] that can deal with noise, incomplete data, and instrumental response functions
* ''Parametric techniques'' (an incomplete list):
** [[Autoregressive model]] (AR) estimation, which assumes that the ''n''th sample is correlated with the previous ''p'' samples.
** [[Moving-average model]] (MA) estimation, which assumes that the ''n''th sample is correlated with noise terms in the previous ''p'' samples.
** [[Autoregressive moving-average model|Autoregressive moving-average]] (ARMA) estimation, which generalizes the AR and MA models.
** [[MUSIC (algorithm)|MUltiple SIgnal Classification]] (MUSIC) is a popular [[Super-resolution imaging|superresolution]] method.
** [[Estimation of signal parameters via rotational invariance techniques]] (ESPRIT) is another superresolution method.
** [[Maximum entropy spectral estimation]] is an ''all-poles'' method useful for SDE when singular spectral features, such as sharp peaks, are expected.
* ''Semi-parametric techniques'' (an incomplete list):
** SParse Iterative Covariance-based Estimation (SPICE) estimation,<ref name=":0" /> and the more generalized <math>(r,q)</math>-SPICE.<ref>{{Cite book|last1=Sward|first1=Johan|last2=Adalbjornsson|first2=Stefan Ingi|last3=Jakobsson|first3=Andreas|title=2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP) |chapter=A generalization of the sparse iterative covariance-based estimator |date=March 2017|chapter-url=http://dx.doi.org/10.1109/icassp.2017.7952898|pages=3954–3958 |publisher=IEEE|doi=10.1109/icassp.2017.7952898|isbn=978-1-5090-4117-6 |s2cid=5640068 }}</ref> 
** Iterative Adaptive Approach (IAA) estimation.<ref>{{Cite journal|last1=Yardibi|first1=Tarik|last2=Li|first2=Jian|last3=Stoica|first3=Petre|last4=Xue|first4=Ming|last5=Baggeroer|first5=Arthur B.|date=January 2010|title=Source Localization and Sensing: A Nonparametric Iterative Adaptive Approach Based on Weighted Least Squares|url=https://ieeexplore.ieee.org/document/5417172|journal=IEEE Transactions on Aerospace and Electronic Systems|volume=46|issue=1|pages=425–443|doi=10.1109/TAES.2010.5417172|bibcode=2010ITAES..46..425Y |hdl=1721.1/59588 |s2cid=18834345 |issn=0018-9251|hdl-access=free}}</ref> 
** [[Lasso (statistics)|Lasso]], similar to [[least-squares spectral analysis]] but with a sparsity enforcing penalty.<ref>{{Cite book|last1=Panahi|first1=Ashkan|last2=Viberg|first2=Mats|title=2011 International ITG Workshop on Smart Antennas |chapter=On the resolution of the LASSO-based DOA estimation method |date=February 2011|chapter-url=http://dx.doi.org/10.1109/wsa.2011.5741938|pages=1–5 |publisher=IEEE|doi=10.1109/wsa.2011.5741938|isbn=978-1-61284-075-8 |s2cid=7013162 |url=https://research.chalmers.se/en/publication/141363 }}</ref>

=== Parametric estimation ===

In parametric spectral estimation, one assumes that the signal is modeled by a [[stationary process]] which has a spectral density function (SDF) <math>S(f; a_1, \ldots, a_p)</math> that is a function of the frequency <math>f</math> and <math>p</math> parameters <math>a_1, \ldots, a_p</math>.<ref name=Percival1993>{{cite book |last1=Percival|first1=Donald B.|last2=Walden|first2=Andrew T.|title=Spectral Analysis for Physical Applications|date=1992|publisher=Cambridge University Press|isbn=9780521435413}}</ref> 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]] <math>\text{AR}(p)</math> of order <math>p</math>.{{r|Percival1993|page1=392}} A signal sequence <math>\{Y_t\}</math> obeying a zero mean <math>\text{AR}(p)</math>  process satisfies the equation

:<math>Y_t = \phi_1Y_{t-1} + \phi_2Y_{t-2} + \cdots + \phi_pY_{t-p} + \epsilon_t,</math>

where the <math>\phi_1,\ldots,\phi_p</math> are fixed coefficients and <math>\epsilon_t</math> is a white noise process with zero mean and ''innovation variance'' <math>\sigma^2_p</math>. The SDF for this process is

:<math>
  S(f; \phi_1, \ldots, \phi_p, \sigma^2_p) =
    \frac{\sigma^2_p\Delta t}{\left| 1 - \sum_{k=1}^p \phi_k e^{-2i\pi f k \Delta t}\right|^2} \qquad |f| < f_N,
</math>

with <math>\Delta t</math> the sampling time interval and <math>f_N</math> the [[Nyquist frequency]].

There are a number of approaches to estimating the parameters <math>\phi_1, \ldots, \phi_p,\sigma^2_p</math> of the <math>\text{AR}(p)</math> process and thus the spectral density:{{r|Percival1993|page1=452-453}}
* The ''[[Autoregressive model#Yule–Walker equations|Yule–Walker estimators]]'' are found by recursively solving the Yule–Walker equations for an <math>\text{AR}(p)</math> 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.{{r|Percival1993|page1=452}} Burg associated these with [[maximum entropy spectral estimation]].<ref name=Burg>Burg, J.P. (1967) "Maximum Entropy Spectral Analysis", ''Proceedings of the 37th Meeting of the Society of Exploration Geophysicists'', Oklahoma City, Oklahoma.</ref>
* The ''forward-backward least-squares estimators'' treat the <math>\text{AR}(p)</math> 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 [[estimation theory|estimating]] the [[frequency]], amplitude, and phase-shift of a [[digital signal processing|signal]] in the presence of [[noise]] given assumptions about the number of the components.<ref>Hayes, Monson H., ''Statistical Digital Signal Processing and Modeling'', John Wiley & Sons, Inc., 1996. {{ISBN|0-471-59431-8}}.</ref> This contrasts with the general methods above, which do not make prior assumptions about the components.

=== Single tone ===
{{See also|Sinusoidal model}}

If one only wants to estimate the frequency of the single loudest [[pure tone|pure-tone signal]], 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 quasiprobability distribution|Wigner–Ville distribution]] and higher order [[ambiguity function]]s.<ref name=Lerga>{{cite web|last=Lerga|first=Jonatan|title=Overview of Signal Instantaneous Frequency Estimation Methods|url=http://www.fer.unizg.hr/_download/repository/Jonatan_Lerga_-_kvalifikacijski_rad.pdf|publisher=University of Rijeka|access-date=22 March 2014}}</ref>

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 <math>x(n)</math> consists of a sum of <math>p</math> complex exponentials in the presence of [[white noise]], <math>w(n)</math>
:<math>x(n) = \sum_{i=1}^p A_i e^{j n \omega_i} + w(n)</math>.
The power spectral density of <math>x(n)</math> is composed of <math>p</math> [[impulse function]]s in addition to the spectral density function due to noise.

The most common methods for frequency estimation involve identifying the noise [[linear subspace|subspace]] to extract these components. These methods are based on [[eigendecomposition of a matrix|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 harmonic decomposition|Pisarenko's method]], the [[multiple signal classification]] (MUSIC) method, the eigenvector method, and the minimum norm method.

; [[Pisarenko harmonic decomposition|Pisarenko's method]]: <math>\hat{P}_\text{PHD}\left(e^{j \omega}\right) = \frac{1}{\left|\mathbf{e}^H \mathbf{v}_\text{min}\right|^2}</math>
; [[Multiple signal classification|MUSIC]]: <math>\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}</math>,
; Eigenvector method: <math>\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}</math>
; Minimum norm method: <math>\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</math>

== Example calculation ==

Suppose <math>x_n</math>, from <math>n=0</math> to <math>N-1</math> 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):

:<math>\begin{align}
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{align}</math>

The variance of <math>x_n</math> is, for a zero-mean function as above, given by

:<math>\frac{1}{N} \sum_{n=0}^{N-1} x_n^2.</math>

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 <math>N\to \infty.</math> 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.

:<math>\lim_{N \to \infty} \frac{1}{N} \sum_{n=0}^{N-1} x_n^2.</math>

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

:<math>x(t) = \sum_k A_k \sin(2\pi\nu_k t + \phi_k)</math>

and

:<math>\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^T x(t)^2 dt.</math>

The root mean square of <math>\sin</math> is <math>1/\sqrt{2}</math>, so the variance of <math>A_k \sin(2\pi\nu_k t + \phi_k)</math> is <math>\tfrac{1}{2} A_k^2.</math> Hence, the contribution to the average power of <math>x(t)</math> coming from the component with frequency <math>\nu_k</math> is <math>\tfrac{1}{2}A_k^2.</math> All these contributions add up to the average power of <math>x(t).</math>

Then the power as a function of frequency is <math>\tfrac{1}{2}A_k^2,</math> and its statistical [[cumulative distribution function]] <math>S(\nu)</math> will be

:<math>S(\nu) = \sum _ {k : \nu_k < \nu} \frac{1}{2} A_k^2.</math>

<math>S</math> is a [[step function]], monotonically non-decreasing. Its jumps occur at the frequencies of the [[periods (ring)|periodic]] components of <math>x</math>, 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 <math>\tau</math>, we can take the [[covariance]] of <math>x(t)</math> with <math>x(t + \tau)</math>, and define this to be the [[autocorrelation function]] <math>c</math> of the signal (or data) <math>x</math>:

:<math>c(\tau) = \lim_{T\to\infty} \frac{1}{2T} \int_{-T}^T x(t) x(t + \tau) dt.</math>

If it exists, it is an even function of <math>\tau.</math> If the average power is bounded, then <math>c</math> exists everywhere, is finite, and is bounded by <math>c(0),</math> which is the average power or variance of the data.

It can be shown that <math>c</math> can be decomposed into periodic components with the same periods as <math>x</math>:

:<math>c(\tau) = \sum_k \frac{1}{2} A_k^2 \cos(2\pi\nu_k\tau).</math>

This is in fact the spectral decomposition of <math>c</math> over the different frequencies, and is related to the distribution of power of <math>x</math> over the frequencies: the amplitude of a frequency component of <math>c</math> 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.

== See also ==
* [[Multidimensional spectral estimation]]
* [[Periodogram]]
* [[SigSpec]]
* [[Spectrogram]]
* [[Time–frequency analysis]]
* [[Time–frequency representation]]
* [[Whittle likelihood]]
* [[Spectral power distribution]]

== References ==
{{reflist}}

== Further reading ==
* {{cite book
  | last = Porat
  | first = B.
  | title = Digital Processing of Random Signals: Theory & Methods
  | year = 1994
  | publisher = Prentice Hall
  | isbn = 978-0-13-063751-2 }}
* {{cite book
  | last = Priestley
  | first = M.B.
  | title = Spectral Analysis and Time Series
  | year = 1991
  | publisher = Academic Press
  | isbn = 978-0-12-564922-3 }}

* {{cite book
  | last1 = Stoica
  | first1 = P.
  | last2 = Moses
  | first2 = R.
  | title = Spectral Analysis of Signals
  | year = 2005
  | publisher = Prentice Hall
  | isbn = 978-0-13-113956-5 }}

* {{Cite journal | doi = 10.1109/PROC.1982.12433| title = Spectrum estimation and harmonic analysis| journal = Proceedings of the IEEE| volume = 70| issue = 9| pages = 1055–1096| year = 1982| last1 = Thomson | first1 = D. J. | bibcode = 1982IEEEP..70.1055T| citeseerx = 10.1.1.471.1278| s2cid = 290772}}

{{Statistics|analysis}}

[[Category:Statistical signal processing]]
[[Category:Signal estimation]]
[[Category:Frequency-domain analysis]]
[[Category:Spectrum (physical sciences)]]