Skip to content

Commit

Permalink
update doctests
Browse files Browse the repository at this point in the history
  • Loading branch information
cokelaer committed Aug 19, 2017
1 parent 63fe8b5 commit aa44c68
Show file tree
Hide file tree
Showing 9 changed files with 147 additions and 124 deletions.
2 changes: 1 addition & 1 deletion doc/source/quickstart.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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)))
Expand Down
3 changes: 1 addition & 2 deletions doc/source/ref_others.rst
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,7 @@ Classes
Correlation
-------------

.. .. module:: correlation
:synopsis: Correlation methods
.. module:: correlation
.. automodule:: spectrum.correlation
:members:
:undoc-members:
Expand Down
6 changes: 5 additions & 1 deletion src/spectrum/datasets.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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()
Expand Down
151 changes: 77 additions & 74 deletions src/spectrum/linalg.py
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
"""
.. topic:: Linear algebra tools
.. topic:: Linear algebra tools
.. autosummary::
csvd
corrmtx
pascal
.. codeauthor:: Thomas Cokelaer 2011
"""
import numpy

Expand All @@ -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
"""
Expand All @@ -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:
Expand All @@ -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)
Expand All @@ -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):
Expand All @@ -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.
::
Expand Down

0 comments on commit aa44c68

Please sign in to comment.