From aa44c68b3bb2df0791eeeef1c8d97420a4cc7222 Mon Sep 17 00:00:00 2001 From: Thomas Cokelaer Date: Sat, 19 Aug 2017 23:06:30 +0200 Subject: [PATCH] update doctests --- doc/source/quickstart.rst | 2 +- doc/source/ref_others.rst | 3 +- src/spectrum/datasets.py | 6 +- src/spectrum/linalg.py | 151 +++++++++++++++--------------- src/spectrum/linear_prediction.py | 93 ++++++++++-------- src/spectrum/mtm.py | 2 +- src/spectrum/tools.py | 9 +- src/spectrum/window.py | 3 +- src/spectrum/yulewalker.py | 2 +- 9 files changed, 147 insertions(+), 124 deletions(-) diff --git a/doc/source/quickstart.rst b/doc/source/quickstart.rst index 2c475d3..8a5af59 100644 --- a/doc/source/quickstart.rst +++ b/doc/source/quickstart.rst @@ -132,7 +132,7 @@ Once the AR and/or MA parameters are found, the :func:`~spectrum.arma.arma2psd` :width: 80% from spectrum import arma_estimate, arma2psd, marple_data - from pylab import * + from pylab import plot, axis, xlabel, ylabel, grid ar, ma, rho = arma_estimate(marple_data, 15, 15, 30) psd = arma2psd(ar, ma, rho=rho, NFFT=4096) plot(10*log10(psd/max(psd))) diff --git a/doc/source/ref_others.rst b/doc/source/ref_others.rst index aa1e880..adfb9a4 100644 --- a/doc/source/ref_others.rst +++ b/doc/source/ref_others.rst @@ -19,8 +19,7 @@ Classes Correlation ------------- -.. .. module:: correlation - :synopsis: Correlation methods +.. module:: correlation .. automodule:: spectrum.correlation :members: :undoc-members: diff --git a/src/spectrum/datasets.py b/src/spectrum/datasets.py index f46eda2..c0ffa9e 100644 --- a/src/spectrum/datasets.py +++ b/src/spectrum/datasets.py @@ -110,7 +110,10 @@ def data_cosine(N=1024, A=0.1, sampling=1024., freq=200): where w[t] is a white noise of variance 1. - >>> a = data_cosine(N=1024, sampling=1024, A=0.5, freq=100) + .. doctest:: + + >>> from spectrum import data_cosine + >>> a = data_cosine(N=1024, sampling=1024, A=0.5, freq=100) """ t = arange(0, float(N)/sampling, 1./sampling) @@ -132,6 +135,7 @@ class TimeSeries(): .. doctest:: + >>> from spectrum import TimeSeries >>> data = [1, 2, 3, 4, 3, 2, 1, 0 ] >>> ts = TimeSeries(data, sampling=1) >>> ts.plot() diff --git a/src/spectrum/linalg.py b/src/spectrum/linalg.py index aceeed7..f42d556 100644 --- a/src/spectrum/linalg.py +++ b/src/spectrum/linalg.py @@ -1,15 +1,15 @@ """ -.. topic:: Linear algebra tools +.. topic:: Linear algebra tools .. autosummary:: - + csvd corrmtx pascal - + .. codeauthor:: Thomas Cokelaer 2011 - - + + """ import numpy @@ -20,16 +20,19 @@ def pascal(n): """Return Pascal matrix - + :param int n: size of the matrix - - >>> pascal(6) - array([[ 1., 1., 1., 1., 1., 1.], - [ 1., 2., 3., 4., 5., 6.], - [ 1., 3., 6., 10., 15., 21.], - [ 1., 4., 10., 20., 35., 56.], - [ 1., 5., 15., 35., 70., 126.], - [ 1., 6., 21., 56., 126., 252.]]) + + .. doctest:: + + >>> from spectrum import pascal + >>> pascal(6) + array([[ 1., 1., 1., 1., 1., 1.], + [ 1., 2., 3., 4., 5., 6.], + [ 1., 3., 6., 10., 15., 21.], + [ 1., 4., 10., 20., 35., 56.], + [ 1., 5., 15., 35., 70., 126.], + [ 1., 6., 21., 56., 126., 252.]]) .. todo:: use the symmetric property to improve computational time if needed """ @@ -44,108 +47,108 @@ def pascal(n): for i in range(1, n): for j in range(1, n): result[i, j] = result[i-1, j] + result[i, j-1] - return result + return result def corrmtx(x_input, m, method='autocorrelation'): r"""Correlation matrix - + This function is used by PSD estimator functions. It generates the correlation matrix from a correlation data set and a maximum lag. - + :param array x: autocorrelation samples (1D) :param int m: the maximum lag - - Depending on the choice of the method, the correlation matrix has different - sizes, but the number of rows is always m+1. + + Depending on the choice of the method, the correlation matrix has different + sizes, but the number of rows is always m+1. Method can be : - - * 'autocorrelation': (default) X is the (n+m)-by-(m+1) rectangular Toeplitz + + * 'autocorrelation': (default) X is the (n+m)-by-(m+1) rectangular Toeplitz matrix derived using prewindowed and postwindowed data. * 'prewindowed': X is the n-by-(m+1) rectangular Toeplitz matrix derived using prewindowed data only. - * 'postwindowed': X is the n-by-(m+1) rectangular Toeplitz matrix that + * 'postwindowed': X is the n-by-(m+1) rectangular Toeplitz matrix that derived using postwindowed data only. - * 'covariance': X is the (n-m)-by-(m+1) rectangular Toeplitz matrix + * 'covariance': X is the (n-m)-by-(m+1) rectangular Toeplitz matrix derived using nonwindowed data. * 'modified': X is the 2(n-m)-by-(m+1) modified rectangular Toeplitz - matrix that generates an autocorrelation estimate for the length n data + matrix that generates an autocorrelation estimate for the length n data vector x, derived using forward and backward prediction error estimates. - - :return: + + :return: * the autocorrelation matrix * R, the (m+1)-by-(m+1) autocorrelation matrix estimate ``R= X'*X``. - + .. rubric:: Algorithm details: - + The **autocorrelation** matrix is a :math:`(N+p) \times (p+1)` rectangular Toeplilz data matrix: - + .. math:: X_p = \begin{pmatrix}L_p\\T_p\\Up\end{pmatrix} - - where the lower triangular :math:`p \times (p+1)` matrix :math:`L_p` is - - .. math:: L_p = + + where the lower triangular :math:`p \times (p+1)` matrix :math:`L_p` is + + .. math:: L_p = \begin{pmatrix} x[1] & \cdots & 0 & 0 \\ \vdots & \ddots & \vdots & \vdots \\ - x[p] & \cdots & x[1] & 0 - \end{pmatrix} - + x[p] & \cdots & x[1] & 0 + \end{pmatrix} + where the rectangular :math:`(N-p) \times (p+1)` matrix :math:`T_p` is - - .. math:: T_p = + + .. math:: T_p = \begin{pmatrix} x[p+1] & \cdots & x[1] \\ \vdots & \ddots & \vdots \\ x[N-p] & \cdots & x[p+1] \\ \vdots & \ddots & \vdots \\ - x[N] & \cdots & x[N-p] + x[N] & \cdots & x[N-p] \end{pmatrix} - + and where the upper triangular :math:`p \times (p+1)` matrix :math:`U_p` is - - .. math:: U_p = + + .. math:: U_p = \begin{pmatrix} 0 & x[N] & \cdots & x[N-p+1] \\ \vdots & \vdots & \ddots & \vdots \\ - 0 & 0 & \cdots & x[N] + 0 & 0 & \cdots & x[N] \end{pmatrix} - From this definition, the prewindowed matrix is - + From this definition, the prewindowed matrix is + .. math:: X_p = \begin{pmatrix}L_p\\T_p\end{pmatrix} - - the postwindowed matrix is - + + the postwindowed matrix is + .. math:: X_p = \begin{pmatrix}T_p\\U_p\end{pmatrix} - + the covariance matrix is: - + .. math:: X_p = \begin{pmatrix}T_p\end{pmatrix} - + and the modified covariance matrix is: - + .. math:: X_p = \begin{pmatrix}T_p\\T_p^*\end{pmatrix} - - + + """ - - valid_methods = ['autocorrelation', 'prewindowed', 'postwindowed', + + valid_methods = ['autocorrelation', 'prewindowed', 'postwindowed', 'covariance', 'modified'] if method not in valid_methods: raise ValueError("Invalid method. Try one of %s" % valid_methods) from scipy.linalg import toeplitz - + # create the relevant matrices that will be useful to create # the correlation matrices N = len(x_input) - # FIXME:do we need a copy ? + # FIXME:do we need a copy ? if isinstance(x_input, list): x = numpy.array(x_input) else: @@ -156,18 +159,18 @@ def corrmtx(x_input, m, method='autocorrelation'): complex_type = True else: complex_type = False - + # Compute the Lp, Up and Tp matrices according to the requested method if method in ['autocorrelation', 'prewindowed']: Lp = toeplitz(x[0:m], [0]*(m+1)) Tp = toeplitz(x[m:N], x[m::-1]) if method in ['autocorrelation', 'postwindowed']: Up = toeplitz([0]*(m+1), numpy.insert(x[N:N-m-1:-1],0,0)) - - - + + + # Create the output matrix - + if method == 'autocorrelation': if complex_type == True: C = numpy.zeros((N+m, m+1), dtype=complex) @@ -184,7 +187,7 @@ def corrmtx(x_input, m, method='autocorrelation'): C = numpy.zeros((N, m+1), dtype=complex) else: C = numpy.zeros((N, m+1)) - + for i in range(0, m): C[i] = Lp[i] for i in range(m, N): @@ -210,28 +213,28 @@ def corrmtx(x_input, m, method='autocorrelation'): Tp = numpy.fliplr(Tp.conj()) for i in range(N-m, 2*(N-m)): C[i] = Tp[i-N+m] - - - + + + return C #autocorrelation n+m by m+1 matrix (Lp, Tp, Up) - + #Rx = toeplitz(r) - + def csvd(A): """SVD decomposition using numpy.linalg.svd :param A: a M by N matrix - - :return: + + :return: * U, a M by M matrix * S the N eigen values * V a N by N matrix - + See :func:`numpy.linalg.svd` for a detailed documentation. - + Should return the same as in [Marple]_ , CSVD routine. :: diff --git a/src/spectrum/linear_prediction.py b/src/spectrum/linear_prediction.py index 1068655..58d9539 100644 --- a/src/spectrum/linear_prediction.py +++ b/src/spectrum/linear_prediction.py @@ -7,7 +7,7 @@ import numpy from scipy.signal import deconvolve -__all__ = ['ac2poly', 'poly2ac', 'ac2rc', 'rc2poly', 'rc2ac', 'is2rc', 'rc2is', +__all__ = ['ac2poly', 'poly2ac', 'ac2rc', 'rc2poly', 'rc2ac', 'is2rc', 'rc2is', 'rc2lar', 'lar2rc', 'poly2rc', 'lsf2poly', 'poly2lsf'] """ @@ -20,7 +20,7 @@ def ac2poly(data): """Convert autocorrelation sequence to prediction polynomial :param array data: input data (list or numpy.array) - :return: + :return: * AR parameters * noise variance @@ -32,16 +32,18 @@ def ac2poly(data): .. doctest:: - >>> r = [5, -2, 1] + >>> from spectrum import ac2poly + >>> from numpy import array + >>> r = [5, -2, 1.01] >>> ar, e = ac2poly(r) >>> ar - array([ 1. , 0. , -0.2]) + array([ 1. , 0.38, -0.05]) >>> e - 4.7999999999999998 + 4.1895000000 """ a, e, _c = LEVINSON(data) - a = numpy.insert(a, 0, 1) + a = numpy.insert(a, 0, 1) return a, e @@ -59,39 +61,44 @@ def ac2rc(data): """ a, e, _c = LEVINSON(data) return a, data[0] - - + + def poly2ac(poly, efinal): """ Convert prediction filter polynomial to autocorrelation sequence :param array poly: the AR parameters :param efinal: an estimate of the final error - :return: the autocorrelation sequence + :return: the autocorrelation sequence in complex format. .. doctest:: - >>> ar = [ 1. , 0. , -0.2] - >>> efinal = 4.8 - >>> r = poly2ac(ar, efinal) + >>> from numpy import array + >>> from spectrum import poly2ac + >>> poly = [ 1. , 0.38 , -0.05] + >>> efinal = 4.1895 + >>> poly2ac(poly, efinal) + array([ 5.00+0.j, -2.00+0.j, 1.01-0.j]) """ results = rlevinson(poly, efinal) - return results[0] - + return results[0] + + def ar2rc(ar): """ Convert autoregressive parameters into reflection coefficients """ raise NotImplementedError - - + + def poly2rc(a, efinal): """Convert prediction filter polynomial to reflection coefficients :param a: AR parameters - :param efinal: + :param efinal: """ results = rlevinson(a, efinal) return results[2] - + + def rc2poly(kr, r0=None): """convert reflection coefficients to prediction filter polynomial @@ -114,7 +121,7 @@ def rc2poly(kr, r0=None): e[0] = e0 * (1. - numpy.conj(numpy.conjugate(kr[0])*kr[0])) - # Continue the recursion for k=2,3,...,p, where p is the order of the + # Continue the recursion for k=2,3,...,p, where p is the order of the # prediction polynomial. for k in range(1, p): @@ -144,17 +151,17 @@ def is2rc(inv_sin): :param inv_sin: inverse sine parameters :return: reflection coefficients - + .. seealso:: :func:`rc2is`, :func:`poly2rc`, :func:`ac2rc`, :func:`lar2rc`. - :Reference: J.R. Deller, J.G. Proakis, J.H.L. Hansen, "Discrete-Time Processing of Speech Signals", Prentice Hall, Section 7.4.5. + :Reference: J.R. Deller, J.G. Proakis, J.H.L. Hansen, + "Discrete-Time Processing of Speech Signals", Prentice Hall, Section 7.4.5. """ return numpy.sin(numpy.array(inv_sin)*numpy.pi/2); - -def rc2is(k): +def rc2is(k): """Convert reflection coefficients to inverse sine parameters. :param k: reflection coefficients @@ -162,7 +169,7 @@ def rc2is(k): .. seealso:: :func:`is2rc`, :func:`rc2poly`, :func:`rc2acC`, :func:`rc2lar`. - Reference: J.R. Deller, J.G. Proakis, J.H.L. Hansen, "Discrete-Time + Reference: J.R. Deller, J.G. Proakis, J.H.L. Hansen, "Discrete-Time Processing of Speech Signals", Prentice Hall, Section 7.4.5. """ @@ -178,7 +185,8 @@ def rc2lar(k): :param k: reflection coefficients :return: inverse sine parameters - The log area ratio is defined by G = log((1+k)/(1-k)) , where the K is the reflection coefficient. + The log area ratio is defined by G = log((1+k)/(1-k)) , where the K + parameter is the reflection coefficient. .. seealso:: :func:`lar2rc`, :func:`rc2poly`, :func:`rc2ac`, :func:`rc2ic`. @@ -220,15 +228,17 @@ def lsf2poly(lsf): .. doctest:: + >>> from spectrum import lsf2poly >>> lsf = [0.7842 , 1.5605 , 1.8776 , 1.8984, 2.3593] >>> a = lsf2poly(lsf) - array([ 1.00000000e+00, 6.14837835e-01, 9.89884967e-01, - 9.31594056e-05, 3.13713832e-03, -8.12002261e-03 ]) + + # array([ 1.00000000e+00, 6.14837835e-01, 9.89884967e-01, + # 9.31594056e-05, 3.13713832e-03, -8.12002261e-03 ]) .. seealso:: poly2lsf, rc2poly, ac2poly, rc2is """ # Reference: A.M. Kondoz, "Digital Speech: Coding for Low Bit Rate Communications - # Systems" John Wiley & Sons 1994 ,Chapter 4 + # Systems" John Wiley & Sons 1994 ,Chapter 4 # Line spectral frequencies must be real. @@ -245,18 +255,18 @@ def lsf2poly(lsf): # Separate the zeros to those belonging to P and Q rQ = z[0::2] rP = z[1::2] - + # Include the conjugates as well rQ = numpy.concatenate((rQ, rQ.conjugate())) rP = numpy.concatenate((rP, rP.conjugate())) - + # Form the polynomials P and Q, note that these should be real Q = numpy.poly(rQ); P = numpy.poly(rP); - + # Form the sum and difference filters by including known roots at z = 1 and - # z = -1 - + # z = -1 + if p%2: # Odd order: z = +1 and z = -1 are roots of the difference filter, P1(z) P1 = numpy.convolve(P, [1, 0, -1]) @@ -270,21 +280,22 @@ def lsf2poly(lsf): # Prediction polynomial is formed by averaging P1 and Q1 a = .5 * (P1+Q1) - return a[0:-1:1] # do not return last element + return a[0:-1:1] # do not return last element def poly2lsf(a): """Prediction polynomial to line spectral frequencies. converts the prediction polynomial specified by A, - into the corresponding line spectral frequencies, LSF. + into the corresponding line spectral frequencies, LSF. normalizes the prediction polynomial by A(1). .. doctest:: - >>> a = [1.0000 , 0.6149 , 0.9899 , 0.0000 , 0.0031, -0.0082 + >>> from spectrum import poly2lsf + >>> a = [1.0000, 0.6149, 0.9899, 0.0000 ,0.0031, -0.0082] >>> lsf = poly2lsf(a) - >>> lsf = array([ 0.7842, 1.5605 , 1.8776 , 1.8984, 2.3593]) + >>> lsf = array([0.7842, 1.5605, 1.8776, 1.8984, 2.3593]) .. seealso:: lsf2poly, poly2rc, poly2qc, rc2is """ @@ -304,10 +315,10 @@ def poly2lsf(a): # Form the sum and differnce filters p = len(a)-1 # The leading one in the polynomial is not used - a1 = numpy.concatenate((a, numpy.array([0]))) + a1 = numpy.concatenate((a, numpy.array([0]))) a2 = a1[-1::-1] P1 = a1 - a2 # Difference filter - Q1 = a1 + a2 # Sum Filter + Q1 = a1 + a2 # Sum Filter # If order is even, remove the known root at z = 1 for P1 and z = -1 for Q1 # If odd, remove both the roots from P1 @@ -315,10 +326,10 @@ def poly2lsf(a): if p%2: # Odd order P, r = deconvolve(P1,[1, 0 ,-1]) Q = Q1 - else: # Even order + else: # Even order P, r = deconvolve(P1, [1, -1]) Q, r = deconvolve(Q1, [1, 1]) - + rP = numpy.roots(P) rQ = numpy.roots(Q) diff --git a/src/spectrum/mtm.py b/src/spectrum/mtm.py index a0509a2..81aa2e0 100644 --- a/src/spectrum/mtm.py +++ b/src/spectrum/mtm.py @@ -87,7 +87,7 @@ def pmtm(x, NW=None, k=None, NFFT=None, e=None, v=None, method='adapt', show=Fal :width: 80% :include-source: - from spectrum import * + from spectrum import data_as, dpss, pmtm data = data_cosine(N=2048, A=0.1, sampling=1024, freq=200) # If you already have the DPSS windows diff --git a/src/spectrum/tools.py b/src/spectrum/tools.py index 4886c98..c499da2 100644 --- a/src/spectrum/tools.py +++ b/src/spectrum/tools.py @@ -73,6 +73,7 @@ def twosided_2_onesided(data): psd[-1] = data[-1] return psd + def onesided_2_twosided(data): """Convert a two-sided PSD to a one-sided PSD @@ -92,6 +93,7 @@ def onesided_2_twosided(data): psd[-1] *= 2. return psd + def twosided_2_centerdc(data): """Convert a two-sided PSD to a center-dc PSD""" N = len(data) @@ -100,6 +102,7 @@ def twosided_2_centerdc(data): newpsd[0] = data[-1] return newpsd + def centerdc_2_twosided(data): """Convert a center-dc PSD to a twosided PSD""" N = len(data) @@ -169,6 +172,7 @@ def pow2db(x): .. doctest:: + >>> from spectrum import pow2db >>> x = pow2db(0.1) >>> x -10.0 @@ -181,6 +185,7 @@ def db2pow(xdb): .. doctest:: + >>> from spectrum import db2pow >>> p = db2pow(-10) >>> p 0.1 @@ -201,6 +206,7 @@ def nextpow2(x): .. doctest:: + >>> from spectrum import nextpow2 >>> x = [255, 256, 257] >>> nextpow2(x) array([8, 8, 9]) @@ -215,6 +221,7 @@ def db2mag(xdb): .. doctest:: + >>> from spectrum import db2mag >>> db2mag(-20) 0.1 @@ -232,11 +239,11 @@ def mag2db(x): .. doctest:: + >>> from spectrum import mag2db >>> mag2db(0.1) -20.0 .. seealso:: :func:`db2mag` """ - return 20. * numpy.log10(x) diff --git a/src/spectrum/window.py b/src/spectrum/window.py index 6a02858..3cc0fdd 100644 --- a/src/spectrum/window.py +++ b/src/spectrum/window.py @@ -454,9 +454,8 @@ def enbw(data): .. doctest:: - >>> from spectrum import create_window + >>> from spectrum import create_window, enbw >>> w = create_window(64, 'rectangular') - >>> enbw(w) 1.0 diff --git a/src/spectrum/yulewalker.py b/src/spectrum/yulewalker.py index 9c93eb3..c320663 100644 --- a/src/spectrum/yulewalker.py +++ b/src/spectrum/yulewalker.py @@ -14,7 +14,7 @@ from .correlation import CORRELATION from .levinson import LEVINSON from .psd import ParametricSpectrum - +from spectrum import tools __all__ = ['aryule', 'pyule']