[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 02:05:13 CST 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