Add timing test script and change imports to work.

deltasigma/_simulateDSM_numba.py

@@ -22,13 +22,11 @@
 
 import numpy as np
 
-from scipy.signal import tf2zpk, zpk2ss
-from scipy.linalg import orth, norm, inv
-from ._utils import carray, _get_zpk
+from _utils import carray
 
 from numba import double, int32, complex64
 from numba.types import pyobject
-from numba.decorators import jit
+from numba.decorators import jit, autojit
 
 # Using object mode, is nopython mode faster?  Array slicing doesn't work in it
 # Input data type assumptions
@@ -37,143 +35,47 @@
 # nlev = 1D int array
 # x0 = 1D double array
 
-#@jit(arg_types=[double[:, :], complex64[:], int32[:], double[:]])
-@jit()
-def simulateDSM(u, arg2, nlev=2, x0=0):
-    """Simulate a Delta Sigma modulator
-
-    **Syntax:**
-
-     * [v, xn, xmax, y] = simulateDSM(u, ABCD, nlev=2, x0=0)
-     * [v, xn, xmax, y] = simulateDSM(u, ntf, nlev=2, x0=0)
-
-    Compute the output of a general delta-sigma modulator with input u,
-    a structure described by ABCD, an initial state x0 (default zero) and
-    a quantizer with a number of levels specified by nlev.
-    Multiple quantizers are implied by making nlev an array,
-    and multiple inputs are implied by the number of rows in u.
-
-    Alternatively, the modulator may be described by an NTF.
-    The NTF is zpk object. (And the STF is assumed to be 1.)
-    The structure that is simulated is the block-diagonal structure used by
-    ``zpk2ss()``.
-
-    **Example:**
-
-    Simulate a 5th-order binary modulator with a half-scale sine-wave input and
-    plot its output in the time and frequency domains.::
-
-        import numpy as np
-        from deltasigma import *
-        OSR = 32
-        H = synthesizeNTF(5, OSR, 1)
-        N = 8192
-        fB = np.ceil(N/(2*OSR))
-        f = 85
-        u = 0.5*np.sin(2*np.pi*f/N*np.arange(N))
-        v = simulateDSM(u, H)[0]
-
-    Graphical display of the results:
-
-    .. plot::
-
-        import numpy as np
-        import pylab as plt
-        from numpy.fft import fft
-        from deltasigma import *
-        OSR = 32
-        H = synthesizeNTF(5, OSR, 1)
-        N = 8192
-        fB = np.ceil(N/(2*OSR))
-        f = 85
-        u = 0.5*np.sin(2*np.pi*f/N*np.arange(N))
-        v = simulateDSM(u, H)[0]
-        plt.figure(figsize=(10, 7))
-        plt.subplot(2, 1, 1)
-        t = np.arange(85)
-        # the equivalent of MATLAB 'stairs' is step in matplotlib
-        plt.step(t, u[t], 'g', label='u(n)')
-        plt.hold(True)
-        plt.step(t, v[t], 'b', label='v(n)')
-        plt.axis([0, 85, -1.2, 1.2]);
-        plt.ylabel('u, v');
-        plt.xlabel('sample')
-        plt.legend()
-        plt.subplot(2, 1, 2)
-        spec = fft(v*ds_hann(N))/(N/4)
-        plt.plot(np.linspace(0, 0.5, N/2 + 1), dbv(spec[:N/2 + 1]))
-        plt.axis([0, 0.5, -120, 0])
-        plt.grid(True)
-        plt.ylabel('dBFS/NBW')
-        snr = calculateSNR(spec[:fB], f)
-        s = 'SNR = %4.1fdB' % snr
-        plt.text(0.25, -90, s)
-        s =  'NBW = %7.5f' % (1.5/N)
-        plt.text(0.25, -110, s)
-        plt.xlabel("frequency $1 \\\\rightarrow f_s$")
-
-    Click on "Source" above to see the source code.
-    """
-
-    #fprintf(1,'Warning: You are running the non-mex version of simulateDSM.\n');
-    #fprintf(1,'Please compile the mex version with "mex simulateDSM.c"\n');
+# Temporarily rewriting to work only the ABCD form
 
-    nlev = carray(nlev)
-    u = np.array(u) if not hasattr(u, 'ndim') else u
-    if not max(u.shape) == np.prod(u.shape):
-        warn("Multiple input delta sigma structures have had little testing.")
-    if u.ndim == 1:
-        u = u.reshape((1, -1))
-    nu = u.shape[0]
-    nq = 1 if np.isscalar(nlev) else nlev.shape[0]
-    # extract poles and zeros
+def simulateDSM(u, arg2, nlev=2, x0=0):
+    if isinstance(arg2, np.ndarray):
+        nlev = carray(nlev)
+        u = np.array(u) if not hasattr(u, 'ndim') else u
+        if not max(u.shape) == np.prod(u.shape):
+            warn("Multiple input delta sigma structures have had little testing.")
+        if u.ndim == 1:
+            u = u.reshape((1, -1))
+        ABCD = np.asarray(arg2, dtype=np.float64)
+        nq = 1 if np.isscalar(nlev) else nlev.shape[0]
+        if ABCD.shape[1] != ABCD.shape[0] + u.shape[0]:
+            raise ValueError('The ABCD argument does not have proper dimensions.')
+        if not isinstance(x0, collections.Iterable):
+            x0 = x0*np.ones((ABCD.shape[0] - nq, 1))
+        else:
+            x0 = np.array(x0).reshape((-1, 1))
+        print('u = {}'.format(u))
+        print('ABCD = {}'.format(ABCD))
+        print('nq = {}'.format(nq))
+        print('nlev = {}'.format(nlev))
+        print('x0 = {}'.format(x0))
+        v, xn, xmax, y = simulateDSM_ABCD(u, ABCD, nq, nlev, x0)
+        return v, xn, xmax, y
+    else:
+        raise ValueError("Stop trying to do things that aren't supported yet!")
 
+@jit(argtypes=([double[:, :], double[:, :], int32, int32[:], double[:]]))
+#@autojit
+def simulateDSM_ABCD(u, ABCD, nq, nlev=2, x0=0):
     # Removing type checking that is conflicting with Numba
     #if (hasattr(arg2, 'inputs') and not arg2.inputs == 1) or \
     #   (hasattr(arg2, 'outputs') and not arg2.outputs == 1):
     #        raise TypeError("The supplied TF isn't a SISO transfer function.")
-    if isinstance(arg2, np.ndarray):
-        ABCD = np.asarray(arg2, dtype=np.float64)
-    #    if ABCD.shape[1] != ABCD.shape[0] + nu:
-    #        raise ValueError('The ABCD argument does not have proper dimensions.')
-        form = 1
-    else:
-        zeros, poles, k = _get_zpk(arg2)
-        form = 2
-    #raise TypeError('%s: Unknown transfer function %s' % (__name__, str(arg2)))
-
-    # need to set order and form now.
-    order = carray(zeros).shape[0] if form == 2 else ABCD.shape[0] - nq
-
-    if not isinstance(x0, collections.Iterable):
-        x0 = x0*np.ones((order, 1))
-    else:
-        x0 = np.array(x0).reshape((-1, 1))
-
-    if form == 1:
-        A = ABCD[:order, :order]
-        B = ABCD[:order, order:order+nu+nq]
-        C = ABCD[order:order+nq, :order]
-        D1 = ABCD[order:order+nq, order:order+nu]
-    else:
-        A, B2, C, D2 = zpk2ss(poles, zeros, -1)    # A realization of 1/H
-        # Transform the realization so that C = [1 0 0 ...]
-        C, D2 = np.real_if_close(C), np.real_if_close(D2)
-        Sinv = orth(np.hstack((np.transpose(C), np.eye(order)))) / norm(C)
-        S = inv(Sinv)
-        C = np.dot(C, Sinv)
-        if C[0, 0] < 0:
-            S = -S
-            Sinv = -Sinv
-        A = np.dot(np.dot(S, A), Sinv)
-        B2 = np.dot(S, B2)
-        C = np.hstack((np.ones((1, 1)), np.zeros((1, order-1)))) # C=C*Sinv;
-        D2 = np.zeros((0,))
-        # !!!! Assume stf=1
-        B1 = -B2
-        D1 = 1
-        B = np.hstack((B1, B2))
-
+    nu = u.shape[0]
+    order = ABCD.shape[0] - nq
+    A = ABCD[:order, :order]
+    B = ABCD[:order, order:order+nu+nq]
+    C = ABCD[order:order+nq, :order]
+    D1 = ABCD[order:order+nq, order:order+nu]
     N = u.shape[1]
     v = np.empty((nq, N), dtype=np.float64)
     y = np.empty((nq, N), dtype=np.float64)     # to store the quantizer input
@@ -193,6 +95,8 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
 
     return v.squeeze(), xn.squeeze(), xmax, y.squeeze()
 
+#@jit(arg_types=[double[:, :], int32[:, :]])
+@autojit
 def ds_quantize(y, n):
     """v = ds_quantize(y,n)
     Quantize y to:

deltasigma/_simulateDSM_python.py

@@ -24,7 +24,7 @@
 
 from scipy.signal import tf2zpk, zpk2ss
 from scipy.linalg import orth, norm, inv
-from ._utils import carray, _get_zpk
+from _utils import carray, _get_zpk
 
 def simulateDSM(u, arg2, nlev=2, x0=0):
     """Simulate a Delta Sigma modulator
@@ -54,7 +54,7 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
         from deltasigma import *
         OSR = 32
         H = synthesizeNTF(5, OSR, 1)
-        N = 8192 
+        N = 8192
         fB = np.ceil(N/(2*OSR))
         f = 85
         u = 0.5*np.sin(2*np.pi*f/N*np.arange(N))
@@ -70,7 +70,7 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
         from deltasigma import *
         OSR = 32
         H = synthesizeNTF(5, OSR, 1)
-        N = 8192 
+        N = 8192
         fB = np.ceil(N/(2*OSR))
         f = 85
         u = 0.5*np.sin(2*np.pi*f/N*np.arange(N))
@@ -126,15 +126,15 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
         zeros, poles, k = _get_zpk(arg2)
         form = 2
     #raise TypeError('%s: Unknown transfer function %s' % (__name__, str(arg2)))
-        
+
     # need to set order and form now.
     order = carray(zeros).shape[0] if form == 2 else ABCD.shape[0] - nq
-    
+
     if not isinstance(x0, collections.Iterable):
         x0 = x0*np.ones((order, 1))
     else:
         x0 = np.array(x0).reshape((-1, 1))
-    
+
     if form == 1:
         A = ABCD[:order, :order]
         B = ABCD[:order, order:order+nu+nq]
@@ -150,9 +150,9 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
         if C[0, 0] < 0:
             S = -S
             Sinv = -Sinv
-        A = np.dot(np.dot(S, A), Sinv) 
-        B2 = np.dot(S, B2) 
-        C = np.hstack((np.ones((1, 1)), np.zeros((1, order-1)))) # C=C*Sinv; 
+        A = np.dot(np.dot(S, A), Sinv)
+        B2 = np.dot(S, B2)
+        C = np.hstack((np.ones((1, 1)), np.zeros((1, order-1)))) # C=C*Sinv;
         D2 = np.zeros((0,))
         # !!!! Assume stf=1
         B1 = -B2
@@ -181,15 +181,15 @@ def simulateDSM(u, arg2, nlev=2, x0=0):
 def ds_quantize(y, n):
     """v = ds_quantize(y,n)
     Quantize y to:
-     
+
     * an odd integer in [-n+1, n-1], if n is even, or
     * an even integer in [-n, n], if n is odd.
 
     This definition gives the same step height for both mid-rise
     and mid-tread quantizers.
     """
     v = np.zeros(y.shape)
-    for qi in range(n.shape[0]): 
+    for qi in range(n.shape[0]):
         if n[qi] % 2 == 0: # mid-rise quantizer
             v[qi, 0] = 2*np.floor(0.5*y[qi, 0]) + 1
         else: # mid-tread quantizer

deltasigma/_utils.py

@@ -24,8 +24,8 @@
 import numpy as np
 from scipy.signal import lti, ss2tf, ss2zpk, zpk2tf
 
-from ._constants import eps
-from ._partitionABCD import partitionABCD
+from _constants import eps
+from _partitionABCD import partitionABCD
 
 
 def rat(x, tol):

deltasigma/numba_profile_script.py

@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+# numba_timing_tests.py
+# This script is to help time changes to _simulateDSM_numba.py
+# Copyright 2014 Giuseppe Venturini & Shayne Hodge
+# This file is part of python-deltasigma.
+#
+# python-deltasigma is a 1:1 Python replacement of Richard Schreier's
+# MATLAB delta sigma toolbox (aka "delsigma"), upon which it is heavily based.
+# The delta sigma toolbox is (c) 2009, Richard Schreier.
+#
+# python-deltasigma is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# LICENSE file for the licensing terms.
+
+"""This script is to help time changes to _simulateDSM_numba.py
+"""
+
+import numpy as np
+from deltasigma import synthesizeNTF, realizeNTF, stuffABCD
+import _simulateDSM_numba
+import _simulateDSM_python
+import cProfile
+import pstats
+
+def setup():
+    OSR = 32
+    H = synthesizeNTF(5, OSR, 1)
+    N = 8192
+    f = 85
+    u = 0.5*np.sin(2*np.pi*f/N*np.arange(N))
+    a, g, b, c = realizeNTF(H, 'CRFB')
+    ABCD = stuffABCD(a, g, b, c, form='CRFB')
+    return OSR, H, N, f, u, ABCD
+
+def test_script_numba():
+    OSR, H, N, f, u, ABCD = setup()
+    #_simulateDSM_numba.simulateDSM(u, H, 2, 0)
+    _simulateDSM_numba.simulateDSM(u, ABCD, 2, 0)
+
+def test_script_cpython():
+    OSR, H, N, f, u, ABCD = setup()
+    #_simulateDSM_python.simulateDSM(u, H, 2, 0)
+    _simulateDSM_python.simulateDSM(u, ABCD, 2, 0)
+
+cProfile.run('test_script_numba()', 'numba_stats')
+pn = pstats.Stats('numba_stats')
+print('Numba stats')
+pn.sort_stats('time').print_stats(10)
+pn.sort_stats('cumtime').print_stats(10)
+
+cProfile.run('test_script_cpython()', 'cpython_stats')
+pc = pstats.Stats('cpython_stats')
+print('CPython stats')
+pc.sort_stats('time').print_stats(10)