5.3. Parametric methods ¶

5.3.1. Power Spectrum Density based on Parametric Methods ¶

5.3.1.1. ARMA and MA estimates (yule-walker)¶

ARMA and MA estimates, ARMA and MA PSD estimates.

arma2psd(A=None, B=None, rho=1.0, T=1.0, NFFT=4096, sides='default', norm=False)[source]¶

Computes power spectral density given ARMA values.

This function computes the power spectral density values given the ARMA parameters of an ARMA model. It assumes that the driving sequence is a white noise process of zero mean and variance $\rho_w$ . The sampling frequency and noise variance are used to scale the PSD output, which length is set by the user with the NFFT parameter.

Parameters:

A (array) – Array of AR parameters (complex or real)
B (array) – Array of MA parameters (complex or real)
rho (float) – White noise variance to scale the returned PSD
T (float) – Sampling frequency in Hertz to scale the PSD.
NFFT (int) – Final size of the PSD
sides (str) – Default PSD is two-sided, but sides can be set to centerdc.

Warning

By convention, the AR or MA arrays does not contain the A0=1 value.

If B is None, the model is a pure AR model. If A is None, the model is a pure MA model.

Returns:: two-sided PSD

Details:

AR case: the power spectral density is:

$P_{ARMA}(f) = T \rho_w \left|\frac{B(f)}{A(f)}\right|^2$

where:

$A(f) = 1 + \sum_{k=1}^q b(k) e^{-j2\pi fkT}$

$B(f) = 1 + \sum_{k=1}^p a(k) e^{-j2\pi fkT}$

Example:

import spectrum.arma
from pylab import plot, log10, legend
plot(10*log10(spectrum.arma.arma2psd([1,0.5],[0.5,0.5])), label='ARMA(2,2)')
plot(10*log10(spectrum.arma.arma2psd([1,0.5],None)), label='AR(2)')
plot(10*log10(spectrum.arma.arma2psd(None,[0.5,0.5])), label='MA(2)')
legend()

(Source code, png, hires.png, pdf)

References:: [Marple]

arma_estimate(X, P, Q, lag)[source]¶

Autoregressive and moving average estimators.

This function provides an estimate of the autoregressive parameters, the moving average parameters, and the driving white noise variance of an ARMA(P,Q) for a complex or real data sequence.

The parameters are estimated using three steps:

Estimate the AR parameters from the original data based on a least squares modified Yule-Walker technique,

Produce a residual time sequence by filtering the original data with a filter based on the AR parameters,

Estimate the MA parameters from the residual time sequence.

Parameters:

X (array) – Array of data samples (length N)
P (int) – Desired number of AR parameters
Q (int) – Desired number of MA parameters
lag (int) – Maximum lag to use for autocorrelation estimates

Returns:

A - Array of complex P AR parameter estimates
B - Array of complex Q MA parameter estimates
RHO - White noise variance estimate

Note

lag must be >= Q (MA order)

dependencies:

from spectrum import arma_estimate, arma2psd, marple_data
import pylab

a,b, rho = arma_estimate(marple_data, 15, 15, 30)
psd = arma2psd(A=a, B=b, rho=rho, sides='centerdc', norm=True)
pylab.plot(10 * pylab.log10(psd))
pylab.ylim([-50,0])

(Source code, png, hires.png, pdf)

Reference:: [Marple]

ma(X, Q, M)[source]¶

Moving average estimator.

This program provides an estimate of the moving average parameters and driving noise variance for a data sequence based on a long AR model and a least squares fit.

Parameters:

X (array) – The input data array
Q (int) – Desired MA model order (must be >0 and <M)
M (int) – Order of “long” AR model (suggest at least 2*Q )

Returns:

MA - Array of Q complex MA parameter estimates
RHO - Real scalar of white noise variance estimate

from spectrum import arma2psd, ma, marple_data
import pylab

# Estimate 15 Ma parameters
b, rho = ma(marple_data, 15, 30)
# Create the PSD from those MA parameters
psd = arma2psd(B=b, rho=rho, sides='centerdc')
# and finally plot the PSD
pylab.plot(pylab.linspace(-0.5, 0.5, 4096), 10 * pylab.log10(psd/max(psd)))
pylab.axis([-0.5, 0.5, -30, 0])

(Source code, png, hires.png, pdf)

Reference:: [Marple]

class parma(data, P, Q, lag, NFFT=None, sampling=1.0, scale_by_freq=False)[source]¶

Class to create PSD using ARMA estimator.

See arma_estimate() for description.

from spectrum import parma, marple_data
p = parma(marple_data, 4, 4, 30, NFFT=4096)
p.plot(sides='centerdc')

(Source code, png, hires.png, pdf)

Constructor:

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

Parameters:

data (array) – input data (list or numpy.array)
P (int)
Q (int)
lag (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 pma(data, Q, M, NFFT=None, sampling=1.0, scale_by_freq=False)[source]¶

Class to create PSD using MA estimator.

See ma() for description.

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

(Source code, png, hires.png, pdf)

Constructor:

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

Parameters:

data (array) – input data (list or numpy.array)
Q (int) – MA order
M (int) – AR model used to estimate the MA parameters
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.3.1.2. AR estimate based on Burg algorithm ¶

BURG method of AR model estimate

arburg(X, order, criteria=None)[source]¶

Estimate the complex autoregressive parameters by the Burg algorithm.

$x(n) = \sqrt{(v}) e(n) + \sum_{k=1}^{P+1} a(k) x(n-k)$

Parameters:

x – Array of complex data samples (length N)
order – Order of autoregressive process (0<order<N)
criteria – select a criteria to automatically select the order

Returns:

A Array of complex autoregressive parameters A(1) to A(order). First value (unity) is not included !!
P Real variable representing driving noise variance (mean square of residual noise) from the whitening operation of the Burg filter.
reflection coefficients defining the filter of the model.

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

AR, P, k = arburg(marple_data, 15)
PSD = arma2psd(AR, sides='centerdc')
plot(linspace(-0.5, 0.5, len(PSD)), 10*log10(PSD/max(PSD)))
axis([-0.5,0.5,-60,0])

(Source code, png, hires.png, pdf)

Note

no detrend. Should remove the mean trend to get PSD. Be careful if presence of large mean.
If you don’t know what the order value should be, choose the criterion=’AKICc’, which has the least bias and best resolution of model-selection criteria.

Note

real and complex results double-checked versus octave using complex 64 samples stored in marple_data. It does not agree with Marple fortran routine but this is due to the simplex precision of complex data in fortran.

Reference:: [Marple] [octave]

class pburg(data, order, criteria=None, NFFT=None, sampling=1.0, scale_by_freq=False)[source]¶

Class to create PSD based on Burg algorithm

See arburg() for description.

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

(Source code, png, hires.png, pdf)

Another example based on a real data set is shown here below. Note here that we set the scale_by_freq value to False and True. False should give results equivalent to octave or matlab convention while setting to True implies that the data is multiplied by $2\pi df$ where $df = sampling / N$ .

from spectrum import data_two_freqs, pburg
p = pburg(data_two_freqs(), 7, NFFT=4096)
p.plot()
p = pburg(data_two_freqs(), 7, NFFT=4096, scale_by_freq=True)
p.plot()
from pylab import legend
legend(["un-scaled", "scaled by 2 pi df"])

(Source code, png, hires.png, pdf)

Constructor

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

Parameters:

data (array) – input data (list or np.array)
order (int)
criteria (str)
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.3.1.3. AR estimate based on YuleWalker ¶

Yule Walker method to estimate AR values.

aryule(X, order, norm='biased', allow_singularity=True)[source]¶

Compute AR coefficients using Yule-Walker method

Parameters:

X – Array of complex data values, X(1) to X(N)
order (int) – Order of autoregressive process to be fitted (integer)
norm (str) – Use a biased or unbiased correlation.
allow_singularity (bool)

Returns:

AR coefficients (complex)
variance of white noise (Real)
reflection coefficients for use in lattice filter

Description:

The Yule-Walker method returns the polynomial A corresponding to the AR parametric signal model estimate of vector X using the Yule-Walker (autocorrelation) method. The autocorrelation may be computed using a biased or unbiased estimation. In practice, the biased estimate of the autocorrelation is used for the unknown true autocorrelation. Indeed, an unbiased estimate may result in nonpositive-definite autocorrelation matrix. So, a biased estimate leads to a stable AR filter. The following matrix form represents the Yule-Walker equations. The are solved by means of the Levinson-Durbin recursion:

$\left( \begin{array}{cccc} r(1) & r(2)^* & \dots & r(n)^*\\ r(2) & r(1)^* & \dots & r(n-1)^*\\ \dots & \dots & \dots & \dots\\ r(n) & \dots & r(2) & r(1) \end{array} \right) \left( \begin{array}{cccc} a(2)\\ a(3) \\ \dots \\ a(n+1) \end{array} \right) = \left( \begin{array}{cccc} -r(2)\\ -r(3) \\ \dots \\ -r(n+1) \end{array} \right)$

The outputs consists of the AR coefficients, the estimated variance of the white noise process, and the reflection coefficients. These outputs can be used to estimate the optimal order by using criteria.

Examples:

From a known AR process or order 4, we estimate those AR parameters using the aryule function.

>>> from scipy.signal import lfilter
>>> from spectrum import *
>>> from numpy.random import randn
>>> A  =[1, -2.7607, 3.8106, -2.6535, 0.9238]
>>> noise = randn(1, 1024)
>>> y = lfilter([1], A, noise);
>>> #filter a white noise input to create AR(4) process
>>> [ar, var, reflec] = aryule(y[0], 4)
>>> # ar should contains values similar to A

The PSD estimate of a data samples is computed and plotted as follows:

from spectrum import *
from pylab import *

ar, P, k = aryule(marple_data, 15, norm='biased')
psd = arma2psd(ar)
plot(linspace(-0.5, 0.5, 4096), 10 * log10(psd/max(psd)))
axis([-0.5, 0.5, -60, 0])

(Source code, png, hires.png, pdf)

Note

The outputs have been double checked against (1) octave outputs (octave has norm=’biased’ by default) and (2) Marple test code.

5.3.2. Criteria ¶

Criteria for parametric methods.

This module provides criteria to automatically select order in parametric PSD estimate or pseudo spectrum estimates (e.g, music).

Some criteria such as the AIC criterion helps to chose the order of PSD models such as the ARMA model. Nevertheless, it is difficult to estimate correctly the order of an ARMA model even by using these criteria. The reason being that even the Akaike criteria (AIC) does not provide the proper order with a probability of 1 with infinite samples.

The order choice is related to an expertise of the signal. There is no exact criteria. However, they may provide useful information.

AIC, AICc, KIC and AKICc are based on information theory. They attempt to balance the complexity (or length) of the model against how well the model fits the data. AIC and KIC are biased estimates of the asymmetric and the symmetric Kullback-Leibler divergence respectively. AICc and AKICc attempt to correct the bias.

There are also criteria related to eigen analysis, which takes as input the eigen values of any PSD estimate method.

Example

from spectrum import aryule, AIC, marple_data
from pylab import plot, arange
order = arange(1, 25)
rho = [aryule(marple_data, i, norm='biased')[1] for i in order]
plot(order, AIC(len(marple_data), rho, order), label='AIC')

(Source code, png, hires.png, pdf)

References:: bd-Krim Seghouane and Maiza Bekara “A small sample model selection criterion based on Kullback’s symmetric divergence”, IEEE Transactions on Signal Processing, Vol. 52(12), pp 3314-3323, Dec. 2004

AIC(N, rho, k)[source]¶

Akaike Information Criterion

Parameters:

rho – rho at order k
N – sample size
k – AR order.

If k is the AR order and N the size of the sample, then Akaike criterion is

$AIC(k) = \log(\rho_k) + 2\frac{k+1}{N}$

AIC(64, [0.5,0.3,0.2], [1,2,3])

Validation:: double checked versus octave.

AICc(N, rho, k, norm=True)[source]¶

corrected Akaike information criterion

$AICc(k) = log(\rho_k) + 2 \frac{k+1}{N-k-2}$

Validation:: double checked versus octave.

AKICc(N, rho, k)[source]¶: approximate corrected Kullback information

$AKICc(k) = log(rho_k) + \frac{p}{N*(N-k)} + (3-\frac{k+2}{N})*\frac{k+1}{N-k-2}$

CAT(N, rho, k)[source]¶: Criterion Autoregressive Transfer Function :

$CAT(k) = \frac{1}{N} \sum_{i=1}^k \frac{1}{\rho_i} - \frac{\rho_i}{\rho_k}$

Todo

validation

class Criteria(name, N)[source]¶

Criteria class for an automatic selection of ARMA order.

Available criteria are


AIC	see `AIC()`
AICc	see `AICc()`
KIC	see `KIC()`
AKICc	see `AKICc()`
FPE	see `FPE()`
MDL	see `MDL()`
CAT	see `_CAT()`

Create a criteria object

Parameters:

name – a string or list of strings containing valid criteria method’s name
N (int) – size of the data sample.

property N¶: Getter/Setter for N

property data¶: Getter/Setter for the criteria output

error_incorrect_name = "Invalid name provided. Correct names are ['AIC', 'AICc', 'KIC', 'FPE', 'AKICc', 'MDL'] "¶

error_no_criteria_found = "No names match the valid criteria names (['AIC', 'AICc', 'KIC', 'FPE', 'AKICc', 'MDL'])"¶

property k¶: Getter for k the order of evaluation

property name¶: Getter/Setter for the criteria name

property old_data¶: Getter/Setter for the previous value

property rho¶: Getter/Setter for rho

valid_criteria_names = ['AIC', 'AICc', 'KIC', 'FPE', 'AKICc', 'MDL']¶

FPE(N, rho, k=None)[source]¶

Final prediction error criterion

$FPE(k) = \frac{N + k + 1}{N - k - 1} \rho_k$

Validation:: double checked versus octave.

KIC(N, rho, k)[source]¶

Kullback information criterion

$KIC(k) = log(\rho_k) + 3 \frac{k+1}{N}$

Validation:: double checked versus octave.

MDL(N, rho, k)[source]¶

Minimum Description Length

$MDL(k) = N log \rho_k + p \log N$

Validation:: results

aic_eigen(s, N)[source]¶

AIC order-selection using eigen values

Parameters:

s – a list of p sorted eigen values
N – the size of the input data. To be defined precisely.

Returns:

an array containing the AIC values

Given $n$ sorted eigen values $\lambda_i$ with $0 <= i < n$ , the proposed criterion from Wax and Kailath (1985) is:

$AIC(k) = -2(n-k)N \ln \frac{g(k)}{a(k)} + 2k(2n-k)$

where the arithmetic sum $a(k)$ is:

$a(k) = \sum_{i=k+1}^{n}\lambda_i$

and the geometric sum $g(k)$ is:

$g(k) = \prod_{i=k+1}^{n} \lambda_i^{-(n-k)}$

The number of relevant sinusoids in the signal subspace is determined by selecting the minimum of AIC.

`aryule`(X, order[, norm, allow_singularity])	Compute AR coefficients using Yule-Walker method
`pyule`(data, order[, norm, NFFT, sampling, ...])	Class to create PSD based on the Yule Walker method

5.3. Parametric methods ¶

5.3.1. Power Spectrum Density based on Parametric Methods ¶

5.3.1.1. ARMA and MA estimates (yule-walker)¶

5.3.1.2. AR estimate based on Burg algorithm ¶

5.3.1.3. AR estimate based on YuleWalker ¶

5.3.2. Criteria ¶

Table of Contents

Previous topic

Next topic

`arma_estimate`	Autoregressive and moving average estimators.
`ma`	Moving average estimator.
`arma2psd`	Computes power spectral density given ARMA values.
`parma`	Class to create PSD using ARMA estimator.
`pma`	Class to create PSD using MA estimator.

`arburg`(X, order[, criteria])	Estimate the complex autoregressive parameters by the Burg algorithm.
`pburg`(data, order[, criteria, NFFT, ...])	Class to create PSD based on Burg algorithm

5.3. Parametric methods¶

5.3.1. Power Spectrum Density based on Parametric Methods¶

5.3.1.1. ARMA and MA estimates (yule-walker)¶

5.3.1.2. AR estimate based on Burg algorithm¶

5.3.1.3. AR estimate based on YuleWalker¶

5.3.2. Criteria¶

5.3. Parametric methods ¶

5.3.1. Power Spectrum Density based on Parametric Methods ¶

5.3.1.2. AR estimate based on Burg algorithm ¶

5.3.1.3. AR estimate based on YuleWalker ¶

5.3.2. Criteria ¶