Source code for deltasigma._utils

# -*- coding: utf-8 -*-
# _utils.py
# Miscellaneous functions and stdlib wrappers for MATLAB functions
# that do not find a direct replacement in numpy/scipy.
# Copyright 2013 Giuseppe Venturini
# 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.

""" Miscellaneous functions and wrappers for MATLAB functions
 that do not find a direct replacement in numpy/scipy.
"""

import collections
import fractions
from fractions import Fraction as Fr

import numpy as np
from scipy.signal import lti, ss2tf, ss2zpk, zpk2tf

from ._constants import eps
from ._partitionABCD import partitionABCD


[docs]def rat(x, tol): """Rational fraction approximation. Calculate A and B such that: .. math:: x = \\frac{A}{B} + \\epsilon where: .. math:: |\\epsilon| < tol .. note:: A, B are of type 'int' """ return Fr(float(x)).limit_denominator(int(1 / float(tol))).numerator, \ Fr(float(x)).limit_denominator(int(1 / float(tol))).denominator
gcd = fractions.gcd lcm = lambda a, b: int(a * b / float(gcd(a, b))) lcm.__doc__ = """Calculate the Least Common Multiple of ``a`` and ``b``.\n""" class empty(object): """An empty function used to hold attributes""" def __init__(self): pass
[docs]def mfloor(x): """Round ``x`` towards -Inf. This is a MATLAB-compatible floor function, numpy's ``floor()`` behaves differently. If the elements of ``x`` are complex, real and imaginary parts are rounded separately. """ iform = save_input_form(x) x = carray(x) def _mfloor(z): """Base function to generate the ufunc floor""" z = np.real_if_close(z) if np.iscomplex(z): return _mfloor(np.real(z)) + 1j * _mfloor(np.imag(z)) return np.floor(z) if np.sign(z) >= 0 else -np.ceil(-z) _internal = np.frompyfunc(_mfloor, 1, 1) xf = np.array(_internal(x), dtype=x.dtype) return restore_input_form(xf, iform)
def carray(x): """Check that x is an ndarray. If not, try to convert it to ndarray. """ if not isinstance(x, np.ndarray): if not isinstance(x, collections.Iterable): x = np.array((x,)) else: x = np.array(x) elif not len(x.shape): x = x.reshape((1,)) else: pass # nothing to do here return x
[docs]def cplxpair(x, tol=100): """ Sort complex numbers into complex conjugate pairs. This function replaces MATLAB's cplxpair for vectors. """ x = carray(x) x = np.atleast_1d(x.squeeze()) x = x.tolist() x = [np.real_if_close(i, tol) for i in x] xreal = np.array(list(filter(np.isreal, x))) xcomplex = np.array(list(filter(np.iscomplex, x))) xreal = np.sort_complex(xreal) xcomplex = np.sort_complex(xcomplex) xcomplex_ipos = xcomplex[xcomplex.imag > 0.] xcomplex_ineg = xcomplex[xcomplex.imag <= 0.] if len(xcomplex_ipos) != len(xcomplex_ineg): raise ValueError("Complex numbers can't be paired.") res = [] for i, j in zip(xcomplex_ipos, xcomplex_ineg): if not abs(i - np.conj(j)) < tol * eps: raise ValueError("Complex numbers can't be paired.") res += [j, i] return np.hstack((np.array(res), xreal))
def minreal(tf, tol=None): """Remove pole/zero pairs from a transfer function when the two match within the tolerance ``tol``. **Parameters:** tf : supported TF representation or list ``tf`` may be a transfer function (LTI object) or a list of transfer functions, each of them expressed in a supported representation. tol : float, optional The tolerance to be accepted when simplifying pole-zero pairs. It defaults to the system epsilon if unset. **Returns:** tf_simplified : tuple or list of tuples A list of TFs in zpk format or a TF (aain in zpk format), depending whether a single TF or a list of multiple TFs were passed to the function. """ # initially based on python-control # which is in turn based on octave minreal # then modified considerably # recursively handle multiple tfs if not _is_zpk(tf) and not _is_num_den(tf) and not _is_A_B_C_D(tf) \ and (isinstance(tf, list) or isinstance(tf, tuple)): ret = [] for tfi in tf: ret += [minreal(tfi, tol)] return ret # default accuracy sqrt_eps = np.sqrt(eps) if (hasattr(tf, 'inputs') and not tf.inputs == 1) or \ (hasattr(tf, 'outputs') and not tf.outputs == 1): raise TypeError("Only SISO transfer functions can be evaluated.") if hasattr(tf, 'zeros') and hasattr(tf, 'poles'): # LTI objects have poles and zeros, zeros = tf.zeros poles = tf.poles if hasattr(tf, 'k'): k = tf.k elif hasattr(tf, 'gain'): k = tf.gain else: # k = num[0] / den[0] zeros, poles, k = _get_zpk(tf) zeros = carray(zeros) poles = carray(poles) zeros.sort() poles.sort() reducedzeros = [] for z in zeros: t = tol or 1000 * max(eps, abs(z) * sqrt_eps) idx = np.where(abs(poles - z) < t)[0] if len(idx): # cancel this zero against one of the poles # remove the pole and do not add the zero to the new poles = np.delete(poles, idx[0]) else: # no matching pole reducedzeros.append(z) newzeros = carray(reducedzeros) return (newzeros, poles, k) def diagonal_indices(a, offset=0): """The indices to the diagonal of a 2D array ``a`` The indices are those to the main diagonal (if ``offset`` is 0), or to a secondary diagonal, having the specified offset from the main one. The array ``A`` does not need to be square. **Parameters:** a : ndarray The 2D ndarray for which the diagonal indices should be calculated. offset : int, optional The diagonal offset from the main one. Note that the sup-diagonal is at offset +1, the sub-diagonal at offset -1, and so on. Defaults to 0, which corresponds to the main diagonal. **Returns:** xs, ys : tuples The indices in the two coordinates. Thanks to ``numpy``'s advanced slicing, the diagonal may be accessed with ``A[(xs, ys)]``. """ di, dj = np.diag_indices_from(a[:min(a.shape), :min(a.shape)]) if offset > 0: di, dj = zip(*[(i, j) for i, j in zip(di, dj + offset) if 0 <= j < a.shape[1]]) elif offset == 0: pass else: di, dj = zip(*[(i, j) for i, j in zip(di - offset, dj) if 0 <= i < a.shape[0]]) return di, dj
[docs]def circshift(a, shift): """Shift an array circularly. The ``circshift(a, shift)`` function circularly shifts the values in the array ``a`` by ``shift`` elements. **Parameters:** a : ndarray the array to be shifted. Notice that a should have a greater or equal number of dimensions than ``shift`` (``shift`` being a scalar is equal to ``shift`` being a one-dimension array.) shift : int or ndarray-like of int. the N-th element specifies the shift amount for the N-th dimension of the input array ``a``. If an element in ``shift`` is positive, the values of ``a`` are shifted to higher-index rows (ie down) or to higher-index columns (ie to the right). If the element is negative, the values of ``a`` are shifted in the opposite directions, towards lower-index rows (ie up) or to lower-index columns (ie right). If ``shift`` is an integer, the shift happens along axis 0. All dimensions that do not have a corresponding shift value in ``shift`` are left untouched (ie ``shift=(1, 0, 0)`` is equal to ``shift=(1,)``, with the exception that the former will trigger an ``IndexError`` if ``a.ndim < 3``). **Returns:** The shifted array. """ if np.isscalar(shift): shift = [shift] for axis, ashift in enumerate(shift): a = np.roll(a, ashift, axis=axis) return a
def save_input_form(a): """Save the form of `a` so that it can be restored later on Returns: an object representing the form of `a`, to be passed to restore_input_form(a, form) """ if np.isscalar(a): ret = 'scalar' elif isinstance(a, np.ndarray): ret = a.shape elif type(a) == tuple: ret = 'tuple' elif type(a) == list: ret = 'list' else: raise TypeError("Unsupported input %s" % repr(a)) return ret def restore_input_form(a, form): """Restore the form of `a` according to `form`. Returns: the object `a`, in the correct `form`. Note: use `save_input_form(a)` to get the object `form` """ if form == 'scalar': a = np.real_if_close(a) if not np.isscalar(a): a = a.reshape((1, ))[0] elif form == 'tuple': if not type(a) == tuple: a = [np.real_if_close(i).reshape((1,))[0] for i in a] a = tuple(a) elif form == 'list': if not type(a) == list: a = [np.real_if_close(i).reshape((1,))[0] for i in a] a = list(a) else: a = a.reshape(form) return a
[docs]def pretty_lti(arg): """Given the lti object ``arg`` return a *pretty* representation.""" z, p, k = _get_zpk(arg) z = np.atleast_1d(z) p = np.atleast_1d(p) z = np.round(np.real_if_close(z), 4) p = np.round(np.real_if_close(p), 4) k = np.round(np.real_if_close(k), 4) signs = {1:'+', -1:'-'} if not len(z) and not len(p): return "%g" % k ppstr = ["", "", ""] if np.allclose(k, 0., atol=1e-5): return "0" if k != 1: if np.isreal(k): ppstr[1] = "%g " % k else: # quadrature modulators support ppstr[1] += "(%g %s %gj) " % (np.real(k), signs[np.sign(np.imag(k))], np.abs(np.imag(k))) for i, s in zip((0, 2), (z, p)): rz = None m = 1 try: sorted_singularities = cplxpair(s) quadrature = False except ValueError: # quadrature modulator sorted_singularities = np.sort_complex(s) quadrature = True for zindex, zi in enumerate(sorted_singularities): zi = np.round(np.real_if_close(zi), 4) if np.isreal(zi) or quadrature: if len(sorted_singularities) > zindex + 1 and \ sorted_singularities[zindex + 1] == zi: m += 1 continue if zi == 0.: ppstr[i] += "z" elif np.isreal(zi): ppstr[i] += "(z %s %g)" % (signs[np.sign(-zi)], np.abs(zi)) else: ppstr[i] += "(z %s %g %s %gj)" % (signs[np.sign(np.real(-zi))], np.abs(np.real(zi)), signs[np.sign(np.imag(-zi))], np.abs(np.imag(zi))) if m == 1: ppstr[i] += " " else: ppstr[i] += "^%d " % m m = 1 else: if len(sorted_singularities) > zindex + 2 and \ sorted_singularities[zindex + 2] == zi: m += .5 continue if rz is None: rz = zi continue ppstr[i] += "(z^2 %s %gz %s %g)" % \ (signs[np.sign(np.real_if_close(np.round(-rz - zi, 3)))], np.abs(np.real_if_close(np.round(-rz - zi, 3))), signs[np.sign(np.real_if_close(np.round(rz * zi, 4)))], np.abs(np.real_if_close(np.round(rz * zi, 4)))) if m == 1: ppstr[i] += " " else: ppstr[i] += "^%d " % m rz = None m = 1 ppstr[i] = ppstr[i][:-1] if len(ppstr[i]) else "1" if ppstr[2] == '1': return ppstr[1] + ppstr[0] else: if ppstr[0] == '1' and len(ppstr[1]) and float(ppstr[1]) != 1.: ppstr[0] = ppstr[1][:-1] ppstr[1] = "" space_pad_ln = len(ppstr[1]) fraction_line = "-" * (max(len(ppstr[0]), len(ppstr[2])) + 2) ppstr[1] += fraction_line ppstr[0] = " "*space_pad_ln + ppstr[0].center(len(fraction_line)) ppstr[2] = " "*space_pad_ln + ppstr[2].center(len(fraction_line)) return "\n".join(ppstr)
def _get_zpk(arg, input=0): """Utility method to convert the input arg to a z, p, k representation. **Parameters:** arg, which may be: * ZPK tuple, * num, den tuple, * A, B, C, D tuple, * a scipy LTI object, * a sequence of the tuples of any of the above types. input : scalar In case the system has multiple inputs, which input is to be considered. Input `0` means first input, and so on. **Returns:** The sequence of ndarrays z, p and a scalar k **Raises:** TypeError, ValueError .. warn: support for MISO transfer functions is experimental. """ z, p, k = None, None, None if isinstance(arg, np.ndarray): # ABCD matrix A, B, C, D = partitionABCD(arg) z, p, k = ss2zpk(A, B, C, D, input=input) elif isinstance(arg, lti): z, p, k = arg.zeros, arg.poles, arg.gain elif _is_zpk(arg): z, p, k = np.atleast_1d(arg[0]), np.atleast_1d(arg[1]), arg[2] elif _is_num_den(arg): sys = lti(*arg) z, p, k = sys.zeros, sys.poles, sys.gain elif _is_A_B_C_D(arg): z, p, k = ss2zpk(*arg, input=input) elif isinstance(arg, collections.Iterable): ri = 0 for i in arg: # Note we do not check if the user has assembled a list with # mismatched lti representations. if hasattr(i, 'B'): iis = i.B.shape[1] if input < ri + iis: z, p, k = ss2zpk(i.A, i.B, i.C, i.D, input=input - ri) break else: ri += iis continue elif _is_A_B_C_D(arg): iis = arg[1].shape[1] if input < ri + iis: z, p, k = ss2zpk(*arg, input=input - ri) break else: ri += iis continue else: if ri == input: sys = lti(*i) z, p, k = sys.zeros, sys.poles, sys.gain break else: ri += 1 continue ri += 1 if (z, p, k) == (None, None, None): raise ValueError("The LTI representation does not have enough" + "inputs: max %d, looking for input %d" % (ri - 1, input)) else: raise TypeError("Unknown LTI representation: %s" % arg) return z, p, k def _get_num_den(arg, input=0): """Utility method to convert the input arg to a (num, den) representation. **Parameters:** arg, which may be: * ZPK tuple, * num, den tuple, * A, B, C, D tuple, * a scipy LTI object, * a sequence of the tuples of any of the above types. input : scalar In case the system has multiple inputs, which input is to be considered. Input `0` means first input, and so on. **Returns:** The sequence of ndarrays num, den **Raises:** TypeError, ValueError .. warn: support for MISO transfer functions is experimental. """ num, den = None, None if isinstance(arg, np.ndarray): # ABCD matrix A, B, C, D = partitionABCD(arg) num, den = ss2tf(A, B, C, D, input=input) elif isinstance(arg, lti): num, den = arg.num, arg.den elif _is_num_den(arg): num, den = carray(arg[0]).squeeze(), carray(arg[1]).squeeze() elif _is_zpk(arg): num, den = zpk2tf(*arg) elif _is_A_B_C_D(arg): num, den = ss2tf(*arg, input=input) elif isinstance(arg, collections.Iterable): ri = 0 for i in arg: # Note we do not check if the user has assembled a list with # mismatched representations. if hasattr(i, 'B'): # lti iis = i.B.shape[1] if input < ri + iis: num, den = ss2tf(i.A, i.B, i.C, i.D, input=input - ri) break else: ri += iis else: sys = lti(*i) iis = sys.B.shape[1] if input < ri + iis: num, den = ss2tf( sys.A, sys.B, sys.C, sys.D, input=input - ri) break else: ri += iis if (num, den) == (None, None): raise ValueError("The LTI representation does not have enough" + "inputs: max %d, looking for input %d" % (ri - 1, input)) else: raise TypeError("Unknown LTI representation: %s" % arg) if len(num.shape) > 1: num = num.squeeze() if len(den.shape) > 1: den = den.squeeze() # default accuracy: sqrt_ps sqrt_eps = np.sqrt(eps) while len(num.shape) and len(num): if abs(num[0]) < sqrt_eps: num = num[1:] else: break while len(den.shape) and len(den): if abs(den[0]) < sqrt_eps: den = den[1:] else: break den = np.atleast_1d(den) num = np.atleast_1d(num) return num, den def _getABCD(arg): """Utility method to convert the input arg to an A, B, C, D representation. **Parameters:** arg, which may be: * ZPK tuple, * num, den tuple, * A, B, C, D tuple, * a scipy LTI object, * a sequence of the tuples of any of the above types. **Returns:** The sequence of ndarrays A, B, C, D **Raises:** TypeError, ValueError """ if isinstance(arg, np.ndarray): # ABCD matrix A, B, C, D = partitionABCD(arg) elif isinstance(arg, lti): A, B, C, D = arg.A, arg.B, arg.C, np.atleast_2d(arg.D) elif _is_zpk(arg) or _is_num_den(arg) or _is_A_B_C_D(arg): sys = lti(*arg) A, B, C, D = sys.A, sys.B, sys.C, sys.D elif isinstance(arg, collections.Iterable): A, B, C, D = None, None, None, None for i in arg: # Note we do not check if the user has assembled a list with # mismatched lti representations. sys = lti(*i) if not hasattr(i, 'A') else i if A is None: A = sys.A elif not np.allclose(sys.A, A, atol=1e-8, rtol=1e-5): raise ValueError("Mismatched lti list, A matrix disagreement.") else: pass if B is None: B = sys.B else: B = np.hstack((B, sys.B)) if C is None: C = sys.C elif not np.allclose(sys.C, C, atol=1e-8, rtol=1e-5): raise ValueError("Mismatched lti list, C matrix disagreement.") if D is None: D = sys.D else: D = np.hstack((D, sys.D)) else: raise TypeError("Unknown LTI representation: %s" % arg) return A, B, C, D def _is_zpk(arg): """Can the argument be safely assumed to be a zpk tuple?""" return isinstance(arg, collections.Iterable) and len(arg) == 3 and \ isinstance(arg[0], collections.Iterable) and \ isinstance(arg[1], collections.Iterable) and np.isscalar(arg[2]) def _is_num_den(arg): """Can the argument be safely assumed to be a num, den tuple?""" return isinstance(arg, collections.Iterable) and len(arg) == 2 and \ isinstance(arg[0], collections.Iterable) and \ isinstance(arg[1], collections.Iterable) def _is_A_B_C_D(arg): """Can the argument be safely assumed to be an (A, B, C, D) tuple?""" return isinstance(arg, collections.Iterable) and len(arg) == 4 and \ (isinstance(arg[0], collections.Iterable) or np.is_scalar(arg[0])) and \ (isinstance(arg[1], collections.Iterable) or np.is_scalar(arg[1])) and \ (isinstance(arg[2], collections.Iterable) or np.is_scalar(arg[2])) and \ (isinstance(arg[3], collections.Iterable) or np.is_scalar(arg[3]))
[docs]def mround(x): """Round ``x`` to the nearest integers. This is a MATLAB-compatible round function, numpy's ``round()`` behaves differently. Behaviour with a fractional part of 0.5: * Positive elements with a fractional part of 0.5 round up to the nearest positive integer. * Negative elements with a fractional part of -0.5 round down to the nearest negative integer. If the elements of ``x`` are complex, real and imaginary parts are rounded separately. **Example:** >>> mround([-1.9, -0.5, -0.2, 3.4, 4.5, 5.6, 7.0, 2.4+3.6j]) [-2.0, -1.0, 0.0, 3.0, 5.0, 6.0, 7.0, 2.0+4.0j] """ iform = save_input_form(x) x = carray(x) def _mround(z): """Base function to generate the ufunc round""" z = np.real_if_close(z) if np.iscomplex(z): return _mround(np.real(z)) + 1j*_mround(np.imag(z)) s = 1 if z >= 0 else -1 res = z - s*(abs(z) % 1) if abs(z) % 1 < .5 \ else z + s*(1 - (abs(z)%1)) return res _internal = np.frompyfunc(_mround, 1, 1) xf = np.array(_internal(x), dtype=x.dtype) return restore_input_form(xf, iform)