Skip to content

Commit

Permalink
cleaning and mre tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@cokelaer
cokelaer committed Aug 26, 2017
1 parent 8388e24 commit 002216d
Show file tree
Hide file tree
Showing 9 changed files with 115 additions and 81 deletions.
2 changes: 0 additions & 2 deletions src/spectrum/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
from __future__ import absolute_import



import logging
def spectrum_set_level(level):
assert level in ['DEBUG', 'INFO', 'CRITICAL', 'ERROR', 'WARNING']
Expand Down
23 changes: 13 additions & 10 deletions src/spectrum/correlation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@
"""#from numpy.fft import fft, ifft
import numpy
import numpy as np
from numpy import arange, isrealobj
# from pylab import rms_flat

Expand All @@ -33,7 +33,7 @@ def pylab_rms_flat(a):
Return the root mean square of all the elements of *a*, flattened out.
(Copied 1:1 from matplotlib.mlab.)
"""
return numpy.sqrt(numpy.mean(numpy.absolute(a) ** 2))
return np.sqrt(np.mean(np.absolute(a) ** 2))


def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
Expand Down Expand Up @@ -87,15 +87,18 @@ def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
"""
assert norm in ['unbiased','biased', 'coeff', None]
#transform lag into list if it is an integer
x = np.array(x)
if y is None:
y = x
else:
y = np.array(y)

# N is the max of x and y
N = max(len(x), len(y))
if len(x)<N:
y = y.copy()
y.resize(N)
if len(y)<N:
if len(x) < N:
x = y.copy()
x.resize(N)
if len(y) < N:
y = y.copy()
y.resize(N)

Expand All @@ -107,9 +110,9 @@ def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
realdata = isrealobj(x) and isrealobj(y)
#create an autocorrelation array with same length as lag
if realdata == True:
r = numpy.zeros(maxlags, dtype=float)
r = np.zeros(maxlags, dtype=float)
else:
r = numpy.zeros(maxlags, dtype=complex)
r = np.zeros(maxlags, dtype=complex)

if norm == 'coeff':
rmsx = pylab_rms_flat(x)
Expand Down Expand Up @@ -143,7 +146,7 @@ def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
elif norm == 'coeff':
r[k-1] = sum/(rmsx*rmsy)/float(N)

r = numpy.insert(r, 0, r0)
r = np.insert(r, 0, r0)
return r


Expand Down Expand Up @@ -207,7 +210,7 @@ def xcorr(x, y=None, maxlags=None, norm='biased'):
assert maxlags <= N, 'maxlags must be less than data length'
lags = arange(N-maxlags-1, N+maxlags)

res = numpy.correlate(x, y, mode='full')
res = np.correlate(x, y, mode='full')

if norm == 'biased':
Nf = float(N)
Expand Down
131 changes: 66 additions & 65 deletions src/spectrum/covar.py
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,13 @@
import logging

import numpy
from .psd import ParametricSpectrum


__all__ = ["arcovar", "arcovar_marple", "pcovar", 'arcovar_marple']

def arcovar_marple(x, order, debug=False):

def arcovar_marple(x, order):
r"""Estimate AR model parameters using covariance method
This implementation is based on [Marple]_. This code is far more
Expand Down Expand Up @@ -75,29 +78,26 @@ def arcovar_marple(x, order, debug=False):
#ip +1 because we want to enter in the loop to run the first part of the code.
pbv = []
for m in range(0, order+1):
if debug:
print(('----------------------------m=', m))
print((c[0:2]))
print((d[0:2]))
logging.debug('----------------------------m=', m)
logging.debug(c[0:2])
logging.debug(d[0:2])
r1 = 1./pf
r2 = 1./pb
r3 = 1./delta
r4 = 1./gamma
if debug:
print(('starting r1r2r3r4=', r1, r2, r3, r4, pf, pb, delta, gamma))

#logging.debug('starting r1r2r3r4=', r1, r2, r3, r4, pf, pb, delta, gamma)

#Order update: AF and AB vectors ; time update: C and D vectors
temp = 0.+0.j
for k in range(m+1, N):
temp = temp + x[k]*x[k-m-1].conjugate()

if debug:
print(('temp=', temp))
r[m] = temp.conjugate()
theta = x[0] * c[m]
if debug:
print(('theta', theta))
print(('cccccccccc', c[0:2]))
print(('dddddd', d[0:2]))
#print(('theta', theta))
# print(('c=', c[0:2]))
# print(('d=', d[0:2]))
if m == 0:
pass
else:
Expand All @@ -108,31 +108,32 @@ def arcovar_marple(x, order, debug=False):
#print 'loop1 k=', k
#print ' theta=',theta, 'r[k]=',r[k], 'temp=', temp
#print ' c=',c[k], 'af=',af[m-k-1]
if m>0:
"""if m > 0:
if debug:
print((m, N-m))
print(('Xk=0',x[m-0],x[N-m-1], x[N-m+0]))
if m>1:
if m > 1:
if debug:
print('Xk=1',x[m-1],x[N-m-1], x[N-m+1])
"""
c1 = -temp * r2
c2 = -r1 * temp.conjugate()
c3 = theta * r3
c4 = r4 *theta.conjugate()
if debug:
print('c1c2c3c4 before af=',c1 ,c2 ,c3 ,c4)
#if debug:
# print('c1 c2 c3 c4 before af=',c1 ,c2 ,c3 ,c4)
af[m] = c1 # ! Eq. (8.C.19)
ab[m] = c2 # ! Eq. (8.C.22)
save = c[m]
c[m] = save + c3*d[m]
d[m] = d[m] + c4*save
if debug:
print('res',m,'af[m]=',af[m], ab[m], save, 'temp=',temp)
#if debug:
# print('res',m,'af[m]=',af[m], ab[m], save, 'temp=',temp)

if m == 0:
pass
else:
if debug:print('af before', af[0:2])
#if debug:print('af before', af[0:2])
for k in range(0, m):
save = af[k]
af[k] = save + c1 * ab[m-k-1] # Eq. (8.C.18)
Expand All @@ -141,36 +142,36 @@ def arcovar_marple(x, order, debug=False):
save = c[k]
c[k] = save + c3*d[k] # Eq. (8.C.37)
d[k] = d[k] + c4*save # Eq. (8.C.38)
if debug:
print('loop2 k=', k)
print(' af[k]=', af[k])
print(' ab[m-k-1]=', ab[m-k-1])
print(' c[k]=', c[k])
print(' d[k]=', d[k])
if debug:
if debug:
print('af after=', af[0:2])
print('ab=', ab[0:2])
#if debug:
# print('loop2 k=', k)
# print(' af[k]=', af[k])
# print(' ab[m-k-1]=', ab[m-k-1])
# print(' c[k]=', c[k])
# print(' d[k]=', d[k])
#if debug:
# print('af after=', af[0:2])
# print('ab=', ab[0:2])

r5 = temp.real**2 + temp.imag**2
pf = pf - r5*r2 # Eq. (8.C.20)
pb = pb - r5*r1 # Eq. (8.C.23)
r5 = theta.real**2 + theta.imag**2
delta = delta - r5*r4 # Eq. (8.C.39)
gamma = gamma - r5*r3 # Eq. (8.C.40)
if debug:
print('r5r2r1deltagamma', r5, r2, r1 , delta, gamma)
print('pf before norm', pf, pb, N-m-1)
#if debug:
# print('r5r2r1deltagamma', r5, r2, r1 , delta, gamma)
# print('pf before norm', pf, pb, N-m-1)
if m != order-1:
pass
else:
pf = pf / float(N-m-1)
pb = pb / float(N-m-1)
if debug:
print('ENDING', N-m-1)
#if debug:
# print('ENDING', N-m-1)
break
if debug:
print('pf and pb', pf, pb)

#if debug:
# print('pf and pb', pf, pb)
if pf > 0 and pb > 0:
pass
else:
Expand All @@ -186,14 +187,14 @@ def arcovar_marple(x, order, debug=False):
r2 = 1./pb
r3 = 1./delta
r4 = 1./gamma
if debug:
print('--------time update', r1, r2, r3, r4, m+1, N-m-1, x[m+1], x[N-m-2])
#if debug:
# print('--------time update', r1, r2, r3, r4, m+1, N-m-1, x[m+1], x[N-m-2])
ef = x[m+1]
eb = x[(N-1)-m-1]
if debug:
print('delta, gamma=', delta, gamma)
if debug:
print('before eb=', eb, ' ef=', ef)
#if debug:
# print('delta, gamma=', delta, gamma)
#if debug:
# print('before eb=', eb, ' ef=', ef)


for k in range(0,m+1):
Expand All @@ -204,28 +205,28 @@ def arcovar_marple(x, order, debug=False):

#ef = sum(af)

if debug:
print('efweb', ef , eb)
#if debug:
# print('efweb', ef , eb)
c1 = ef*r3
c2 = eb*r4
c3 = eb.conjugate() * r2
c4 = ef.conjugate() * r1
if debug:
print('c1c2c3c4', c1, c2, c3, c4)
print('af before', af[0:2])
#if debug:
# print('c1c2c3c4', c1, c2, c3, c4)
# print('af before', af[0:2])
for k in range(m, -1, -1):
save = af[k]
af[k] = save + c1 * d[k] # Eq. (8.C.33)
d[k+1] = d[k] + c4 * save # Eq. (8.C.25)
save = ab[k]
ab[k] = save + c2 * c[m-k] # Eq. (8.C.35)
c[m-k] = c[m-k] + c3 * save # Eq. (8.C.24)
if debug:
print('af after', af[0:2])
print('d', d[0:2])
print('ab', ab[0:2])
print('c', c[0:2])
if debug:print('Pb before', pf, pb)
#if debug:
# print('af after', af[0:2])
# print('d', d[0:2])
# print('ab', ab[0:2])
# print('c', c[0:2])
#if debug:print('Pb before', pf, pb)
c[m+1] = c3
d[0] = c4
#r5 = abs(ef)**2
Expand All @@ -236,24 +237,24 @@ def arcovar_marple(x, order, debug=False):
#r5 = abs(eb)**2
r5 = eb.real**2 + eb.imag**2
pb = pb - r5 * r4 # Eq. (8.C.36)
if debug:
print('Pb---------------------', m, pb, r5, r4)
#if debug:
# print('Pb---------------------', m, pb, r5, r4)
gamma = gamma-r5*r2 # Eq. (8.C.31)
if debug:
print('Gamma', gamma)
print('r5 r3,r1,and delta', r5, r3, r1, delta)
print('eb=', eb)
print('ef=', ef)
print('pf=', pf)
print('pb=', pb)
#if debug:
# print('Gamma', gamma)
# print('r5 r3,r1,and delta', r5, r3, r1, delta)
# print('eb=', eb)
# print('ef=', ef)
# print('pf=', pf)
# print('pb=', pb)
pbv.append(pb)

if (pf > 0. and pb > 0.):
pass
else:
ValueError("Negative PF or PB value")
if debug:
print(delta, gamma)
#if debug:
# print(delta, gamma)
if (delta > 0. and delta <= 1.) and (gamma > 0. and gamma <= 1.):
pass
else:
Expand Down
5 changes: 3 additions & 2 deletions src/spectrum/criteria.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,8 @@ def __init__(self, name, N):
"""
#valid attributes
self.__name = name
self.__name = None
self.name = name
self.__N = N
self.__rho = 0
self.__k = None
Expand Down Expand Up @@ -116,7 +117,7 @@ def _getK(self):
def _getN(self):
return self.__N
def _setN(self, N):
assert N>0, 'N must be positive'
assert N > 0, 'N must be positive'
self.__N = N
N = property(fget=_getN, fset=_setN, doc="Getter/Setter for N")

Expand Down
7 changes: 7 additions & 0 deletions test/test_correlation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,17 @@
from numpy.testing import assert_array_almost_equal, assert_almost_equal
from numpy import array


def test_CORRELATION():
R15 = CORRELATION(marple_data, maxlags=15, norm='biased')
R15 = CORRELATION(marple_data, maxlags=15, norm='unbiased')
R15 = CORRELATION(marple_data, maxlags=15, norm='coeff')
R15 = CORRELATION(marple_data, maxlags=15, norm=None)

def test_correlation_others():
CORRELATION(marple_data, y=array([1,2,3,4]), maxlags=5)
CORRELATION(array([1,2,3,4]), marple_data, maxlags=5)
CORRELATION(array([1,2,3,4]), marple_data)


def test_CORRELATION_biased():
Expand All @@ -23,6 +28,7 @@ def test_CORRELATION_biased():
assert_almost_equal(R30[1], R15[1])
assert_almost_equal(R30[2], R15[2])


def test_CORRELATION_unbiased():
R15 = CORRELATION(marple_data, maxlags=15, norm='unbiased')
assert_almost_equal(R15[0], 1.7804598944893049+0j)
Expand Down Expand Up @@ -75,6 +81,7 @@ def test_xcorr_versus_CORRELATION_real_data():
corr2 = CORRELATION(x, x, maxlags=4, norm=norm)
assert_array_almost_equal(corr1, func(corr2))


def _test_xcorr_versus_CORRELATION_imag_data():
from spectrum.tools import twosided as func
x = array([1,2,3,4,5+1.j])
Expand Down
7 changes: 7 additions & 0 deletions test/test_correlog.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
from spectrum import CORRELOGRAMPSD, CORRELATION, pcorrelogram, marple_data
from spectrum import data_two_freqs

from pylab import log10, plot, savefig, linspace
from numpy.testing import assert_array_almost_equal, assert_almost_equal
Expand Down Expand Up @@ -33,6 +34,12 @@ def test_pcorrelogram_class():
p = pcorrelogram(marple_data, lag=16)
p()
print(p)
p = pcorrelogram(data_two_freqs(), lag=16)
p.plot()
print(p)

def test_CORRELOGRAMPSD_others():
p = CORRELOGRAMPSD(marple_data, marple_data, lag=16, NFFT=None)

def create_figure():
psd = test_correlog()
Expand Down
Loading

0 comments on commit 002216d

Please sign in to comment.