# peig

Pseudospectrum using eigenvector method

## Syntax

`[S,w] = peig(x,p)[S,w] = peig(x,p,w)[S,w] = peig(...,nfft)[S,f] = peig(x,p,nfft,fs)[S,f] = peig(x,p,f,fs)[S,f] = peig(...,'corr')[S,f] = peig(x,p,nfft,fs,nwin,noverlap)[...] = peig(...,freqrange)[...,v,e] = peig(...)peig(...)`

## Description

`[S,w] = peig(x,p)` implements the eigenvector spectral estimation method and returns `S`, the pseudospectrum estimate of the input signal `x`, and `w`, a vector of normalized frequencies (in rad/sample) at which the pseudospectrum is evaluated. The pseudospectrum is calculated using estimates of the eigenvectors of a correlation matrix associated with the input data `x`, where `x` is specified as either:

• A row or column vector representing one observation of the signal

• A rectangular array for which each row of `x` represents a separate observation of the signal (for example, each row is one output of an array of sensors, as in array processing), such that `x'*x` is an estimate of the correlation matrix

 Note   You can use the output of `corrmtx` to generate such an array `x`.

You can specify the second input argument `p` as either:

• A scalar integer. In this case, the signal subspace dimension is `p`.

• A two-element vector. In this case, `p(2)`, the second element of `p`, represents a threshold that is multiplied by λmin, the smallest estimated eigenvalue of the signal's correlation matrix. Eigenvalues below the threshold λmin × `p(2)` are assigned to the noise subspace. In this case, `p(1)` specifies the maximum dimension of the signal subspace.

 Note:   If the inputs to `peig` are real sinusoids, set the value of `p` to double the number of input signals. If the inputs are complex sinusoids, set `p` equal to the number of inputs.

The extra threshold parameter in the second entry in `p` provides you more flexibility and control in assigning the noise and signal subspaces.

`S` and `w` have the same length. In general, the length of the FFT and the values of the input `x` determine the length of the computed `S` and the range of the corresponding normalized frequencies. The following table indicates the length of `S` (and `w`) and the range of the corresponding normalized frequencies for this syntax.

S Characteristics for an FFT Length of 256 (Default)

Real/Complex Input DataLength of S and wRange of the Corresponding Normalized Frequencies

Real-valued

129

[0, π]

Complex-valued

256

[0, 2π)

`[S,w] = peig(x,p,w)` returns the pseudospectrum in the vector `S` computed at the normalized frequencies specified in vector `w`, which has two or more elements

`[S,w] = peig(...,nfft)` specifies the integer length of the FFT `nfft` used to estimate the pseudospectrum. The default value for `nfft` (entered as an empty vector `[]`) is 256.

The following table indicates the length of `S` and `w`, and the frequency range for `w` for this syntax.

S and Frequency Vector Characteristics

Real/Complex Input Datanfft Even/OddLength of S and wRange of w

Real-valued

Even

(`nfft`/2 + 1)

[0, π]

Real-valued

Odd

(`nfft` + 1)/2

[0, π)

Complex-valued

Even or odd

`nfft`

[0, 2π)

`[S,f] = peig(x,p,nfft,fs)` returns the pseudospectrum in the vector `S` evaluated at the corresponding vector of frequencies `f` (in Hz). You supply the sampling frequency `fs` in Hz. If you specify `fs` with the empty vector `[]`, the sampling frequency defaults to 1 Hz.

The frequency range for `f` depends on `nfft`, `fs`, and the values of the input `x`. The length of `S` (and `f`) is the same as in the S and Frequency Vector Characteristics above. The following table indicates the frequency range for `f` for this syntax.

S and Frequency Vector Characteristics with fs Specified

Real/Complex Input Datanfft Even/OddRange of f

Real-valued

Even

[0, `fs`/2]

Real-valued

Odd

[0, `fs`/2)

Complex-valued

Even or odd

[0, `fs`)

`[S,f] = peig(x,p,f,fs)` returns the pseudospectrum in the vector `S` computed at the frequencies specified in vector `f`, which has two or more elements

`[S,f] = peig(...,'corr')` forces the input argument `x` to be interpreted as a correlation matrix rather than matrix of signal data. For this syntax `x` must be a square matrix, and all of its eigenvalues must be nonnegative.

`[S,f] = peig(x,p,nfft,fs,nwin,noverlap)` allows you to specify `nwin`, a scalar integer indicating a rectangular window length, or a real-valued vector specifying window coefficients. Use the scalar integer `noverlap` in conjunction with `nwin` to specify the number of input sample points by which successive windows overlap. `noverlap` is not used if `x` is a matrix. The default value for `nwin` is 2 × `p(1)` and `noverlap` is `nwin` – 1.

With this syntax, the input data `x` is segmented and windowed before the matrix used to estimate the correlation matrix eigenvalues is formulated. The segmentation of the data depends on `nwin`, `noverlap`, and the form of `x`. Comments on the resulting windowed segments are described in the following table.

Windowed Data Depending on x and nwin

Input data xForm of nwinWindowed Data

Data vector

Scalar

Length is `nwin`

Data vector

Vector of coefficients

Length is `length(nwin)`

Data matrix

Scalar

Data is not windowed.

Data matrix

Vector of coefficients

`length(nwin)` must be the same as the column length of `x`, and `noverlap` is not used.

See the table, Eigenvector Length Depending on Input Data and Syntax, for related information on this syntax.

 Note   The arguments `nwin` and `noverlap` are ignored when you include the string `'corr'` in the syntax.

`[...] = peig(...,freqrange)` specifies the range of frequency values to include in `f` or `w`. This syntax is useful when `x` is real. `freqrange` can be either:

• `'onesided'` — returns the one-sided PSD of a real input signal, `x`. If `nfft` is even, `Pxx` has length`nfft`/2+1 and is computed over the interval [0,π]. If `nfft` is odd, the length of `Pxx` is (`nfft` + 1)/2 and the frequency interval is [0,π). When your specify `fs` , the intervals are [0,`fs`/2) and [0,`fs`/2] for even and odd length`nfft`respectively.

• `'twosided'` — returns the two-sided PSD for either real or complex input, `x`. In this case, `Pxx` has length `nfft` and is computed over the interval [0,2π). When you specify `fs`, the frequency interval is [0,`fs`).

• `'centered'` — returns the centered two-sided PSD for either real or complex input, `x`. In this case, `Pxx` has length `nfft` and is computed over the interval (-π, π] for even length `nfft` and (–π,π) for odd length `nfft`. When you specify `fs`, the frequency intervals are (`-fs`/2, `fs`/2] and (`-fs`/2,`fs`/2) for even and odd length `nfft` respectively.

 Note   You can put the string arguments `freqrange` or `'corr'` anywhere in the input argument list after `p`.

`[...,v,e] = peig(...)` returns the matrix `v` of noise eigenvectors, along with the associated eigenvalues in the vector `e`. The columns of `v` span the noise subspace of dimension `size(v,2)`. The dimension of the signal subspace is `size(v,1)-size(v,2)`. For this syntax, `e` is a vector of estimated eigenvalues of the correlation matrix.

`peig(...)` with no output arguments plots the pseudospectrum in the current figure window.

## Examples

collapse all

### Pseudospectrum of Sum of Sinusoids

Implement the eigenvector method to find the pseudospectrum of the sum of three sinusoids in noise. Use the default FFT length of 256. The inputs are complex sinusoids so you set `p` equal to the number of inputs. Use the modified covariance method for the correlation matrix estimate.

```n = 0:99; s = exp(1i*pi/2*n)+2*exp(1i*pi/4*n)+exp(1i*pi/3*n)+randn(1,100); X = corrmtx(s,12,'mod'); peig(X,3,'whole') ```

collapse all

### Tips

In the process of estimating the pseudospectrum, `peig` computes the noise and signal subspaces from the estimated eigenvectors vj and eigenvalues λj of the signal's correlation matrix. The smallest of these eigenvalues is used in conjunction with the threshold parameter `p(2)` to affect the dimension of the noise subspace in some cases.

The length n of the eigenvectors computed by `peig` is the sum of the dimensions of the signal and noise subspaces. This eigenvector length depends on your input (signal data or correlation matrix) and the syntax you use.

The following table summarizes the dependency of the eigenvector length on the input argument.

Eigenvector Length Depending on Input Data and Syntax

Form of Input Data xComments on the SyntaxLength n of Eigenvectors

Row or column vector

`nwin` is specified as a scalar integer.

`nwin`

Row or column vector

`nwin` is specified as a vector.

`length(nwin)`

Row or column vector

`nwin` is not specified.

2 × `p(1)`

l-by-m matrix

If `nwin` is specified as a scalar, it is not used. If `nwin` is specified as a vector, `length(nwin)` must equal m.

m

m-by-m nonnegative definite matrix

The string `'corr'` is specified and `nwin` is not used.

m

You should specify `nwin` > `p(1)` or `length(nwin)` > `p(1)` if you want `p(2)` > `1` to have any effect.

### Algorithms

The eigenvector method estimates the pseudospectrum from a signal or a correlation matrix using a weighted version of the MUSIC algorithm derived from Schmidt's eigenspace analysis method [1] [2]. The algorithm performs eigenspace analysis of the signal's correlation matrix in order to estimate the signal's frequency content. The eigenvalues and eigenvectors of the signal's correlation matrix are estimated using `svd` if you don't supply the correlation matrix. This algorithm is particularly suitable for signals that are the sum of sinusoids with additive white Gaussian noise.

The eigenvector method produces a pseudospectrum estimate given by

${P}_{ev}\left(f\right)=\frac{1}{\sum _{k=p+1}^{N}|{\text{v}}_{k}^{H}e\left(f\right){|}^{2}/{\lambda }_{k}}$

where N is the dimension of the eigenvectors and vk is the kth eigenvector of the correlation matrix of the input signal. The integer p is the dimension of the signal subspace, so the eigenvectors vk used in the sum correspond to the smallest eigenvalues λk of the correlation matrix. The eigenvectors used span the noise subspace. The vector e(f) consists of complex exponentials, so the inner product

${\text{v}}_{k}^{H}e\left(f\right)$

amounts to a Fourier transform. This is used for computation of the pseudospectrum. The FFT is computed for each vk and then the squared magnitudes are summed and scaled.

## References

[1] Marple, S. Lawrence. Digital Spectral Analysis. Englewood Cliffs, NJ: Prentice-Hall, 1987, pp. 373–378.

[2] Schmidt, R. O. "Multiple Emitter Location and Signal Parameter Estimation." IEEE® Transactions on Antennas and Propagation. Vol. AP-34, March, 1986, pp. 276–280.

[3] Stoica, Petre, and Randolph L. Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.