Source code for spectrum.minvar

.. topic:: Minimum Variance Spectral Estimators

    .. autosummary::


    .. codeauthor:: Thomas Cokelaer, 2011
import numpy as np
from numpy.fft import fft

from spectrum.burg import arburg
from spectrum.psd import ParametricSpectrum
from spectrum import default_NFFT
from spectrum import errors

__all__ = ["minvar", "pminvar"]

[docs]def minvar(X, order, sampling=1., NFFT=default_NFFT): r"""Minimum Variance Spectral Estimation (MV) This function computes the minimum variance spectral estimate using the Musicus procedure. The Burg algorithm from :func:`~spectrum.burg.arburg` is used for the estimation of the autoregressive parameters. The MV spectral estimator is given by: .. math:: P_{MV}(f) = \frac{T}{e^H(f) R^{-1}_p e(f)} where :math:`R^{-1}_p` is the inverse of the estimated autocorrelation matrix (Toeplitz) and :math:`e(f)` is the complex sinusoid vector. :param X: Array of complex or real data samples (length N) :param int order: Dimension of correlation matrix (AR order = order - 1 ) :param float T: Sample interval (PSD scaling) :param int NFFT: length of the final PSD :return: * 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 :func:`minvar` The output PSD is transformed to a ``centerdc`` PSD and plotted. .. plot:: :width: 80% :include-source: 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 ) .. seealso:: * External functions used are :meth:`~spectrum.burg.arburg` and numpy.fft.fft * :class:`pminvar`, a Class dedicated to MV method. :Reference: [Marple]_ """ errors.is_positive_integer(order) errors.is_positive_integer(NFFT) psi = np.zeros(NFFT, dtype=complex) # First, we need to compute the AR values (note that order-1) A, P, k = arburg (X, order - 1) # add the order 0 A = np.insert(A, 0, 1.+0j) # We cannot compare the output with those of MARPLE in a precise way. # Indeed the burg algorithm is only single precision in fortram code # So, the AR values are slightly differnt. # The followign values are those from Marple """A[1] = 2.62284255-0.701703191j A[2] = 4.97930574-2.32781982j A[3] = 6.78445101-5.02477741j A[4] =7.85207081-8.01284409j A[5] =7.39412165-10.7684202j A[6] =6.03175116-12.7067814j A[7] =3.80106878-13.6808891j A[8] =1.48207295-13.2265558j A[9] =-0.644280195-11.4574194j A[10] =-2.02386642-8.53268814j A[11] =-2.32437634-5.25636244j A[12] =-1.75356281-2.46820402j A[13] =-0.888899028-0.781434655j A[14] =-0.287197977-0.0918145925j P = 0.00636525545 """ # if we use exactly the same AR coeff and P from Marple Burg output, then # we can compare the following code. This has been done and reveals that # the FFT in marple is also slightly different (precision) from this one. # However, the results are sufficiently close (when NFFT is small) that # we are confident the following code is correct. # Compute the psi coefficients for K in range(0, order): SUM = 0. MK = order-K # Correlate the autoregressive parameters for I in range(0, order - K): SUM = SUM + float(MK-2*I) * A[I].conjugate()*A[I+K] # Eq. (12.25) SUM = SUM/P if K != 0: psi[NFFT-K] = SUM.conjugate() psi[K] = SUM # Compute FFT of denominator psi = fft(psi, NFFT) # Invert the psi terms at this point to get PSD values PSD = sampling / np.real(psi) return PSD, A, k
[docs]class pminvar(ParametricSpectrum): """Class to create PSD based on the Minimum variance spectral estimation See :func:`minvar` for description. .. plot:: :width: 80% :include-source: from spectrum import * p = pminvar(marple_data, 15, NFFT=4096) p.plot(sides='centerdc') """ def __init__(self, data, order, NFFT=None, sampling=1., scale_by_freq=False): """**Constructor** For a detailled description of the parameters, see :func:`minvar`. :param array data: input data (list or numpy.array) :param int order: :param int NFFT: total length of the final data sets (padded with zero if needed; default is 4096) :param float sampling: sampling frequency of the input :attr:`data`. """ super(pminvar, self).__init__(data, ar_order=order, sampling=sampling, NFFT=NFFT, scale_by_freq=scale_by_freq) def __call__(self): res = minvar(, self.ar_order, sampling=self.sampling, NFFT=self.NFFT) psd = res[0] # save the AR and reflection coefficients. = res[1] self.reflection = res[2] # save the PSD if self.datatype == 'real': if self.NFFT % 2 == 0: newpsd = psd[0:int(self.NFFT/2+1)] * 2 else: newpsd = psd[0:int((self.NFFT+1)/2)] * 2 self.psd = newpsd else: self.psd = psd self.scale() def _str_title(self): return "Minimum Variance spectral estimation\n" def __str__(self): return super(pminvar, self).__str__()