Front page|Spectrum - Spectral Analysis in Python (0.5.2)

previous | next | modules | index

5.4. Other Power Spectral Density estimates

5.4.1. Covariance method

arcovar(x, order)

Simple and fast implementation of the covariance AR estimate

This code is 10 times faster than arcovar_marple() and more importantly only 10 lines of code, compared to a 200 loc for arcovar_marple()

Parameters:
  • X (array) – Array of complex data samples
  • oder (int) – Order of linear prediction model
Returns:

  • a - Array of complex forward linear prediction coefficients
  • e - error

The covariance method fits a Pth order autoregressive (AR) model to the input signal, which is assumed to be the output of an AR system driven by white noise. This method minimizes the forward prediction error in the least-squares sense. The output vector contains the normalized estimate of the AR system parameters

The white noise input variance estimate is also returned.

If is the power spectral density of y(n), then:

\frac{e}{\left| A(e^{jw}) \right|^2} = \frac{e}{\left| 1+\sum_{k-1}^P a(k)e^{-jwk}\right|^2}

Because the method characterizes the input data using an all-pole model, the correct choice of the model order p is important.

from spectrum import *
from pylab import plot, log10, linspace, axis
ar_values, error = arcovar(marple_data, 15)
psd = arma2psd(ar_values, sides='centerdc')
plot(linspace(-0.5, 0.5, len(psd)), 10*log10(psd/max(psd)))
axis([-0.5, 0.5, -60, 0])

[hires.png, pdf]

../_images/1961b7fcb6.png

See also

pcovar

Validation :the AR parameters are the same as those returned by a completely different function arcovar_marple().
References :[Mathworks]
arcovar_marple(x, order)

Estimate AR model parameters using covariance method

This implementation is based on [Marple]. This code is far more complicated and slower than arcovar() function, which is now the official version. See arcovar() for a detailed description of Covariance method.

This function should be used in place of arcovar only if order<=4, for which arcovar() does not work.

Fast algorithm for the solution of the covariance least squares normal equations from Marple.

Parameters:
  • X (array) – Array of complex data samples
  • oder (int) – Order of linear prediction model
Returns:

  • AF - Array of complex forward linear prediction coefficients
  • PF - Real forward linear prediction variance at order IP
  • AB - Array of complex backward linear prediction coefficients
  • PB - Real backward linear prediction variance at order IP
  • PV - store linear prediction coefficients

Note

this code and the original code in Marple diverge for ip>10. it seems that this is related to single precision used with complex type in fortran whereas numpy uses double precision for complex type.

Validation :the AR parameters are the same as those returned by a completely different function arcovar().
References :[Marple]
class pcovar(data, order, NFFT=None, sampling=1.0, scale_by_freq=False)

Class to create PSD based on covariance algorithm

See arcovar() for description.

from spectrum import *
p = pcovar(marple_data, 15, NFFT=4096)
p()
p.plot(sides='centerdc')

[hires.png, pdf]

../_images/731df46588.png

See also

arcovar

Constructor

For a detailled description of the parameters, see arcovar().

Parameters:
  • data (array) – input data (list or numpy.array)
  • order (int) –
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
  • sampling (float) – sampling frequency of the input data.
arcovar_marple(x, order)

Estimate AR model parameters using covariance method

This implementation is based on [Marple]. This code is far more complicated and slower than arcovar() function, which is now the official version. See arcovar() for a detailed description of Covariance method.

This function should be used in place of arcovar only if order<=4, for which arcovar() does not work.

Fast algorithm for the solution of the covariance least squares normal equations from Marple.

Parameters:
  • X (array) – Array of complex data samples
  • oder (int) – Order of linear prediction model
Returns:

  • AF - Array of complex forward linear prediction coefficients
  • PF - Real forward linear prediction variance at order IP
  • AB - Array of complex backward linear prediction coefficients
  • PB - Real backward linear prediction variance at order IP
  • PV - store linear prediction coefficients

Note

this code and the original code in Marple diverge for ip>10. it seems that this is related to single precision used with complex type in fortran whereas numpy uses double precision for complex type.

Validation :the AR parameters are the same as those returned by a completely different function arcovar().
References :[Marple]

5.4.2. Eigen-analysis methods

eigen(X, P, NSIG=None, method='music', threshold=None, NFFT=4096, criteria='aic', verbose=False)

Pseudo spectrum using eigenvector method (EV or Music)

This function computes either the Music or EigenValue (EV) noise subspace frequency estimator.

First, an autocorrelation matrix of order P is computed from the data. Second, this matrix is separated into vector subspaces, one a signal subspace and the other a noise subspace using a SVD method to obtain the eigen values and vectors. From the eigen values \lambda_i, and eigen vectors v_k, the pseudo spectrum (see note below) is computed as follows:

P_{ev}(f) = \frac{1}{e^H(f)\left(\sum\limits_{k=M+1}^{p} \frac{1}{\lambda_k}v_kv_k^H\right)e(f)}

The separation of the noise and signal subspaces requires expertise of the signal. However, AIC and MDL criteria may be used to automatically perform this task.

You still need to provide the parameter P to indicate the maximum number of eigen values to be computed. The criteria will just select a subset to estimate the pseudo spectrum (see aic_eigen() and mdl_eigen() for details.

Note

pseudo spectrum. func:eigen does not compute a PSD estimate. Indeed, the method does not preserve the measured process power.

Parameters:
  • X – Array data samples
  • P (int) – maximum number of eigen values to compute. NSIG (if specified) must therefore be less than P.
  • method (str) – ‘music’ or ‘ev’.
  • NSIG (int) – If specified, the signal sub space uses NSIG eigen values.
  • threshold (float) – If specified, the signal sub space is made of the eigen values larger than \rm{threshold} \times \lambda_{min}, where \lambda_{min} is the minimum eigen values.
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
Returns:

  • PSD: Array of real frequency estimator values (two sided for

    complex data and one sided for real data)

  • S, the eigen values

 
from spectrum import *
from pylab import plot, log10, linspace, legend, axis

psd, ev = eigen(marple_data, 15, NSIG=11)
f = linspace(-0.5, 0.5, len(psd))
plot(f, 10 * log10(psd/max(psd)), label='User defined')
psd, ev = eigen(marple_data, 15, threshold=2)
plot(f, 10 * log10(psd/max(psd)), label='threshold method (100)')
psd, ev = eigen(marple_data, 15)
plot(f, 10 * log10(psd/max(psd)), label='AIC method (8)')
legend()
axis([-0.5, 0.5, -120, 0])

[hires.png, pdf]

../_images/d379f09bbc.png

See also

pev(), pmusic(), aic_eigen()

References :[Marple], Chap 13

Todo

for developers:

  • what should be the second argument of the criteria N, N-P, P...?
  • what should be the max value of NP
class pev(data, IP, NSIG=None, NFFT=None, sampling=1.0)

Class to create PSD using ARMA estimator.

See eigenfre() for description.

from spectrum import *
p = pev(marple_data, 15, NFFT=4096)
p()
p.plot()

[hires.png, pdf]

../_images/87cae96006.png

Constructor:

For a detailed description of the parameters, see arma_estimate().

Parameters:
  • data (array) – input data (list or numpy.array)
  • IP (int) –
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
  • sampling (float) – sampling frequency of the input data.
class pmusic(data, IP, NSIG=None, NFFT=None, sampling=1.0)

Class to create PSD using ARMA estimator.

See pmusic() and eigenfre() for description.

from spectrum import *
p = pmusic(marple_data, 15, NFFT=4096)
p()
p.plot()

[hires.png, pdf]

../_images/4e458ad31d.png

Constructor:

For a detailed description of the parameters, see arma_estimate().

Parameters:
  • data (array) – input data (list or numpy.array)
  • IP (int) –
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
  • sampling (float) – sampling frequency of the input data.

5.4.3. Minimum Variance estimator

Minimum Variance Spectral Estimators

pminvar(data, order[, NFFT, sampling]) Class to create PSD based on the Minimum variance spectral estimation

Code author: Thomas Cokelaer, 2011

minvar(X, order, sampling=1.0, NFFT=4096)

Minimum Variance Spectral Estimation (MV)

This function computes the minimum variance spectral estimate using the Musicus procedure. The Burg algorithm from arburg() is used for the estimation of the autoregressive parameters. The MV spectral estimator is given by:

P_{MV}(f) = \frac{T}{e^H(f) R^{-1}_p e(f)}

where R^{-1}_p is the inverse of the estimated autocorrelation matrix (Toeplitz) and e(f) is the complex sinusoid vector.

Parameters:
  • X – Array of complex or real data samples (length N)
  • order (int) – Dimension of correlation matrix (AR order = order - 1 )
  • T (float) – Sample interval (PSD scaling)
  • NFFT (int) – length of the final PSD
Returns:

  • PSD - Power spectral density values (two-sided)
  • AR - AR coefficients (Burg algorithm)
  • k - Reflection coefficients (Burg algorithm)

Note

The MV spectral estimator is not a true PSD function because the area under the MV estimate does not represent the total power in the measured process. MV minimises the variance of the output of a narrowband filter and adpats itself to the spectral content of the input data at each frequency.

Example :The following example computes a PSD estimate using minvar() The output PSD is transformed to a centerdc PSD and plotted. 
from spectrum import *
from pylab import plot, log10, linspace, xlim
psd, A, k = minvar(marple_data, 15)
psd = twosided_2_centerdc(psd) # switch positive and negative freq
f = linspace(-0.5, 0.5, len(psd))
plot(f, 10 * log10(psd/max(psd)))
xlim(-0.5, 0.5 )

[hires.png, pdf]

../_images/f4fb59ef13.png

See also

  • External functions used are arburg() and numpy.fft.fft
  • pminvar, a Class dedicated to MV method.
Reference :[Marple]
class pminvar(data, order, NFFT=None, sampling=1.0)

Class to create PSD based on the Minimum variance spectral estimation

See minvar() for description.

from spectrum import *
p = pminvar(marple_data, 15, NFFT=4096)
p()
p.plot(sides='centerdc')

[hires.png, pdf]

../_images/82298c7bfa.png

Constructor

For a detailled description of the parameters, see minvar().

Parameters:
  • data (array) – input data (list or numpy.array)
  • order (int) –
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
  • sampling (float) – sampling frequency of the input data.

5.4.4. Modified Covariance method

modcovar(x, order)

Simple and fast implementation of the covariance AR estimate

This code is 10 times faster than modcovar_marple() and more importantly only 10 lines of code, compared to a 200 loc for modcovar_marple()

Parameters:
  • X – Array of complex data samples
  • order (int) – Order of linear prediction model
Returns:

  • P - Real linear prediction variance at order IP
  • A - Array of complex linear prediction coefficients

 
from spectrum import *
from pylab import *

a, p = modcovar(marple_data, 15)
PSD = arma2psd(a)
PSD = cshift(PSD, len(PSD)/2) # switch positive and negative freq
plot(linspace(-0.5, 0.5, 4096), 10*log10(PSD/max(PSD)))
axis([-0.5,0.5,-60,0])

[hires.png, pdf]

../_images/d55b2332e8.png

See also

pmodcovar

Validation :the AR parameters are the same as those returned by a completely different function modcovar_marple().
References :Mathworks
modcovar_marple(X, IP)

Fast algorithm for the solution of the modified covariance least squares normal equations.

This implementation is based on [Marple]. This code is far more complicated and slower than modcovar() function, which is now the official version. See modcovar() for a detailed description of Modified Covariance method.

Parameters:
  • X –
    • Array of complex data samples X(1) through X(N)
  • IP (int) –
    • Order of linear prediction model (integer)
Returns:

  • P - Real linear prediction variance at order IP

  • A - Array of complex linear prediction coefficients

  • ISTAT - Integer status indicator at time of exit
    1. for normal exit (no numerical ill-conditioning)
    2. if P is not a positive value
    3. if DELTA’ and GAMMA’ do not lie in the range 0 to 1
    4. if P’ is not a positive value
    5. if DELTA and GAMMA do not lie in the range 0 to 1
Validation :

the AR parameters are the same as those returned by a completely different function modcovar().

Note

validation. results similar to test example in Marple but starts to differ for ip~8. with ratio of 0.975 for ip=15 probably due to precision.

References :[Marple]
class pmodcovar(data, order, NFFT=None, sampling=1.0, scale_by_freq=False)

Class to create PSD based on modified covariance algorithm

See modcovar() for description.

Examples

from spectrum import *
p = pmodcovar(marple_data, 15, NFFT=4096)
p()
p.plot(sides='centerdc')

[hires.png, pdf]

../_images/b0bd0d17ce.png

See also

modcovar

Constructor

For a detailled description of the parameters, see modcovar().

Parameters:
  • data (array) – input data (list or numpy.array)
  • order (int) –
  • NFFT (int) – total length of the final data sets (padded with zero if needed; default is 4096)
  • sampling (float) – sampling frequency of the input data.

Search

Enter search terms or a module, class or function name.

Authors: thomas-cokelaer.info
This page: Show source.

Table Of Contents

© Copyright Thomas Cokelaer. Last updated on Mar 03, 2012.
Created using Sphinx 1.1.2.
This page uses Google Analytics to collect statistics. You can disable it by blocking the JavaScript coming from www.google-analytics.com.