[Scipy-svn] r2332 - in trunk/Lib/sandbox: . cdavid cdavid/src cdavid/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Nov 28 03:05:13 EST 2006
Author: cdavid
Date: 2006-11-28 02:04:47 -0600 (Tue, 28 Nov 2006)
New Revision: 2332
Added:
trunk/Lib/sandbox/cdavid/
trunk/Lib/sandbox/cdavid/Changelog
trunk/Lib/sandbox/cdavid/TODO
trunk/Lib/sandbox/cdavid/__init__.py
trunk/Lib/sandbox/cdavid/autocorr.py
trunk/Lib/sandbox/cdavid/info.py
trunk/Lib/sandbox/cdavid/lpc.py
trunk/Lib/sandbox/cdavid/setup.py
trunk/Lib/sandbox/cdavid/src/
trunk/Lib/sandbox/cdavid/src/Makefile
trunk/Lib/sandbox/cdavid/src/README
trunk/Lib/sandbox/cdavid/src/autocorr_nofft.c
trunk/Lib/sandbox/cdavid/src/autocorr_nofft.def
trunk/Lib/sandbox/cdavid/src/autocorr_nofft.h
trunk/Lib/sandbox/cdavid/src/autocorr_nofft.tpl
trunk/Lib/sandbox/cdavid/src/common.h
trunk/Lib/sandbox/cdavid/src/levinson.c
trunk/Lib/sandbox/cdavid/src/levinson.def
trunk/Lib/sandbox/cdavid/src/levinson.h
trunk/Lib/sandbox/cdavid/src/levinson.tpl
trunk/Lib/sandbox/cdavid/src/lpc.c
trunk/Lib/sandbox/cdavid/src/lpc.def
trunk/Lib/sandbox/cdavid/src/lpc.h
trunk/Lib/sandbox/cdavid/src/lpc.tpl
trunk/Lib/sandbox/cdavid/tests/
trunk/Lib/sandbox/cdavid/tests/test_autocorr.py
trunk/Lib/sandbox/cdavid/tests/test_lpc.py
Modified:
trunk/Lib/sandbox/setup.py
Log:
Intial commit of cdavid, for lpc
Added: trunk/Lib/sandbox/cdavid/Changelog
===================================================================
--- trunk/Lib/sandbox/cdavid/Changelog 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/Changelog 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,5 @@
+pyem (0.1) Tue, 28 Nov 2006 16:56:35 +0900
+
+ * first release
+
+-- David Cournapeau <david at ar.media.kyoto-u.ac.jp>
Added: trunk/Lib/sandbox/cdavid/TODO
===================================================================
--- trunk/Lib/sandbox/cdavid/TODO 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/TODO 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,6 @@
+# Last Change: Tue Nov 28 05:00 PM 2006 J
+
+ - there is no doc.
+ - the handling of non contiguous arrays is not really
+ elegant, and the code is difficult to maintain
+ - rank > 2: must code in C ? (yuk)
Added: trunk/Lib/sandbox/cdavid/__init__.py
===================================================================
--- trunk/Lib/sandbox/cdavid/__init__.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/__init__.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,9 @@
+# Last Change: Tue Nov 28 04:00 PM 2006 J
+from info import __doc__
+
+from lpc import lpc2 as lpc
+from autocorr import autocorr_oneside_nofft as autocorr
+
+from numpy.testing import NumpyTest
+test = NumpyTest().test
+
Added: trunk/Lib/sandbox/cdavid/autocorr.py
===================================================================
--- trunk/Lib/sandbox/cdavid/autocorr.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/autocorr.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,322 @@
+#! /usr/bin/env python
+# Last Change: Tue Nov 28 03:00 PM 2006 J
+
+# TODO: - proper test
+# TODO: - proper profiling
+
+from numpy.fft import fft, ifft
+from numpy import correlate, log2, floor, conj, real, \
+ concatenate, sum, max
+
+from warnings import warn
+
+# use ctype to have one sided c imp of autocorr
+import ctypes
+from ctypes import c_uint, c_int
+from numpy.ctypeslib import ndpointer, load_library
+
+ctypes_major = int(ctypes.__version__.split('.')[0])
+if ctypes_major < 1:
+ msg = "version of ctypes is %s, expected at least %s" \
+ % (ctypes.__version__, '1.0.1')
+ raise importerror(msg)
+
+import numpy as N
+
+# load autocorr lib
+_autocorr = load_library('gabsig.so', __file__)
+
+#===============================
+# define the functions with args
+#===============================
+
+# contiguous 1d
+arg1 = ndpointer(dtype = N.float64, flags='CONTIGUOUS,ALIGNED')
+arg2 = c_uint
+arg3 = ndpointer(dtype = N.float64, flags='CONTIGUOUS,ALIGNED')
+arg4 = c_uint
+_autocorr.dbl_xcorr_nofft_1d.argtypes = [arg1, arg2, arg3, arg4]
+_autocorr.dbl_xcorr_nofft_1d.restype = c_int
+
+arg1 = ndpointer(dtype = N.float32, flags='CONTIGUOUS,ALIGNED')
+arg2 = c_uint
+arg3 = ndpointer(dtype = N.float32, flags='CONTIGUOUS,ALIGNED')
+arg4 = c_uint
+_autocorr.flt_xcorr_nofft_1d.argtypes = [arg1, arg2, arg3, arg4]
+_autocorr.flt_xcorr_nofft_1d.restype = c_int
+
+# non contiguous 1d
+arg1 = ndpointer(dtype = N.float64, flags = 'ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = N.float64, flags = 'ALIGNED')
+arg5 = c_uint
+arg6 = c_uint
+_autocorr.dbl_xcorr_nofft_1d_noncontiguous.argtypes = [arg1, \
+ arg2, arg3, arg4, arg5, arg6]
+_autocorr.dbl_xcorr_nofft_1d_noncontiguous.restype = c_int
+
+arg1 = ndpointer(dtype = N.float32, flags = 'ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = N.float32, flags = 'ALIGNED')
+arg5 = c_uint
+arg6 = c_uint
+_autocorr.flt_xcorr_nofft_1d_noncontiguous.argtypes = [arg1, \
+ arg2, arg3, arg4, arg5, arg6]
+_autocorr.flt_xcorr_nofft_1d_noncontiguous.restype = c_int
+
+# contiguous 2d
+arg1 = ndpointer(dtype = N.float64, flags='ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = N.float64, flags='ALIGNED')
+arg5 = c_uint
+_autocorr.dbl_xcorr_nofft_2d.argtypes = [arg1, arg2, arg3, arg4, arg5]
+_autocorr.dbl_xcorr_nofft_2d.restype = c_int
+
+arg1 = ndpointer(dtype = N.float32, flags='ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = N.float32, flags='ALIGNED')
+arg5 = c_uint
+_autocorr.flt_xcorr_nofft_2d.argtypes = [arg1, arg2, arg3, arg4, arg5]
+_autocorr.flt_xcorr_nofft_2d.restype = c_int
+
+# non contiguous 2d
+arg1 = ndpointer(dtype = N.float64, flags='ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = c_uint
+arg5 = c_uint
+arg6 = ndpointer(dtype = N.float64, flags='ALIGNED')
+arg7 = c_uint
+arg8 = c_uint
+arg9 = c_uint
+_autocorr.dbl_xcorr_nofft_2d_noncontiguous.argtypes = [arg1, arg2, \
+ arg3, arg4, arg5, arg6, arg7, arg8, arg9]
+_autocorr.dbl_xcorr_nofft_2d_noncontiguous.restype = c_int
+
+arg1 = ndpointer(dtype = N.float32, flags='ALIGNED')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = c_uint
+arg5 = c_uint
+arg6 = ndpointer(dtype = N.float32, flags='ALIGNED')
+arg7 = c_uint
+arg8 = c_uint
+arg9 = c_uint
+_autocorr.flt_xcorr_nofft_2d_noncontiguous.argtypes = [arg1, arg2, \
+ arg3, arg4, arg5, arg6, arg7, arg8, arg9]
+_autocorr.flt_xcorr_nofft_2d_noncontiguous.restype = c_int
+
+#======================================
+# Fonctions to be used for testing only
+#======================================
+def _raw_autocorr_1d(signal, lag):
+ assert signal.ndim == 1
+ assert signal.flags['CONTIGUOUS']
+
+ if lag >= signal.size:
+ raise RuntimeError("lag should be < to input size")
+
+ if signal.dtype == N.float64:
+ res = N.zeros((lag+1), N.float64)
+ _autocorr.dbl_xcorr_nofft_1d(signal, signal.size, res, lag)
+ elif signal.dtype == N.float32:
+ res = N.zeros((lag+1), N.float32)
+ _autocorr.flt_xcorr_nofft_1d(signal, signal.size, res, lag)
+ else:
+ raise TypeError("only float 32 and 64 bits supported for now")
+
+ return res
+
+def _raw_autocorr_1d_noncontiguous(signal, lag):
+ assert signal.ndim == 1
+
+ if lag >= signal.size:
+ raise RuntimeError("lag should be < to input size")
+
+ if signal.dtype == N.float64:
+ res = N.zeros((lag+1), N.float64)
+ _autocorr.dbl_xcorr_nofft_1d_noncontiguous(signal, signal.size,
+ signal.strides[0], res, res.strides[0], lag)
+ elif signal.dtype == N.float32:
+ res = N.zeros((lag+1), N.float32)
+ _autocorr.flt_xcorr_nofft_1d_noncontiguous(signal, signal.size,
+ signal.strides[0], res, res.strides[0], lag)
+ else:
+ raise TypeError("only float 32 and 64 bits supported for now")
+
+ return res
+
+# python implementation of autocorr for rank <= 2
+def _autocorr_oneside_nofft_py(signal, lag, axis = -1):
+ if signal.ndim > 2:
+ raise NotImplemented("only for rank <=2")
+
+ if axis % 2 == 0:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype)
+ center = signal.shape[0] - 1
+ for i in range(signal.shape[1]):
+ #print "compute corr of " + str(signal[:, i])
+ res[:, i] = correlate(signal[:, i], signal[:, i], \
+ 'full')[center:center+lag+1]
+ elif axis % 2 == 1:
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ center = signal.shape[1] - 1
+ for i in range(signal.shape[0]):
+ #print "compute corr of " + str(signal[i])
+ res[i] = correlate(signal[i], signal[i], \
+ 'full')[center:center+lag+1]
+ else:
+ raise RuntimeError("this should bnot happen, please fill a bug")
+
+ return res
+
+#=============
+# Public API
+#=============
+def autocorr_oneside_nofft(signal, lag, axis = -1):
+ """Compute the righ side of autocorrelation along the axis, for lags up to lag.
+
+ This implementation does NOT use FFT."""
+ # TODO For rank < 2, the overhead of python code may be significant. Should
+ # TODO not be difficult to do in C anyway (we can still use ctypes)
+
+ # rank 0, 1
+ if signal.ndim < 2:
+ size = signal.shape[-1]
+ if lag >= size:
+ raise RuntimeError("lag should be < to input size")
+
+ res = N.zeros((lag+1), signal.dtype)
+
+ if signal.flags['CONTIGUOUS']:
+ if signal.dtype == N.float64:
+ _autocorr.dbl_xcorr_nofft_1d(signal, size, res, lag)
+ elif signal.dtype == N.float32:
+ _autocorr.flt_xcorr_nofft_1d(signal, size, res, lag)
+ else:
+ raise TypeError("only float 32 and 64 bits supported for now")
+ else:
+ istride = signal.strides[0]
+ ostride = signal.itemsize
+ if signal.dtype == N.float64:
+ _autocorr.dbl_xcorr_nofft_1d_noncontiguous(signal, size, istride,
+ res, ostride, lag)
+ elif signal.dtype == N.float32:
+ _autocorr.flt_xcorr_nofft_1d_noncontiguous(signal, size, istride,
+ res, ostride, lag)
+ else:
+ raise TypeError("only float 32 and 64 bits supported for now")
+
+ # rank 2 case
+ elif signal.ndim == 2:
+ size = signal.shape[axis]
+ if lag >= size:
+ raise RuntimeError("lag should be < to input size")
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ else:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype)
+
+ if signal.dtype == N.float64:
+ # contiguous case
+ if signal.flags['C'] and axis % 2 == 1:
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ _autocorr.dbl_xcorr_nofft_2d(signal, signal.shape[0], signal.shape[1],
+ res, lag)
+ # contiguous case
+ elif signal.flags['F'] and axis % 2 == 0:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype, order = 'F')
+ _autocorr.dbl_xcorr_nofft_2d(signal, signal.shape[1], signal.shape[0],
+ res, lag)
+ # non contiguous case
+ elif axis % 2 == 0:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype)
+ warn("non contiguous used, this will be slow")
+ _autocorr.dbl_xcorr_nofft_2d_noncontiguous(signal,
+ signal.shape[1], signal.shape[0],
+ signal.strides[1], signal.strides[0],
+ res, res.strides[1], res.strides[0], lag)
+ elif axis % 2 == 1:
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ warn("non contiguous used, this will be slow")
+ _autocorr.dbl_xcorr_nofft_2d_noncontiguous(signal,
+ signal.shape[0], signal.shape[1],
+ signal.strides[0], signal.strides[1],
+ res, res.strides[0], res.strides[1], lag)
+ elif signal.dtype == N.float32:
+ # contiguous case
+ if signal.flags['C'] and axis % 2 == 1:
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ _autocorr.flt_xcorr_nofft_2d(signal, signal.shape[0], signal.shape[1],
+ res, lag)
+ # contiguous case
+ elif signal.flags['F'] and axis % 2 == 0:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype, order = 'F')
+ _autocorr.flt_xcorr_nofft_2d(signal, signal.shape[1], signal.shape[0],
+ res, lag)
+ # non contiguous case
+ elif axis % 2 == 0:
+ res = N.zeros((lag+1, signal.shape[1]), signal.dtype)
+ warn("non contiguous used, this will be slow")
+ _autocorr.flt_xcorr_nofft_2d_noncontiguous(signal,
+ signal.shape[1], signal.shape[0],
+ signal.strides[1], signal.strides[0],
+ res, res.strides[1], res.strides[0], lag)
+ elif axis % 2 == 1:
+ res = N.zeros((signal.shape[0], lag+1), signal.dtype)
+ warn("non contiguous used, this will be slow")
+ _autocorr.flt_xcorr_nofft_2d_noncontiguous(signal,
+ signal.shape[0], signal.shape[1],
+ signal.strides[0], signal.strides[1],
+ res, res.strides[0], res.strides[1], lag)
+ else:
+ raise TypeError("only float 32 and 64 bits supported for now")
+ else:
+ raise RuntimeError("rank > 2 not supported yet")
+
+ return res
+
+def bench():
+ size = 256
+ nframes = 4000
+ lag = 24
+
+ X = N.random.randn(nframes, size)
+ X = N.require(X, requirements = 'C')
+
+ niter = 10
+
+ # Contiguous
+ print "Running optimized with ctypes"
+ def contig(*args, **kargs):
+ return autocorr_oneside_nofft(*args, **kargs)
+ for i in range(niter):
+ Yt = contig(X, lag, axis = 1)
+
+ Yr = _autocorr_oneside_nofft_py(X, lag, axis = 1)
+ N.testing.assert_array_almost_equal(Yt, Yr, 10)
+
+ # Non contiguous
+ print "Running optimized with ctypes (non contiguous)"
+ def ncontig(*args, **kargs):
+ return autocorr_oneside_nofft(*args, **kargs)
+ X = N.require(X, requirements = 'F')
+ for i in range(niter):
+ Yt = ncontig(X, lag, axis = 1)
+
+ Yr = _autocorr_oneside_nofft_py(X, lag, axis = 1)
+ N.testing.assert_array_almost_equal(Yt, Yr, 10)
+
+ print "Benchmark func done"
+
+if __name__ == '__main__':
+ import hotshot, hotshot.stats
+ profile_file = 'autocorr.prof'
+ prof = hotshot.Profile(profile_file, lineevents=1)
+ prof.runcall(bench)
+ p = hotshot.stats.load(profile_file)
+ print p.sort_stats('cumulative').print_stats(20)
+ prof.close()
Added: trunk/Lib/sandbox/cdavid/info.py
===================================================================
--- trunk/Lib/sandbox/cdavid/info.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/info.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,10 @@
+"""
+LPC Routines for speech processing
+
+Copyright: David Cournapeau 2006
+License: BSD-style (see LICENSE.txt in main source directory)
+"""
+version = '0.1'
+
+depends = ['linalg']
+ignore = False
Added: trunk/Lib/sandbox/cdavid/lpc.py
===================================================================
--- trunk/Lib/sandbox/cdavid/lpc.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/lpc.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,249 @@
+#! /usr/bin/env python
+
+# Last Change: Tue Nov 28 03:00 PM 2006 J
+
+# TODO: make lpc work for any rank
+from warnings import warn
+
+import numpy as _N
+import scipy.signal as _sig
+
+from scipy.linalg import toeplitz, inv
+
+from autocorr import autocorr_oneside_nofft
+
+# use ctype to have fast lpc computation
+import ctypes
+from ctypes import c_uint, c_int
+from numpy.ctypeslib import ndpointer, load_library
+ctypes_major = int(ctypes.__version__.split('.')[0])
+if ctypes_major < 1:
+ msg = "version of ctypes is %s, expected at least %s" \
+ % (ctypes.__version__, '1.0.1')
+ raise importerror(msg)
+
+# load lpc
+_lpc = load_library('gabsig.so', __file__)
+
+#===============================
+# define the functions with args
+#===============================
+arg1 = ndpointer(dtype = _N.float64, flags = 'C')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = _N.float64, flags = 'C')
+arg5 = ndpointer(dtype = _N.float64, flags = 'C')
+arg6 = ndpointer(dtype = _N.float64, flags = 'C')
+
+_lpc.dbl_lpc.argtypes = [arg1, arg2, arg3, arg4, arg5, arg6]
+_lpc.dbl_lpc.restype = c_int
+
+arg1 = ndpointer(dtype = _N.float32, flags = 'C')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = _N.float32, flags = 'C')
+arg5 = ndpointer(dtype = _N.float32, flags = 'C')
+arg6 = ndpointer(dtype = _N.float32, flags = 'C')
+
+_lpc.flt_lpc.argtypes = [arg1, arg2, arg3, arg4, arg5, arg6]
+_lpc.flt_lpc.restype = c_int
+
+arg1 = ndpointer(dtype = _N.float64, flags = 'C')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = _N.float64, flags = 'C')
+arg5 = ndpointer(dtype = _N.float64, flags = 'C')
+arg6 = ndpointer(dtype = _N.float64, flags = 'C')
+
+_lpc.dbl_levinson2d.argtypes = [arg1, arg2, arg3, arg4, arg5, arg6]
+_lpc.dbl_levinson2d.restype = c_int
+
+arg1 = ndpointer(dtype = _N.float32, flags = 'C')
+arg2 = c_uint
+arg3 = c_uint
+arg4 = ndpointer(dtype = _N.float32, flags = 'C')
+arg5 = ndpointer(dtype = _N.float32, flags = 'C')
+arg6 = ndpointer(dtype = _N.float32, flags = 'C')
+
+_lpc.flt_levinson2d.argtypes = [arg1, arg2, arg3, arg4, arg5, arg6]
+_lpc.flt_levinson2d.restype = c_int
+
+def lpc_ref(signal, order):
+ """ Return the order + 1 LPC coefficients
+ for the signal. This is just for reference, as it is using
+ the direct inversion of the toeplitz matrix, which is really slow"""
+ if signal.ndim > 1:
+ print "Warning, not tested for rank > 1"
+ p = order + 1
+ r = autocorr_oneside_nofft(signal, order) / signal.shape[0]
+ return _N.concatenate(([1.], _N.dot(inv(toeplitz(r[:-1])), -r[1:])))
+
+def lpc(signal, order):
+ """ Return the order + 1 LPC coefficients
+ for the signal using levinson durbin algorithm """
+ if signal.ndim > 1:
+ warn("Warning, not tested for rank > 1")
+ if signal.size <= order:
+ raise RuntimeError("size is smaller than order !")
+
+ if signal.dtype == _N.float64:
+ coeff = _N.zeros((order+1), _N.float64)
+ kcoeff = _N.zeros((order), _N.float64)
+ err = _N.zeros((1), _N.float64)
+ st = _lpc.dbl_lpc(signal, signal.size, order, coeff, kcoeff, err)
+ elif signal.dtype == _N.float32:
+ coeff = _N.zeros((order+1), _N.float32)
+ kcoeff = _N.zeros((order), _N.float32)
+ err = _N.zeros((1), _N.float32)
+ st = _lpc.flt_lpc(signal, signal.size, order, coeff, kcoeff, err)
+ else:
+ raise TypeError("Sorry, only float32 and float64 supported")
+ if not (st == 0):
+ raise RuntimeError("Error while using levinson algo, returns err is %d", st)
+ return coeff, err, kcoeff
+
+def lpcres_ref(signal, order, axis = 0):
+ return _sig.lfilter(lpc_ref(signal, order), 1., signal, axis = 0)
+
+def lpcres(signal, order, axis = 0):
+ return _sig.lfilter(lpc(signal, order)[0], 1., signal, axis = 0)
+
+def _lpc2_py(signal, order, axis = -1):
+ """python implementation of lpc for rank 2., Do not use, for testing purpose only"""
+ if signal.ndim > 2:
+ raise NotImplemented("only for rank <=2")
+
+ if signal.ndim < 2:
+ return lpc(_N.require(signal, requirements = 'C'), order)
+
+ # For each array of direction axis, compute levinson durbin
+ if axis % 2 == 0:
+ # Prepare output arrays
+ coeff = _N.zeros((order+1, signal.shape[1]), signal.dtype)
+ kcoeff = _N.zeros((order, signal.shape[1]), signal.dtype)
+ err = _N.zeros(signal.shape[1], signal.dtype)
+ for i in range(signal.shape[1]):
+ coeff[:, i], err[i], kcoeff[:, i] = \
+ lpc(_N.require(signal[:, i], requirements = 'C'), order)
+ elif axis % 2 == 1:
+ # Prepare output arrays
+ coeff = _N.zeros((signal.shape[0], order+1), signal.dtype)
+ kcoeff = _N.zeros((signal.shape[0], order), signal.dtype)
+ err = _N.zeros(signal.shape[0], signal.dtype)
+ for i in range(signal.shape[0]):
+ coeff[i], err[i], kcoeff[i] = \
+ lpc(_N.require(signal[i], requirements = 'C'), order)
+ else:
+ raise RuntimeError("this should not happen, please fill a bug")
+
+ return coeff, err, kcoeff
+
+def lpc2(signal, order, axis = -1):
+ """ Returns ar coeff, err and k coeff"""
+ sz = signal.shape[axis]
+ if order >= sz:
+ raise RuntimeError("order should be strictly smaller than the length of signal")
+
+ # rank 1
+ if signal.ndim < 2:
+ if signal.dtype == _N.float64:
+ coeff = _N.zeros((order+1), _N.float64)
+ kcoeff = _N.zeros((order), _N.float64)
+ err = _N.zeros((1), _N.float64)
+ st = _lpc.dbl_lpc(signal, signal.size, order, coeff, kcoeff, err)
+ elif signal.dtype == _N.float32:
+ coeff = _N.zeros((order+1), _N.float32)
+ kcoeff = _N.zeros((order), _N.float32)
+ err = _N.zeros((1), _N.float32)
+ st = _lpc.flt_lpc(signal, signal.size, order, coeff, kcoeff, err)
+ else:
+ raise TypeError("Sorry, only float32 and float64 supported")
+
+ if not (st == 0):
+ raise RuntimeError("Error while using levinson algo, returns err is %d", st)
+
+ return coeff, err, kcoeff
+ # rank 2
+ elif signal.ndim == 2:
+ # Compute biased autocorrelation up to lag = order
+ bias = 1. / sz
+ corr = bias * autocorr_oneside_nofft(signal, order, axis)
+
+ if axis % 2 == 0:
+ # we transpose to have a major row order
+ icorr = corr.T
+ icorr = _N.require(icorr, requirements = 'C')
+ at = _N.zeros((icorr.shape[0], order+1), icorr.dtype)
+ kt = _N.zeros((icorr.shape[0], order), icorr.dtype)
+ et = _N.zeros(icorr.shape[0], icorr.dtype)
+
+ if icorr.dtype == _N.float64:
+ _lpc.dbl_levinson2d(icorr, icorr.shape[0], icorr.shape[1], at, et, kt)
+ elif icorr.dtype == _N.float32:
+ _lpc.flt_levinson2d(icorr, icorr.shape[0], icorr.shape[1], at, et, kt)
+ else:
+ raise TypeError("Only float32 and float64 supported")
+
+ return at.T, et.T, kt.T
+ elif axis % 2 == 1:
+ icorr = _N.require(corr, requirements = 'C')
+ at = _N.zeros((icorr.shape[0], order+1), icorr.dtype)
+ kt = _N.zeros((icorr.shape[0], order), icorr.dtype)
+ et = _N.zeros(icorr.shape[0], icorr.dtype)
+
+ if icorr.dtype == _N.float64:
+ _lpc.dbl_levinson2d(icorr, icorr.shape[0], icorr.shape[1], at, et, kt)
+ elif icorr.dtype == _N.float32:
+ _lpc.flt_levinson2d(icorr, icorr.shape[0], icorr.shape[1], at, et, kt)
+ else:
+ raise TypeError("Only float32 and float64 supported")
+
+ return at, et, kt
+ else:
+ raise RuntimeError("This should not happen, this is a bug")
+ else:
+ raise RuntimeError("Sorry, only rank <= 2 supported for now")
+
+def bench():
+ size = 256
+ nframes = 4000
+ order = 24
+
+ X = _N.random.randn(nframes, size)
+ X = _N.require(X, requirements = 'C')
+
+ niter = 10
+
+ # Contiguous
+ print "Running optimized with ctypes"
+ for i in range(niter):
+ at, et, kt = lpc2(X, order, axis = 1)
+
+ a, e, k = _lpc2_py(X, order, axis = 1)
+ _N.testing.assert_array_almost_equal(a, at, 10)
+ _N.testing.assert_array_almost_equal(e, et, 10)
+ _N.testing.assert_array_almost_equal(k, kt, 10)
+
+ ## Non contiguous
+ #print "Running optimized with ctypes (non contiguous)"
+ #def ncontig(*args, **kargs):
+ # return lpc2(*args, **kargs)
+ #X = _N.require(X, requirements = 'F')
+ #for i in range(niter):
+ # at, et, kt = ncontig(X, order, axis = 1)
+
+ #a, e, k = _lpc2_py(X, order, axis = 1)
+ #_N.testing.assert_array_almost_equal(a, at, 10)
+ #_N.testing.assert_array_almost_equal(e, et, 10)
+ #_N.testing.assert_array_almost_equal(k, kt, 10)
+
+ print "Benchmark func done"
+
+if __name__ == '__main__':
+ import hotshot, hotshot.stats
+ profile_file = 'lpc.prof'
+ prof = hotshot.Profile(profile_file, lineevents=1)
+ prof.runcall(bench)
+ p = hotshot.stats.load(profile_file)
+ print p.sort_stats('cumulative').print_stats(20)
+ prof.close()
Added: trunk/Lib/sandbox/cdavid/setup.py
===================================================================
--- trunk/Lib/sandbox/cdavid/setup.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/setup.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,41 @@
+#! /usr/bin/env python
+# Last Change: Tue Nov 28 03:00 PM 2006 J
+
+""" toolbox of Cournapeau David. Implements various things such as
+autocorrelation with lag argument, lpc coefficients computation.
+
+2006, David Cournapeau
+
+LICENSE: the license of pyaudio is the same than scipy"""
+
+from os.path import join
+import os
+
+from info import version as cdavid_version
+
+DISTNAME = 'cdavid'
+VERSION = cdavid_version
+DESCRIPTION ='A scipy package for various speech processing tools lpc'
+MAINTAINER ='David Cournapeau',
+MAINTAINER_EMAIL='david at ar.media.kyoto-u.ac.jp',
+URL ='http://ar.media.kyoto-u.ac.jp/members/david',
+LICENSE = 'BSD'
+
+def configuration(parent_package='',top_path=None, package_name=DISTNAME):
+ from numpy.distutils.misc_util import Configuration
+ config = Configuration(package_name,parent_package,top_path,
+ version = VERSION)
+ config.add_data_dir('tests')
+ #config.add_data_dir('profile_data')
+ config.add_extension('gabsig',
+ #define_macros=[('LIBSVM_EXPORTS', None),
+ # ('LIBSVM_DLL', None)],
+ sources=[join('src', 'levinson.c'),
+ join('src', 'autocorr_nofft.c'),
+ join('src', 'lpc.c')])
+
+ return config
+
+if __name__ == "__main__":
+ from numpy.distutils.core import setup
+ setup(configuration=configuration)
Added: trunk/Lib/sandbox/cdavid/src/Makefile
===================================================================
--- trunk/Lib/sandbox/cdavid/src/Makefile 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/Makefile 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,64 @@
+CC = colorgcc
+LD = gcc
+
+AUTOGEN = autogen
+
+PYTHONINC = -I/usr/include/python2.4
+NUMPYINC = -I/home/david/local/lib/python2.4/site-packages/numpy/core/include
+
+DEBUG = -g
+OPTIMS = -O3 -funroll-all-loops -march=pentium4 -msse2
+#OPTIMS = $(DEBUG)
+WARN = -W -Wall -Winline -Wstrict-prototypes -Wmissing-prototypes -Waggregate-return -Wcast-align -Wcast-qual -Wnested-externs -Wshadow -Wbad-function-cast -Wwrite-strings
+LANGSTD = -std=gnu9x
+
+CFLAGS = $(PYTHONINC) $(NUMPYINC) $(OPTIMS) $(WARN) $(LANGSTD)
+
+OBJS = autocorr_nofft.o levinson.o lpc.o
+
+SRC = levinson.c lpc.c autocorr_nofft.c
+
+# targets
+lib: libgabsig.so
+
+src: $(SRC)
+
+#=================
+# Shared libraries
+#=================
+libgabsig.so: $(OBJS)
+ $(LD) -shared -o $@ $(OBJS) -Wl,-soname,$@
+
+#=============
+# Object files
+#=============
+lpc.o: lpc.c
+ $(CC) -c $(CFLAGS) -fPIC $<
+
+levinson.o: levinson.c
+ $(CC) -c $(CFLAGS) -fPIC $<
+
+autocorr_nofft.o: autocorr_nofft.c autocorr_nofft.h
+ $(CC) -c $(CFLAGS) -fPIC $<
+
+#=========================
+# Autogenerated c sources
+#=========================
+autocorr_nofft.c: autocorr_nofft.def autocorr_nofft.tpl
+ $(AUTOGEN) autocorr_nofft.def
+
+lpc.c: lpc.def lpc.tpl lpc.h
+ $(AUTOGEN) lpc.def
+
+levinson.c: levinson.def levinson.tpl levinson.h
+ $(AUTOGEN) levinson.def
+
+clean:
+ rm -f libautocorr.so
+ rm -f test_autocorr
+ rm -f *.o
+ rm -f *.so
+ rm -f *.pyc
+ rm -f autocorr_nofft.c
+ rm -f lpc.c
+ rm -f levinson.c
Added: trunk/Lib/sandbox/cdavid/src/README
===================================================================
--- trunk/Lib/sandbox/cdavid/src/README 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/README 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,6 @@
+C sources are generated by autogen. This enables the support of float and double with the
+same sources (ala C++ template).
+
+To generate the sources, just execute make src.
+To modify a source, change the file .tpl instead of the .c, and regenerate the source
+(using make src, for example)
Added: trunk/Lib/sandbox/cdavid/src/autocorr_nofft.c
===================================================================
--- trunk/Lib/sandbox/cdavid/src/autocorr_nofft.c 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/autocorr_nofft.c 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,299 @@
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ * vim:syntax=c
+ *
+ * TODO: is size_t 64 bits long on 64 bits machines ?
+ */
+#include <stddef.h> /* for size_t */
+#include <stdio.h> /* for size_t */
+
+#include "autocorr_nofft.h"
+
+/*
+ * NOFFT auto correlation
+ */
+
+
+/*
+ * float version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int flt_xcorr_nofft_1d(const float *in,
+ const size_t size, float *out, const size_t lag)
+{
+ size_t i, j;
+ float acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < size; ++i) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag; ++i) {
+ acc = 0;
+ for (j = i; j < size; ++j) {
+ acc += in[j-i]*in[j];
+ }
+ out[i] = acc;
+ }
+
+ return 0;
+}
+
+/*
+ * double version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int dbl_xcorr_nofft_1d(const double *in,
+ const size_t size, double *out, const size_t lag)
+{
+ size_t i, j;
+ double acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < size; ++i) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag; ++i) {
+ acc = 0;
+ for (j = i; j < size; ++j) {
+ acc += in[j-i]*in[j];
+ }
+ out[i] = acc;
+ }
+
+ return 0;
+}
+
+
+
+/*
+ * float version for non contiguous arrays; the corresponding
+ * array should have at least in_size elements.
+ *
+ * Constraints:
+ * - lag should be < in_size
+ * - strides in bytes
+ * - TODO: check if should be aligned ?
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int flt_xcorr_nofft_1d_noncontiguous(const float *in, size_t in_size,
+ size_t in_stride, float *out, size_t out_stride, size_t lag)
+{
+ size_t i, j, clag;
+ size_t istride = in_stride / sizeof(float);
+ size_t ostride = out_stride / sizeof(float);
+ float acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < in_size * istride; i+= istride) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag ; ++i) {
+ acc = 0;
+ clag = i * istride;
+ for (j = clag; j < in_size * istride; j += istride) {
+ acc += in[j-clag]*in[j];
+ }
+ out[i * ostride] = acc;
+ }
+
+ return 0;
+}
+
+/*
+ * double version for non contiguous arrays; the corresponding
+ * array should have at least in_size elements.
+ *
+ * Constraints:
+ * - lag should be < in_size
+ * - strides in bytes
+ * - TODO: check if should be aligned ?
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int dbl_xcorr_nofft_1d_noncontiguous(const double *in, size_t in_size,
+ size_t in_stride, double *out, size_t out_stride, size_t lag)
+{
+ size_t i, j, clag;
+ size_t istride = in_stride / sizeof(double);
+ size_t ostride = out_stride / sizeof(double);
+ double acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < in_size * istride; i+= istride) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag ; ++i) {
+ acc = 0;
+ clag = i * istride;
+ for (j = clag; j < in_size * istride; j += istride) {
+ acc += in[j-clag]*in[j];
+ }
+ out[i * ostride] = acc;
+ }
+
+ return 0;
+}
+
+
+
+/*
+ * For rank 2 arrays, contiguous cases
+ * float version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int flt_xcorr_nofft_2d(const float *in,
+ size_t dim0, size_t dim1, float *out, const size_t lag)
+{
+ size_t i;
+ float *coaxis;
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d autocorr, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * (lag + 1);
+ flt_xcorr_nofft_1d(in + i * dim1, dim1, coaxis, lag);
+ }
+
+ return 0;
+}
+
+/*
+ * For rank 2 arrays, contiguous cases
+ * double version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int dbl_xcorr_nofft_2d(const double *in,
+ size_t dim0, size_t dim1, double *out, const size_t lag)
+{
+ size_t i;
+ double *coaxis;
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d autocorr, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * (lag + 1);
+ dbl_xcorr_nofft_1d(in + i * dim1, dim1, coaxis, lag);
+ }
+
+ return 0;
+}
+
+
+
+/*
+ * For rank 2 arrays, non contiguous cases
+ * float version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int flt_xcorr_nofft_2d_noncontiguous(const float *in,
+ size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1,
+ float *out, size_t out_stride0, size_t out_stride1,
+ const size_t lag)
+{
+ size_t i;
+
+ size_t istride0 = in_stride0 / sizeof(float);
+ size_t ostride0 = out_stride0 / sizeof(float);
+
+ float *coaxis;
+#if 0
+ fprintf(stdout, "%s: shape is (%d, %d)\n", __func__, dim0, dim1);
+ fprintf(stdout, "%s: istrides are (%d, %d)\n", __func__, istride0, istride1);
+
+ fprintf(stdout, "%s: ostrides are (%d, %d)\n", __func__, ostride0, ostride1);
+ for(i = 0; i < dim0; ++i) {
+ ciaxis = in + i * istride0;
+ coaxis = out + i * istride0;
+ fprintf(stdout, "%d 1d autocorr, first element is %f, last is %f (%d el)\n",
+ i, ciaxis[0], ciaxis[(dim1-1) * istride1], dim1);
+ }
+#endif
+
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * ostride0;
+ flt_xcorr_nofft_1d_noncontiguous(in + i * istride0, dim1, in_stride1,
+ coaxis, out_stride1, lag);
+ }
+ return 0;
+}
+
+/*
+ * For rank 2 arrays, non contiguous cases
+ * double version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int dbl_xcorr_nofft_2d_noncontiguous(const double *in,
+ size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1,
+ double *out, size_t out_stride0, size_t out_stride1,
+ const size_t lag)
+{
+ size_t i;
+
+ size_t istride0 = in_stride0 / sizeof(double);
+ size_t ostride0 = out_stride0 / sizeof(double);
+
+ double *coaxis;
+#if 0
+ fprintf(stdout, "%s: shape is (%d, %d)\n", __func__, dim0, dim1);
+ fprintf(stdout, "%s: istrides are (%d, %d)\n", __func__, istride0, istride1);
+
+ fprintf(stdout, "%s: ostrides are (%d, %d)\n", __func__, ostride0, ostride1);
+ for(i = 0; i < dim0; ++i) {
+ ciaxis = in + i * istride0;
+ coaxis = out + i * istride0;
+ fprintf(stdout, "%d 1d autocorr, first element is %f, last is %f (%d el)\n",
+ i, ciaxis[0], ciaxis[(dim1-1) * istride1], dim1);
+ }
+#endif
+
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * ostride0;
+ dbl_xcorr_nofft_1d_noncontiguous(in + i * istride0, dim1, in_stride1,
+ coaxis, out_stride1, lag);
+ }
+ return 0;
+}
+
+
Added: trunk/Lib/sandbox/cdavid/src/autocorr_nofft.def
===================================================================
--- trunk/Lib/sandbox/cdavid/src/autocorr_nofft.def 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/autocorr_nofft.def 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,18 @@
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ * vim:syntax=c
+ */
+autogen definitions autocorr_nofft.tpl;
+
+float_type = {
+ type_name = "float" ;
+ short_name = "flt" ;
+ upper_name = "FLOAT" ;
+} ;
+
+float_type = {
+ type_name = "double" ;
+ short_name = "dbl" ;
+ upper_name = "DOUBLE" ;
+} ;
+
Added: trunk/Lib/sandbox/cdavid/src/autocorr_nofft.h
===================================================================
--- trunk/Lib/sandbox/cdavid/src/autocorr_nofft.h 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/autocorr_nofft.h 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,47 @@
+/*
+ * Last Change: Mon Nov 27 07:00 PM 2006 J
+ */
+
+#ifndef _C_AUTOCORR_NOFFT_H_
+#define _C_AUTOCORR_NOFFT_H_
+
+/*
+ * Direct implementation of auto correlation (faster when a few lags only are
+ * necessary). One side only, out should have at least lag+1 elements allocated
+ *
+ * Expect in and out to be contiguous
+ */
+int flt_xcorr_nofft_1d(const float *in, size_t size, float *out, size_t lag);
+int dbl_xcorr_nofft_1d(const double *in, size_t size, double *out, size_t lag);
+
+/*
+ * Direct implementation of auto correlation (faster when a few lags only are
+ * necessary). One side only, out should have at least lag+1 elements allocated
+ *
+ * Expect in and out need not to be contiguous
+ */
+int flt_xcorr_nofft_1d_noncontiguous(const float *in, size_t in_size, size_t in_stride,
+ float *out, size_t out_stride, size_t lag);
+int dbl_xcorr_nofft_1d_noncontiguous(const double *in, size_t in_size, size_t in_stride,
+ double *out, size_t out_stride, size_t lag);
+
+/*
+ * 1d autocorrelation for rank 2 arrays
+ */
+int flt_xcorr_nofft_2d(const float *in, size_t dim0, size_t dim1,
+ float *out, size_t lag);
+int dbl_xcorr_nofft_2d(const double *in, size_t dim0, size_t dim1,
+ double *out, size_t lag);
+
+/*
+ * 1d autocorrelation for rank 2 arrays, non contiguous cases
+ */
+int dbl_xcorr_nofft_2d_noncontiguous(const double *in,
+ size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1,
+ double *out, size_t out_stride0, size_t out_stride1,
+ const size_t lag);
+int flt_xcorr_nofft_2d_noncontiguous(const float *in,
+ size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1,
+ float *out, size_t out_stride0, size_t out_stride1,
+ const size_t lag);
+#endif
Added: trunk/Lib/sandbox/cdavid/src/autocorr_nofft.tpl
===================================================================
--- trunk/Lib/sandbox/cdavid/src/autocorr_nofft.tpl 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/autocorr_nofft.tpl 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,162 @@
+[+ AutoGen5 template c +]
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ * vim:syntax=c
+ *
+ * TODO: is size_t 64 bits long on 64 bits machines ?
+ */
+#include <stddef.h> /* for size_t */
+#include <stdio.h> /* for size_t */
+
+#include "autocorr_nofft.h"
+
+/*
+ * NOFFT auto correlation
+ */
+
+[+ For float_type +]
+/*
+ * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int [+ (get "short_name") +]_xcorr_nofft_1d(const [+ (get "type_name") +] *in,
+ const size_t size, [+ (get "type_name") +] *out, const size_t lag)
+{
+ size_t i, j;
+ [+ (get "type_name") +] acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < size; ++i) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag; ++i) {
+ acc = 0;
+ for (j = i; j < size; ++j) {
+ acc += in[j-i]*in[j];
+ }
+ out[i] = acc;
+ }
+
+ return 0;
+}
+[+ ENDFOR float_type +]
+
+[+ For float_type +]
+/*
+ * [+ (get "type_name") +] version for non contiguous arrays; the corresponding
+ * array should have at least in_size elements.
+ *
+ * Constraints:
+ * - lag should be < in_size
+ * - strides in bytes
+ * - TODO: check if should be aligned ?
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int [+ (get "short_name") +]_xcorr_nofft_1d_noncontiguous(const [+ (get "type_name") +] *in, size_t in_size,
+ size_t in_stride, [+ (get "type_name") +] *out, size_t out_stride, size_t lag)
+{
+ size_t i, j, clag;
+ size_t istride = in_stride / sizeof([+ (get "type_name") +]);
+ size_t ostride = out_stride / sizeof([+ (get "type_name") +]);
+ [+ (get "type_name") +] acc;
+
+ /* lag 0 */
+ acc = 0;
+ for (i = 0; i < in_size * istride; i+= istride) {
+ acc += in[i]*in[i];
+ }
+ out[0] = acc;
+
+ /* lag : 1 -> lag */
+ for (i = 1; i <= lag ; ++i) {
+ acc = 0;
+ clag = i * istride;
+ for (j = clag; j < in_size * istride; j += istride) {
+ acc += in[j-clag]*in[j];
+ }
+ out[i * ostride] = acc;
+ }
+
+ return 0;
+}
+[+ ENDFOR float_type +]
+
+[+ For float_type +]
+/*
+ * For rank 2 arrays, contiguous cases
+ * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int [+ (get "short_name") +]_xcorr_nofft_2d(const [+ (get "type_name") +] *in,
+ size_t dim0, size_t dim1, [+ (get "type_name") +] *out, const size_t lag)
+{
+ size_t i;
+ [+ (get "type_name") +] *coaxis;
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d autocorr, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * (lag + 1);
+ [+ (get "short_name") +]_xcorr_nofft_1d(in + i * dim1, dim1, coaxis, lag);
+ }
+
+ return 0;
+}
+[+ ENDFOR float_type +]
+
+[+ For float_type +]
+/*
+ * For rank 2 arrays, non contiguous cases
+ * [+ (get "type_name") +] version; out must have a size of lag+1, already pre allocated
+ *
+ * lag should be < size
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int [+ (get "short_name") +]_xcorr_nofft_2d_noncontiguous(const [+ (get "type_name") +] *in,
+ size_t dim0, size_t dim1, size_t in_stride0, size_t in_stride1,
+ [+ (get "type_name") +] *out, size_t out_stride0, size_t out_stride1,
+ const size_t lag)
+{
+ size_t i;
+
+ size_t istride0 = in_stride0 / sizeof([+ (get "type_name") +]);
+ size_t ostride0 = out_stride0 / sizeof([+ (get "type_name") +]);
+
+ [+ (get "type_name") +] *coaxis;
+#if 0
+ fprintf(stdout, "%s: shape is (%d, %d)\n", __func__, dim0, dim1);
+ fprintf(stdout, "%s: istrides are (%d, %d)\n", __func__, istride0, istride1);
+
+ fprintf(stdout, "%s: ostrides are (%d, %d)\n", __func__, ostride0, ostride1);
+ for(i = 0; i < dim0; ++i) {
+ ciaxis = in + i * istride0;
+ coaxis = out + i * istride0;
+ fprintf(stdout, "%d 1d autocorr, first element is %f, last is %f (%d el)\n",
+ i, ciaxis[0], ciaxis[(dim1-1) * istride1], dim1);
+ }
+#endif
+
+ for(i = 0; i < dim0; ++i) {
+ coaxis = out + i * ostride0;
+ [+ (get "short_name") +]_xcorr_nofft_1d_noncontiguous(in + i * istride0, dim1, in_stride1,
+ coaxis, out_stride1, lag);
+ }
+ return 0;
+}
+[+ ENDFOR float_type +]
+
Added: trunk/Lib/sandbox/cdavid/src/common.h
===================================================================
--- trunk/Lib/sandbox/cdavid/src/common.h 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/common.h 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,38 @@
+/*
+ * Last Change: Tue Nov 28 04:00 PM 2006 J
+ *
+ * Implements FP macros missing in C89
+ */
+#ifndef _GABSIG_C_COMMON_H
+ #define _GABSIG_C_COMMON_H
+
+#include <math.h>
+#if defined(fpclassify)
+
+ #if !defined(isnan)
+ #define isnan(x) (fpclassify((x)) == FP_NAN)
+ #endif
+ #if !defined(isinf)
+ #define isinf(x) (fpclassify((x)) == FP_INFINITE)
+ #endif
+
+#else /* check to see if already have a function like this */
+ #if !defined(HAVE_ISNAN)
+ #if !defined(isnan)
+ #define isnan(x) ((x) == (x))
+ #endif
+ #endif /* HAVE_ISNAN */
+
+ #if !defined(HAVE_ISINF)
+ #if !defined(isinf)
+ #define isinf(x) (!isnan((x)) && isnan((x)-(x)))
+ #endif
+ #endif /* HAVE_ISINF */
+
+#endif /* defined(fpclassify) */
+
+#if !defined(isfinite)
+ #define isfinite(x) (!isnan((x)) && !isinf((x)))
+#endif
+
+#endif /* nef of recursive header inclusion protection */
Added: trunk/Lib/sandbox/cdavid/src/levinson.c
===================================================================
--- trunk/Lib/sandbox/cdavid/src/levinson.c 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/levinson.c 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,313 @@
+/*
+ * Last Change: Tue Nov 28 04:00 PM 2006 J
+ *
+ * vim:syntax=c
+ */
+#include <float.h>
+#include <math.h> /* for isfinite */
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "common.h"
+
+#include "levinson.h"
+
+
+/* The actual computation :
+ * - in : the input vector which defines the toeplitz matrix
+ * - size : size of in (ie number of elements)
+ * - order : size of the system to solve. order must be < size -1
+ * - acoeff: solution (ie ar coefficients). Size must be at last order+1
+ * - err : *prediction* error (scalar)
+ * - kcoeff: reflexion coefficients. Size must be at last equal to equal to order.
+ * - tmp : cache, mnust have at least order elements
+ */
+
+/*
+ * this function assume all arrays are allocated with the
+ * right size, and that the parameters make sense. No checking
+ * is done, must be done before calling this function
+ *
+ * Returns 0 on success, -1 if a compuation error happened (overflow, underflow
+ * for error calculation)
+ */
+
+int flt_levinson1d(const float* in,
+ size_t order, float* acoeff,
+ float* err, float* kcoeff,
+ float* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ float acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = (float)1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+int flt_levinson1d_check(const float* in,
+ size_t order, float* acoeff,
+ float* err, float* kcoeff,
+ float* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ float acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = (float)1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ if (!isfinite(kcoeff[i-1])) {
+ fprintf(stderr, "%s:%s, kcoeff is not finite, err is %e\n",
+ __FILE__, __func__, *err);
+ ret = -1;
+ }
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+/*
+ * For rank 2 arrays, contiguous cases
+ * float version; out must have a size of dim0 * (lag + 1),
+ * already pre allocated
+ *
+ * order should be < dim1
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int flt_levinson2d(const float *in,
+ size_t dim0, size_t dim1,
+ float *acoeff,
+ float *err,
+ float *kcoeff)
+{
+ size_t i;
+ size_t order = dim1 - 1;
+ float *caaxis, *ceaxis, *ckaxis, *buff;
+
+ buff = malloc(sizeof(*buff) * order);
+ if (buff == NULL) {
+ goto fail_malloc;
+ }
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d levinson, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ caaxis = acoeff + i * (order + 1);
+ ckaxis = kcoeff + i * (order);
+ ceaxis = err + i;
+ flt_levinson1d(in + i * dim1, order, caaxis, ceaxis, ckaxis, buff);
+ }
+
+ free(buff);
+ return 0;
+
+fail_malloc:
+ return -1;
+}
+
+/* The actual computation :
+ * - in : the input vector which defines the toeplitz matrix
+ * - size : size of in (ie number of elements)
+ * - order : size of the system to solve. order must be < size -1
+ * - acoeff: solution (ie ar coefficients). Size must be at last order+1
+ * - err : *prediction* error (scalar)
+ * - kcoeff: reflexion coefficients. Size must be at last equal to equal to order.
+ * - tmp : cache, mnust have at least order elements
+ */
+
+/*
+ * this function assume all arrays are allocated with the
+ * right size, and that the parameters make sense. No checking
+ * is done, must be done before calling this function
+ *
+ * Returns 0 on success, -1 if a compuation error happened (overflow, underflow
+ * for error calculation)
+ */
+
+int dbl_levinson1d(const double* in,
+ size_t order, double* acoeff,
+ double* err, double* kcoeff,
+ double* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ double acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = (double)1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+int dbl_levinson1d_check(const double* in,
+ size_t order, double* acoeff,
+ double* err, double* kcoeff,
+ double* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ double acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = (double)1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ if (!isfinite(kcoeff[i-1])) {
+ fprintf(stderr, "%s:%s, kcoeff is not finite, err is %e\n",
+ __FILE__, __func__, *err);
+ ret = -1;
+ }
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+/*
+ * For rank 2 arrays, contiguous cases
+ * double version; out must have a size of dim0 * (lag + 1),
+ * already pre allocated
+ *
+ * order should be < dim1
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int dbl_levinson2d(const double *in,
+ size_t dim0, size_t dim1,
+ double *acoeff,
+ double *err,
+ double *kcoeff)
+{
+ size_t i;
+ size_t order = dim1 - 1;
+ double *caaxis, *ceaxis, *ckaxis, *buff;
+
+ buff = malloc(sizeof(*buff) * order);
+ if (buff == NULL) {
+ goto fail_malloc;
+ }
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d levinson, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ caaxis = acoeff + i * (order + 1);
+ ckaxis = kcoeff + i * (order);
+ ceaxis = err + i;
+ dbl_levinson1d(in + i * dim1, order, caaxis, ceaxis, ckaxis, buff);
+ }
+
+ free(buff);
+ return 0;
+
+fail_malloc:
+ return -1;
+}
+
Added: trunk/Lib/sandbox/cdavid/src/levinson.def
===================================================================
--- trunk/Lib/sandbox/cdavid/src/levinson.def 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/levinson.def 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,18 @@
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ * vim:syntax=c
+ */
+autogen definitions levinson.tpl;
+
+float_type = {
+ type_name = "float" ;
+ short_name = "flt" ;
+ upper_name = "FLOAT" ;
+} ;
+
+float_type = {
+ type_name = "double" ;
+ short_name = "dbl" ;
+ upper_name = "DOUBLE" ;
+} ;
+
Added: trunk/Lib/sandbox/cdavid/src/levinson.h
===================================================================
--- trunk/Lib/sandbox/cdavid/src/levinson.h 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/levinson.h 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,24 @@
+#ifndef LEVINSON1D_H
+#define LEVINSON1D_H
+
+#include <stddef.h>
+
+int dbl_levinson1d(const double* in, size_t order,
+ double* acoeff, double *err, double* kcoeff, double* tmp);
+
+int dbl_levinson1d_check(const double* in, size_t order,
+ double* acoeff, double *err, double* kcoeff, double* tmp);
+
+int flt_levinson1d(const float* in, size_t order,
+ float* acoeff, float *err, float* kcoeff, float* tmp);
+
+int flt_levinson1d_check(const float* in, size_t order,
+ float* acoeff, float *err, float* kcoeff, float* tmp);
+
+int dbl_levinson2d(const double* in, size_t dim0, size_t dim1,
+ double* acoeff, double *err, double* kcoeff);
+
+int flt_levinson2d(const float* in, size_t dim0, size_t dim1,
+ float* acoeff, float *err, float* kcoeff);
+
+#endif
Added: trunk/Lib/sandbox/cdavid/src/levinson.tpl
===================================================================
--- trunk/Lib/sandbox/cdavid/src/levinson.tpl 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/levinson.tpl 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,165 @@
+[+ AutoGen5 template c +]
+/*
+ * Last Change: Tue Nov 28 04:00 PM 2006 J
+ *
+ * vim:syntax=c
+ */
+#include <float.h>
+#include <math.h> /* for isfinite */
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "common.h"
+
+#include "levinson.h"
+
+[+ For float_type +]
+/* The actual computation :
+ * - in : the input vector which defines the toeplitz matrix
+ * - size : size of in (ie number of elements)
+ * - order : size of the system to solve. order must be < size -1
+ * - acoeff: solution (ie ar coefficients). Size must be at last order+1
+ * - err : *prediction* error (scalar)
+ * - kcoeff: reflexion coefficients. Size must be at last equal to equal to order.
+ * - tmp : cache, mnust have at least order elements
+ */
+
+/*
+ * this function assume all arrays are allocated with the
+ * right size, and that the parameters make sense. No checking
+ * is done, must be done before calling this function
+ *
+ * Returns 0 on success, -1 if a compuation error happened (overflow, underflow
+ * for error calculation)
+ */
+
+int [+ (get "short_name") +]_levinson1d(const [+ (get "type_name") +]* in,
+ size_t order, [+ (get "type_name") +]* acoeff,
+ [+ (get "type_name") +]* err, [+ (get "type_name") +]* kcoeff,
+ [+ (get "type_name") +]* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ [+ (get "type_name") +] acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = ([+ (get "type_name") +])1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+int [+ (get "short_name") +]_levinson1d_check(const [+ (get "type_name") +]* in,
+ size_t order, [+ (get "type_name") +]* acoeff,
+ [+ (get "type_name") +]* err, [+ (get "type_name") +]* kcoeff,
+ [+ (get "type_name") +]* tmp)
+{
+ /* TODO: to check if first element of corr is 0*/
+
+ size_t i, j;
+ [+ (get "type_name") +] acc;
+ int ret = 0;
+
+ /*
+ * order 0
+ */
+ acoeff[0] = ([+ (get "type_name") +])1.0;
+ *err = in[0];
+
+ /*
+ * order >= 1
+ */
+ for ( i = 1; i <= order; ++i) {
+ acc = in[i];
+ for ( j = 1; j <= i-1; ++j) {
+ acc += acoeff[j]*in[i-j];
+ }
+ kcoeff[i-1] = -acc/(*err);
+ if (!isfinite(kcoeff[i-1])) {
+ fprintf(stderr, "%s:%s, kcoeff is not finite, err is %e\n",
+ __FILE__, __func__, *err);
+ ret = -1;
+ }
+ acoeff[i] = kcoeff[i-1];
+
+ for ( j = 0; j < order; ++j) {
+ tmp[j] = acoeff[j];
+ }
+
+ for (j = 1; j < i; ++j) {
+ acoeff[j] += kcoeff[i-1]*tmp[i-j];
+ }
+ *err *= (1-kcoeff[i-1]*kcoeff[i-1]);
+ }
+
+ return ret;
+}
+
+/*
+ * For rank 2 arrays, contiguous cases
+ * [+ (get "type_name") +] version; out must have a size of dim0 * (lag + 1),
+ * already pre allocated
+ *
+ * order should be < dim1
+ *
+ * returns 0 is succesfull, other value otherwise.
+ */
+int [+ (get "short_name") +]_levinson2d(const [+ (get "type_name") +] *in,
+ size_t dim0, size_t dim1,
+ [+ (get "type_name") +] *acoeff,
+ [+ (get "type_name") +] *err,
+ [+ (get "type_name") +] *kcoeff)
+{
+ size_t i;
+ size_t order = dim1 - 1;
+ [+ (get "type_name") +] *caaxis, *ceaxis, *ckaxis, *buff;
+
+ buff = malloc(sizeof(*buff) * order);
+ if (buff == NULL) {
+ goto fail_malloc;
+ }
+
+#if 0
+ for(i = 0; i < dim0; ++i) {
+ fprintf(stdout, "%d 1d levinson, first element is %f\n", i, in[i * dim1]);
+ }
+#endif
+ for(i = 0; i < dim0; ++i) {
+ caaxis = acoeff + i * (order + 1);
+ ckaxis = kcoeff + i * (order);
+ ceaxis = err + i;
+ [+ (get "short_name") +]_levinson1d(in + i * dim1, order, caaxis, ceaxis, ckaxis, buff);
+ }
+
+ free(buff);
+ return 0;
+
+fail_malloc:
+ return -1;
+}
+[+ ENDFOR float_type +]
Added: trunk/Lib/sandbox/cdavid/src/lpc.c
===================================================================
--- trunk/Lib/sandbox/cdavid/src/lpc.c 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/lpc.c 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,144 @@
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ *
+ * vim:syntax=c
+ */
+#include <stdlib.h> /* for malloc and co */
+#include <stdio.h>
+
+#include "levinson.h"
+#include "autocorr_nofft.h"
+
+#include "lpc.h"
+
+
+/*
+ * (float version) Compute lpc coeff (order x coefficients) of a
+ * contiguous array
+ *
+ * err is a double, coeff must be able to contain order+1 elements, and kcoeff
+ * order elements
+ */
+int flt_lpc(const float* signal,
+ size_t size, size_t order, float* coeff,
+ float* kcoeff, float* err)
+{
+ size_t i, nbuff, ncache;
+ float *buff, *cache, biasnorm;
+ int status;
+
+ biasnorm = 1.0/size;
+ nbuff = order + 1;
+ ncache = order;
+
+ buff = malloc(sizeof(*buff) * nbuff);
+ if (buff == NULL) {
+ status = -2;
+ goto fail_buff_malloc;
+ }
+
+ cache = malloc(sizeof(*cache) * ncache);
+ if (cache == NULL) {
+ status = -2;
+ goto fail_cache_malloc;
+ }
+
+ /*
+ * Compute the autocorreleation up to lag order, normalized by the
+ * size of the signal
+ */
+ flt_xcorr_nofft_1d(signal, size, buff, order);
+ for(i = 0; i < nbuff; ++i) {
+ buff[i] *= biasnorm;
+ }
+
+ /*
+ * Compute the inverse coefficients using (simple) levinson recursive algo
+ */
+ status = flt_levinson1d(buff, order,
+ coeff, err, kcoeff, cache);
+ if (status) {
+ status = -1;
+ goto fail_levinson;
+ }
+
+ free(cache);
+ free(buff);
+
+ return 0;
+
+fail_levinson:
+ free(cache);
+fail_cache_malloc:
+ free(buff);
+fail_buff_malloc:
+ fprintf(stderr, "Failure\n");
+ return status;
+}
+
+
+/*
+ * (double version) Compute lpc coeff (order x coefficients) of a
+ * contiguous array
+ *
+ * err is a double, coeff must be able to contain order+1 elements, and kcoeff
+ * order elements
+ */
+int dbl_lpc(const double* signal,
+ size_t size, size_t order, double* coeff,
+ double* kcoeff, double* err)
+{
+ size_t i, nbuff, ncache;
+ double *buff, *cache, biasnorm;
+ int status;
+
+ biasnorm = 1.0/size;
+ nbuff = order + 1;
+ ncache = order;
+
+ buff = malloc(sizeof(*buff) * nbuff);
+ if (buff == NULL) {
+ status = -2;
+ goto fail_buff_malloc;
+ }
+
+ cache = malloc(sizeof(*cache) * ncache);
+ if (cache == NULL) {
+ status = -2;
+ goto fail_cache_malloc;
+ }
+
+ /*
+ * Compute the autocorreleation up to lag order, normalized by the
+ * size of the signal
+ */
+ dbl_xcorr_nofft_1d(signal, size, buff, order);
+ for(i = 0; i < nbuff; ++i) {
+ buff[i] *= biasnorm;
+ }
+
+ /*
+ * Compute the inverse coefficients using (simple) levinson recursive algo
+ */
+ status = dbl_levinson1d(buff, order,
+ coeff, err, kcoeff, cache);
+ if (status) {
+ status = -1;
+ goto fail_levinson;
+ }
+
+ free(cache);
+ free(buff);
+
+ return 0;
+
+fail_levinson:
+ free(cache);
+fail_cache_malloc:
+ free(buff);
+fail_buff_malloc:
+ fprintf(stderr, "Failure\n");
+ return status;
+}
+
+
Added: trunk/Lib/sandbox/cdavid/src/lpc.def
===================================================================
--- trunk/Lib/sandbox/cdavid/src/lpc.def 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/lpc.def 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,18 @@
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ * vim:syntax=c
+ */
+autogen definitions lpc.tpl;
+
+float_type = {
+ type_name = "float" ;
+ short_name = "flt" ;
+ upper_name = "FLOAT" ;
+} ;
+
+float_type = {
+ type_name = "double" ;
+ short_name = "dbl" ;
+ upper_name = "DOUBLE" ;
+} ;
+
Added: trunk/Lib/sandbox/cdavid/src/lpc.h
===================================================================
--- trunk/Lib/sandbox/cdavid/src/lpc.h 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/lpc.h 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,20 @@
+/*
+ * Last Change: Tue Nov 28 12:00 PM 2006 J
+ */
+
+/*
+ * Functions to compute lpc coefficients with contiguous arrays.
+ *
+ * input is signal, output are coeff, kcoeff and err.
+ *
+ * requirements:
+ * - signal must have size elements
+ * - order < size
+ * - coeff must have order + 1 elements at least
+ * - kcoeff must have order elements at least
+ * - err must have at least one element
+ */
+int dbl_lpc(const double* signal, size_t size, size_t order, double* coeff,
+ double* kcoeff, double* err);
+int flt_lpc(const float* signal, size_t size, size_t order, float* coeff,
+ float* kcoeff, float* err);
Added: trunk/Lib/sandbox/cdavid/src/lpc.tpl
===================================================================
--- trunk/Lib/sandbox/cdavid/src/lpc.tpl 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/src/lpc.tpl 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,80 @@
+[+ AutoGen5 template c +]
+/*
+ * Last Change: Tue Nov 28 03:00 PM 2006 J
+ *
+ * vim:syntax=c
+ */
+#include <stdlib.h> /* for malloc and co */
+#include <stdio.h>
+
+#include "levinson.h"
+#include "autocorr_nofft.h"
+
+#include "lpc.h"
+
+[+ For float_type +]
+/*
+ * ([+ (get "type_name") +] version) Compute lpc coeff (order x coefficients) of a
+ * contiguous array
+ *
+ * err is a double, coeff must be able to contain order+1 elements, and kcoeff
+ * order elements
+ */
+int [+ (get "short_name") +]_lpc(const [+ (get "type_name") +]* signal,
+ size_t size, size_t order, [+ (get "type_name") +]* coeff,
+ [+ (get "type_name") +]* kcoeff, [+ (get "type_name") +]* err)
+{
+ size_t i, nbuff, ncache;
+ [+ (get "type_name") +] *buff, *cache, biasnorm;
+ int status;
+
+ biasnorm = 1.0/size;
+ nbuff = order + 1;
+ ncache = order;
+
+ buff = malloc(sizeof(*buff) * nbuff);
+ if (buff == NULL) {
+ status = -2;
+ goto fail_buff_malloc;
+ }
+
+ cache = malloc(sizeof(*cache) * ncache);
+ if (cache == NULL) {
+ status = -2;
+ goto fail_cache_malloc;
+ }
+
+ /*
+ * Compute the autocorreleation up to lag order, normalized by the
+ * size of the signal
+ */
+ [+ (get "short_name") +]_xcorr_nofft_1d(signal, size, buff, order);
+ for(i = 0; i < nbuff; ++i) {
+ buff[i] *= biasnorm;
+ }
+
+ /*
+ * Compute the inverse coefficients using (simple) levinson recursive algo
+ */
+ status = [+ (get "short_name") +]_levinson1d(buff, order,
+ coeff, err, kcoeff, cache);
+ if (status) {
+ status = -1;
+ goto fail_levinson;
+ }
+
+ free(cache);
+ free(buff);
+
+ return 0;
+
+fail_levinson:
+ free(cache);
+fail_cache_malloc:
+ free(buff);
+fail_buff_malloc:
+ fprintf(stderr, "Failure\n");
+ return status;
+}
+
+[+ ENDFOR float_type +]
Added: trunk/Lib/sandbox/cdavid/tests/test_autocorr.py
===================================================================
--- trunk/Lib/sandbox/cdavid/tests/test_autocorr.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/tests/test_autocorr.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,385 @@
+#! /usr/bin/env python
+# Last Change: Mon Nov 27 08:00 PM 2006 J
+
+from numpy.testing import *
+from numpy.random import randn, seed
+from numpy import correlate, array, concatenate, require
+
+from numpy.ctypeslib import ndpointer, load_library
+from ctypes import c_uint
+
+set_package_path()
+from autocorr import _raw_autocorr_1d, _raw_autocorr_1d_noncontiguous
+from autocorr import autocorr_oneside_nofft as autocorr
+from autocorr import _autocorr_oneside_nofft_py as autocorr_py
+restore_path()
+
+import numpy
+
+# number of decimals to check
+nd = 20
+# minimum number of correct decimals required
+md = 12
+
+a = array([1, 2, 3.])
+b = a + 3
+
+x = concatenate((a, b)).reshape(2, 3)
+
+# float and double C order
+xc = require(x, dtype = numpy.float64, requirements = 'C')
+xcf = require(x, dtype = numpy.float32, requirements = 'C')
+xc1 = xc[0]
+xcf1 = xcf[0]
+
+# float and double F order
+xf = require(x, dtype = numpy.float64, requirements = 'FORTRAN')
+xff = require(x, dtype = numpy.float32, requirements = 'FORTRAN')
+xf1 = xf[0]
+xff1 = xff[0]
+
+# This class tests the C functions directly. This is more a debugging tool
+# that a test case, as the tested functions are not part of the public API
+class test_ctype_1d(NumpyTestCase):
+ def check_contiguous_double(self):
+ # double test
+ xt = xc1
+ yt = _raw_autocorr_1d(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_contiguous_float(self):
+ # float test
+ xt = xcf1
+
+ yt = _raw_autocorr_1d(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_non_contiguous_double(self):
+ # double test
+ xt = xf1
+ yt = _raw_autocorr_1d_noncontiguous(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_non_contiguous_float(self):
+ # float test
+ xt = xff1
+ yt = _raw_autocorr_1d_noncontiguous(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+# Test autocorrelation for rank 1 arrays
+class test_autocorr_1d(NumpyTestCase):
+ def check_contiguous_double(self):
+ # double test
+ xt = xc1
+ yt = autocorr(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_contiguous_float(self):
+ # float test
+ xt = xcf1
+
+ yt = autocorr(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_non_contiguous_double(self):
+ # double test
+ xt = xf1
+ yt = autocorr(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+ def check_non_contiguous_float(self):
+ # float test
+ xt = xff1
+ yt = autocorr(xt, xt.size - 1)
+
+ yr = correlate(xt, xt, mode = 'full')
+ yr = yr[xt.size-1:]
+
+ assert_array_equal(yt, yr)
+
+# This class is a pure python implementation of autocorrelation
+# with rank 2 arrays. This will be used in the above test cases;
+# this function implements the expected behaviour of the public
+# autocorr function.
+class test_autocorr_py(NumpyTestCase):
+ def check_full(self):
+ xt = xc
+ axis = -1
+ lag = xt.shape[axis] - 1
+ yt = autocorr_py(xt, lag, axis = axis)
+
+ yr = yt.copy()
+ for i in range(xt.shape[(axis +1) % 2]):
+ tmp = correlate(xt[i], xt[i], 'full')
+ center = xt[i].size - 1
+ assert_array_equal(tmp[center:center+1+lag], yt[i])
+
+ xt = xc
+ axis = 0
+ lag = xt.shape[axis] - 1
+ yt = autocorr_py(xt, lag, axis = axis)
+
+ yr = yt.copy()
+ for i in range(xt.shape[(axis +1) % 2]):
+ tmp = correlate(xt[:, i], xt[:, i], 'full')
+ center = xt[:,i].size - 1
+ assert_array_equal(tmp[center:center+1+lag], yt[:, i])
+
+ def check_partial(self):
+ xt = xc
+ axis = -1
+ lag = 1
+ yt = autocorr_py(xt, lag, axis = axis)
+
+ yr = yt.copy()
+ for i in range(xt.shape[(axis +1) % 2]):
+ tmp = correlate(xt[i], xt[i], 'full')
+ center = xt[i].size - 1
+ assert_array_equal(tmp[center:center+1+lag], yt[i])
+
+ xt = xc
+ axis = 0
+ lag = 1
+ yt = autocorr_py(xt, lag, axis = axis)
+
+ yr = yt.copy()
+ for i in range(xt.shape[(axis +1) % 2]):
+ tmp = correlate(xt[:, i], xt[:, i], 'full')
+ center = xt[:,i].size - 1
+ assert_array_equal(tmp[center:center+1+lag], yt[:, i])
+
+# Test autocorrelation for rank 2 arrays
+class test_autocorr_2d(NumpyTestCase):
+ def check_double_full(self):
+ # C, axis 1 test
+ xt = xc
+ axis = -1
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # C, axis 0 test
+ xt = xc
+ axis = 0
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 0 test
+ xt = xf
+ axis = 0
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 1 test
+ xt = xf
+ axis = -1
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ def check_float(self):
+ # C, axis 1 test
+ xt = xcf
+ axis = -1
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # C, axis 0 test
+ xt = xcf
+ axis = 0
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 0 test
+ xt = xff
+ axis = 0
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 1 test
+ xt = xff
+ axis = -1
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ def check_double_partial(self):
+ # C, axis 1 test
+ xt = xc
+ axis = -1
+ lag = 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # C, axis 0 test
+ xt = xc
+ axis = 0
+ lag = 0
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 0 test
+ xt = xf
+ axis = 1
+ lag = xt.shape[axis] - 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+ # F, axis 1 test
+ xt = xf
+ axis = -1
+ lag = 1
+ yt = autocorr(xt, lag, axis = axis)
+
+ yr = autocorr_py(xt, lag, axis = axis)
+ assert_array_equal(yt, yr)
+
+if __name__ == "__main__":
+ ScipyTest().run()
+
+#class test_autocorr_2d(NumpyTestCase):
+# def check_double(self):
+# # C, axis 1 test
+# xt = xc
+# axis = -1
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[i], xt[i], 'full')
+# assert_array_equal(tmp[lag:], yt[i])
+#
+# # C, axis 0 test
+# xt = xc
+# axis = 0
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[:, i], xt[:, i], 'full')
+# assert_array_equal(tmp[lag:], yt[:, i])
+#
+# # F, axis 0 test
+# xt = xf
+# axis = 0
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[:, i], xt[:, i], 'full')
+# assert_array_equal(tmp[lag:], yt[:, i])
+#
+# # F, axis 1 test
+# xt = xf
+# axis = -1
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[i], xt[i], 'full')
+# assert_array_equal(tmp[lag:], yt[i])
+#
+# def check_float(self):
+# # C, axis 1 test
+# xt = xcf
+# axis = -1
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[i], xt[i], 'full')
+# assert_array_equal(tmp[lag:], yt[i])
+#
+# # C, axis 0 test
+# xt = xcf
+# axis = 0
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[:, i], xt[:, i], 'full')
+# assert_array_equal(tmp[lag:], yt[:, i])
+#
+# # F, axis 0 test
+# xt = xff
+# axis = 0
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[:, i], xt[:, i], 'full')
+# assert_array_equal(tmp[lag:], yt[:, i])
+#
+# # F, axis 1 test
+# xt = xff
+# axis = -1
+# lag = xt.shape[axis] - 1
+# yt = autocorr(xt, lag, axis = axis)
+#
+# yr = yt.copy()
+# for i in range(xt.shape[(axis +1) % 2]):
+# tmp = correlate(xt[i], xt[i], 'full')
+# assert_array_equal(tmp[lag:], yt[i])
+#
Added: trunk/Lib/sandbox/cdavid/tests/test_lpc.py
===================================================================
--- trunk/Lib/sandbox/cdavid/tests/test_lpc.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/cdavid/tests/test_lpc.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -0,0 +1,187 @@
+#! /usr/bin/env python
+# Last Change: Tue Nov 28 03:00 PM 2006 J
+
+from numpy.testing import *
+from numpy.random import randn, seed
+from numpy import correlate, array, concatenate, require
+
+from numpy.ctypeslib import ndpointer, load_library
+from ctypes import c_uint
+
+set_package_path()
+from lpc import _lpc2_py as lpc_py
+from lpc import lpc_ref, lpc2
+from autocorr import autocorr_oneside_nofft
+restore_path()
+
+import numpy
+
+# number of decimals to check
+nd = 20
+# minimum number of correct decimals required
+md = 12
+
+a = array([1, 2, 3.])
+b = a + 3
+
+x = concatenate((a, b)).reshape(2, 3)
+
+# float and double C order
+xc = require(x, dtype = numpy.float64, requirements = 'C')
+xcf = require(x, dtype = numpy.float32, requirements = 'C')
+xc1 = xc[0]
+xcf1 = xcf[0]
+
+# float and double F order
+xf = require(x, dtype = numpy.float64, requirements = 'FORTRAN')
+xff = require(x, dtype = numpy.float32, requirements = 'FORTRAN')
+xf1 = xf[0]
+xff1 = xff[0]
+
+# This class uses lpc in 1 dimension and loop on the axis. Is tested against
+# a direct matrix inversion of the autocorrelation matrix (using matrix inverse
+# instead of levinson durbin)
+class test_lpc_py(NumpyTestCase):
+ def check_float(self):
+ # Axis -1
+ xt = xcf
+ axis = -1
+ order = 1
+
+ a, k, e = lpc_py(xt, order, axis)
+ assert a.dtype == k.dtype == e.dtype == numpy.float32
+
+ tmp = numpy.zeros((xt.shape[0], order+1), xt.dtype)
+ for i in range(xt.shape[0]):
+ tmp[i] = lpc_ref(xt[i], order)
+
+ assert_array_almost_equal(tmp, a)
+
+ # Axis 0
+ xt = xcf
+ axis = 0
+ order = 1
+
+ a, k, e = lpc_py(xt, order, axis)
+ assert a.dtype == k.dtype == e.dtype == numpy.float32
+
+ tmp = numpy.zeros((order + 1, xt.shape[1]), xt.dtype)
+ for i in range(xt.shape[1]):
+ tmp[:, i] = lpc_ref(xt[:, i], order)
+
+ assert_array_almost_equal(tmp, a)
+
+ def check_double(self):
+ # Axis -1
+ xt = xc
+ axis = -1
+ order = 1
+
+ a, e, k = lpc_py(xt, order, axis)
+ assert a.dtype == k.dtype == e.dtype == numpy.float64
+
+ tmp = numpy.zeros((xt.shape[0], order+1), xt.dtype)
+ for i in range(xt.shape[0]):
+ tmp[i] = lpc_ref(xt[i], order)
+
+ assert_array_almost_equal(tmp, a)
+
+ # Axis 0
+ xt = xc
+ axis = 0
+ order = 1
+
+ a, e, k = lpc_py(xt, order, axis)
+ assert a.dtype == k.dtype == e.dtype == numpy.float64
+
+ tmp = numpy.zeros((order + 1, xt.shape[1]), xt.dtype)
+ for i in range(xt.shape[1]):
+ tmp[:, i] = lpc_ref(xt[:, i], order)
+
+ assert_array_almost_equal(tmp, a)
+
+class test_lpc(NumpyTestCase):
+ def check_float(self):
+ # Axis -1
+ xt = xcf
+ axis = -1
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert a.dtype == e.dtype == k.dtype == numpy.float32
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+ # Axis 0
+ xt = xcf
+ axis = 0
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert a.dtype == e.dtype == k.dtype == numpy.float32
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+ def check_float_rank1(self):
+ # test rank 1
+ xt = xcf[0]
+ axis = 0
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert a.dtype == e.dtype == k.dtype == numpy.float32
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+ def check_double(self):
+ # Axis -1
+ xt = xc
+ axis = -1
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+ # Axis 0
+ xt = xc
+ axis = 0
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+ def check_double_rank1(self):
+ # test rank 1
+ xt = xc[0]
+ axis = 0
+ order = 1
+
+ a, e, k = lpc2(xt, order, axis)
+ at, et, kt = lpc_py(xt, order, axis)
+
+ assert_array_almost_equal(a, at)
+ assert_array_almost_equal(e, et)
+ assert_array_almost_equal(k, kt)
+
+if __name__ == "__main__":
+ ScipyTest().run()
Modified: trunk/Lib/sandbox/setup.py
===================================================================
--- trunk/Lib/sandbox/setup.py 2006-11-26 05:13:41 UTC (rev 2331)
+++ trunk/Lib/sandbox/setup.py 2006-11-28 08:04:47 UTC (rev 2332)
@@ -76,6 +76,9 @@
# Package for Gaussian Mixture Models
#config.add_subpackage('pyem')
+ # David Cournapeau's corner: autocorrelation, lpc, lpc residual
+ config.add_subpackage('cdavid')
+
# New spline package (based on scipy.interpolate)
#config.add_subpackage('spline')
More information about the Scipy-svn
mailing list