From 475b2f72eeb961d53043f8808ae53fd1e15a4527 Mon Sep 17 00:00:00 2001 From: Shayne Hodge Date: Wed, 13 Aug 2014 17:40:53 -0700 Subject: [PATCH 1/3] Add timing test script and change imports to work. --- deltasigma/_simulateDSM_numba.py | 2 +- deltasigma/_simulateDSM_python.py | 22 +++++------ deltasigma/_utils.py | 4 +- deltasigma/numba_timing_tests.py | 63 +++++++++++++++++++++++++++++++ 4 files changed, 77 insertions(+), 14 deletions(-) create mode 100644 deltasigma/numba_timing_tests.py diff --git a/deltasigma/_simulateDSM_numba.py b/deltasigma/_simulateDSM_numba.py index deada0f..398c225 100644 --- a/deltasigma/_simulateDSM_numba.py +++ b/deltasigma/_simulateDSM_numba.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 from numba import double, int32, complex64 from numba.types import pyobject diff --git a/deltasigma/_simulateDSM_python.py b/deltasigma/_simulateDSM_python.py index 6dd20dc..22bb947 100644 --- a/deltasigma/_simulateDSM_python.py +++ b/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,7 +181,7 @@ 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. @@ -189,7 +189,7 @@ def ds_quantize(y, n): 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 diff --git a/deltasigma/_utils.py b/deltasigma/_utils.py index 3ccad42..28a063b 100644 --- a/deltasigma/_utils.py +++ b/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): diff --git a/deltasigma/numba_timing_tests.py b/deltasigma/numba_timing_tests.py new file mode 100644 index 0000000..8781eef --- /dev/null +++ b/deltasigma/numba_timing_tests.py @@ -0,0 +1,63 @@ +# -*- 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 +import timeit +from deltasigma import synthesizeNTF, realizeNTF, stuffABCD + +# setup variables, wrap in functino for timeit + +def setup_env(): + 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 u, H, ABCD + +def test1(): + import _simulateDSM_numba + u, H, ABCD = setup_env() + _simulateDSM_numba.simulateDSM(u, H, 2, 0) + +def test2(): + import _simulateDSM_numba + u, H, ABCD = setup_env() + _simulateDSM_numba.simulateDSM(u, ABCD, 2, 0) + +def test3(): + import _simulateDSM_python + u, H, ABCD = setup_env() + _simulateDSM_python.simulateDSM(u, H, 2, 0) + +def test4(): + import _simulateDSM_python + u, H, ABCD = setup_env() + _simulateDSM_python.simulateDSM(u, ABCD, 2, 0) + +if __name__ == '__main__': + print(timeit.timeit("test1()", setup="from __main__ import test1", + number=1)) + print(timeit.timeit("test2()", setup="from __main__ import test2", + number=1)) + print(timeit.timeit("test3()", setup="from __main__ import test3", + number=1)) + print(timeit.timeit("test4()", setup="from __main__ import test4", + number=1)) From f5b7e79ae63df94a17ea7db2103813124cfc4760 Mon Sep 17 00:00:00 2001 From: Shayne Hodge Date: Wed, 13 Aug 2014 18:19:56 -0700 Subject: [PATCH 2/3] Change jit to autojit, add profiling script. --- deltasigma/_simulateDSM_numba.py | 83 ++---------------------------- deltasigma/numba_profile_script.py | 44 ++++++++++++++++ 2 files changed, 48 insertions(+), 79 deletions(-) create mode 100644 deltasigma/numba_profile_script.py diff --git a/deltasigma/_simulateDSM_numba.py b/deltasigma/_simulateDSM_numba.py index 398c225..f744df9 100644 --- a/deltasigma/_simulateDSM_numba.py +++ b/deltasigma/_simulateDSM_numba.py @@ -28,7 +28,7 @@ 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 @@ -38,85 +38,8 @@ # x0 = 1D double array #@jit(arg_types=[double[:, :], complex64[:], int32[:], double[:]]) -@jit() +@autojit() 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'); nlev = carray(nlev) u = np.array(u) if not hasattr(u, 'ndim') else u @@ -193,6 +116,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: diff --git a/deltasigma/numba_profile_script.py b/deltasigma/numba_profile_script.py new file mode 100644 index 0000000..aa3d915 --- /dev/null +++ b/deltasigma/numba_profile_script.py @@ -0,0 +1,44 @@ +# -*- 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 cProfile +import pstats + +def test_script(): + 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') + + _simulateDSM_numba.simulateDSM(u, H, 2, 0) + + #_simulateDSM_numba.simulateDSM(u, ABCD, 2, 0) + + #_simulateDSM_python.simulateDSM(u, H, 2, 0) + + #_simulateDSM_python.simulateDSM(u, ABCD, 2, 0) + +cProfile.run('test_script()', 'numba_stats') +p = pstats.Stats('numba_stats') +p.sort_stats('time').print_stats(25) From 1629db0856c111bd8fda1eb41e31fedbee85433d Mon Sep 17 00:00:00 2001 From: Shayne Hodge Date: Fri, 15 Aug 2014 12:58:19 -0700 Subject: [PATCH 3/3] Latest (failed) attempts to speed up. Only using ABCD form. --- deltasigma/_simulateDSM_numba.py | 93 +++++++++------------- deltasigma/numba_profile_script.py | 37 ++++++--- deltasigma/tests/test_simulateDSM_numba.py | 14 ++-- 3 files changed, 68 insertions(+), 76 deletions(-) diff --git a/deltasigma/_simulateDSM_numba.py b/deltasigma/_simulateDSM_numba.py index f744df9..6ff6939 100644 --- a/deltasigma/_simulateDSM_numba.py +++ b/deltasigma/_simulateDSM_numba.py @@ -22,9 +22,7 @@ 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 @@ -37,66 +35,47 @@ # nlev = 1D int array # x0 = 1D double array -#@jit(arg_types=[double[:, :], complex64[:], int32[:], double[:]]) -@autojit() -def simulateDSM(u, arg2, nlev=2, x0=0): +# 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 diff --git a/deltasigma/numba_profile_script.py b/deltasigma/numba_profile_script.py index aa3d915..302dde9 100644 --- a/deltasigma/numba_profile_script.py +++ b/deltasigma/numba_profile_script.py @@ -19,10 +19,11 @@ import numpy as np from deltasigma import synthesizeNTF, realizeNTF, stuffABCD import _simulateDSM_numba +import _simulateDSM_python import cProfile import pstats -def test_script(): +def setup(): OSR = 32 H = synthesizeNTF(5, OSR, 1) N = 8192 @@ -30,15 +31,25 @@ def test_script(): 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') - - _simulateDSM_numba.simulateDSM(u, H, 2, 0) - - #_simulateDSM_numba.simulateDSM(u, ABCD, 2, 0) - - #_simulateDSM_python.simulateDSM(u, H, 2, 0) - - #_simulateDSM_python.simulateDSM(u, ABCD, 2, 0) - -cProfile.run('test_script()', 'numba_stats') -p = pstats.Stats('numba_stats') -p.sort_stats('time').print_stats(25) + 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) diff --git a/deltasigma/tests/test_simulateDSM_numba.py b/deltasigma/tests/test_simulateDSM_numba.py index 80923b8..63dafe8 100644 --- a/deltasigma/tests/test_simulateDSM_numba.py +++ b/deltasigma/tests/test_simulateDSM_numba.py @@ -43,17 +43,19 @@ def setUp(self): a, g, b, c = realizeNTF(self.H, 'CRFB') self.ABCD = stuffABCD(a, g, b, c, form='CRFB') - def test_simulateDSM_numba(self): + def test_simulateDSM_numba_ABCD(self): """Test function for simulateDSM_numba()""" #print self.ABCD #print self.u #print self.H - v, xn, xmax, y = ds._simulateDSM._simulateDSM_numba( - self.u, self.H, 2, 0) - self.assertTrue(np.allclose( - v.reshape(-1), self.v_ref.reshape(-1), atol=1e-6, rtol=1e-4)) - self.assertTrue(np.allclose(y, self.y_ref, atol=1e-6, rtol=1e-4)) v, xn, xmax, y = ds._simulateDSM._simulateDSM_numba( self.u, self.ABCD, 2, 0) self.assertTrue(np.allclose(v, self.v_ref, atol=1e-6, rtol=1e-4)) self.assertTrue(np.allclose(y, self.y_ref, atol=1e-6, rtol=1e-4)) + +# def test_simulateDSM_numba_H(self): +# v, xn, xmax, y = ds._simulateDSM._simulateDSM_numba( +# self.u, self.H, 2, 0) +# self.assertTrue(np.allclose( +# v.reshape(-1), self.v_ref.reshape(-1), atol=1e-6, rtol=1e-4)) +# self.assertTrue(np.allclose(y, self.y_ref, atol=1e-6, rtol=1e-4))