Source code for deltasigma._stuffABCD

# -*- coding: utf-8 -*-
# _stuffABCD.py
# Module providing the stuffABCD function
# 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.

"""Module providing the stuffABCD() function
"""

from __future__ import division, print_function

import numpy as np

from ._partitionABCD import partitionABCD
from ._utils import carray, diagonal_indices


[docs]def stuffABCD(a, g, b, c, form='CRFB'): """Calculate the ABCD matrix from the parameters of a modulator topology. **Parameters:** a : array_like Feedback/feedforward coefficients from/to the quantizer. Length :math:`n`. g : array_like Resonator coefficients. Length :math:`floor(n/2)`. b : array_like Feed-in coefficients from the modulator input to each integrator. Length :math:`n + 1`. c : array_like Integrator inter-stage coefficients. Length :math:`n`. form : str, optional See :func:`realizeNTF` for a list of supported structures. **Returns:** ABCD : ndarray A state-space description of the modulator loop filter. .. seealso:: :func:`mapABCD`, the inverse function. """ # Code common to all structures. # a, g, b, and c are internally row vectors a = carray(a).reshape((1, -1)) g = carray(g).reshape((1, -1)) b = carray(b).reshape((1, -1)) c = carray(c).reshape((1, -1)) order = max(a.shape) odd = order % 2 even = 1 - odd # OMG ABCD = np.zeros((order + 1, order + 2)) if np.isscalar(b) or max(b.shape) == 1: b = np.atleast_2d(b) b = np.hstack((b, np.zeros((1, order)))) # mutually exclusive matrix stuffing if form == 'CRFB': # C=(0 0...c_n) # This is done as part of the construction of A, below # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b # B2 = -(a_1 a_2... a_n) ABCD[:order, order + 1] = -a diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order,)) subdiag = [i[even:order:2] for i in diagonal_indices(ABCD, -1)] ABCD[subdiag] = c[0, even:order:2] if order > odd: supdiag = [i[odd:order:2] for i in diagonal_indices(ABCD[:order, :order], +1)] ABCD[supdiag] = -g.reshape((-1,)) # row numbers of delaying integrators dly = np.arange(odd + 1, order, 2) ABCD[dly, :] = ABCD[dly, :] + np.dot(np.diag(c[0, dly - 1]), ABCD[dly - 1, :]) elif form == 'CRFF': # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b # B2 = -(c_1 0... 0) ABCD[0, order + 1] = -c[0, 0] # number of elements = order diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order, )) # subdiag = diagonal[:order - 1:2] + 1 subdiag = [i[:order - 1:2] for i in diagonal_indices(ABCD, -1)] ABCD[subdiag] = c[0, 1:order:2] if even: # rows to have g*(following row) subtracted. multg = np.arange(0, order, 2) ABCD[multg, :] = ABCD[multg, :] \ - np.dot(np.diag(g.reshape((-1,))), ABCD[multg + 1, :]) else: if order >= 3: # supdiag = diagonal[2:order:2] - 1 supdiag = [i[1:order:2] for i in diagonal_indices(ABCD[:order, :order], +1)] ABCD[supdiag] = -g.reshape((-1,)) # rows to have c*(preceding row) added. multc = np.arange(2, order, 2) ABCD[multc, :] = ABCD[multc, :] + np.dot(np.diag(c[0, multc]), ABCD[multc - 1, :]) ABCD[order, :order:2] = a[0, :order:2] for i in range(1, order, 2): ABCD[order, :] = ABCD[order, :] + a[0, i]*ABCD[i, :] elif form == 'CIFB': # C=(0 0...c_n) # This is done as part of the construction of A, below # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b # B2 = -(a_1 a_2... a_n) ABCD[:order, order + 1] = -a diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order,)) subdiag = diagonal_indices(ABCD, -1) ABCD[subdiag] = c.reshape((-1,)) supdiag = [i[odd:order+odd:2] for i in diagonal_indices(ABCD[:order, :order], +1)] ABCD[supdiag] = -g.reshape((-1,)) elif form == 'CIFF': # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b # B2 = -(c_1 0... 0) ABCD[0, order + 1] = -c[0, 0] diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order, )) subdiag = diagonal_indices(ABCD[:order, :order], -1) ABCD[subdiag] = c[0, 1:] # C = (a_1 a_2... a_n) ABCD[order, :order] = a[0, :order] # supdiag = diagonal[odd + 1:order:2] - 1 supdiag = [i[odd:order+odd:2] for i in diagonal_indices(ABCD[:order, :order], +1)] ABCD[supdiag] = -g.reshape((-1,)) elif form == 'CRFBD': # C=(0 0...c_n) ABCD[order, order - 1] = c[0, order - 1] # B1 = (b_1 b_2... b_n), D=(b_n+1 0) ABCD[:, order] = b # B2 = -(a_1 a_2... a_n) ABCD[:order, order + 1] = -a diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order, )) # row numbers of delaying integrators dly = np.arange(odd, order, 2) subdiag = [np.atleast_1d(i)[dly] for i in diagonal_indices(ABCD[:order, :order], -1)] ABCD[subdiag] = c[0, dly] supdiag = [np.atleast_1d(i)[dly] for i in diagonal_indices(ABCD[:order, :order], +1)] supdiag = [np.atleast_1d(i[odd:]) for i in supdiag] ABCD[dly, :] = ABCD[dly, :] - np.dot(np.diag(g.reshape((-1))), \ ABCD[dly + 1, :]) if order > 2: coupl = np.arange(even + 1, order, 2) ABCD[coupl, :] = ABCD[coupl, :] + np.dot(np.diag(c[0, coupl - 1]), \ ABCD[coupl - 1, :]) elif form == 'CRFFD': diagonal = diagonal_indices(ABCD[:order, :order]) subdiag = diagonal_indices(ABCD[:order, :order], -1) supdiag = diagonal_indices(ABCD[:order, :order], +1) # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b.reshape((-1,)) # B2 = -(c_1 0... 0) ABCD[0, order + 1] = -c[0, 0] ABCD[diagonal] = np.ones((order, )) # rows to have c*(preceding row) added. multc = np.arange(1, order, 2) if order > 2: ABCD[[i[1::2] for i in subdiag]] = c[0, 2::2] if even: ABCD[[i[::2] for i in supdiag]] = -g.reshape((-1,)) else: # subtract g*(following row) from the multc rows ABCD[multc, :] = ABCD[multc, :] - np.dot(np.diag(g[0, :]), ABCD[multc + 1, :]) ABCD[multc, :] = ABCD[multc, :] + np.dot(np.diag(c[0, multc]), ABCD[multc - 1, :]) # C ABCD[order, 1:order:2] = a[0, 1:order:2] for i in range(0, order, 2): ABCD[order, :] = ABCD[order, :] + a[0, i]*ABCD[i, :] # The above gives y(n+1); need to add a delay to get y(n). # Do this by augmenting the states. Note: this means that # the apparent order of the NTF is one higher than it actually is. A, B, C, D = partitionABCD(ABCD, 2) A = np.vstack((np.hstack((A, np.zeros((order, 1)))), np.hstack((C, np.zeros((1, 1)))))) B = np.vstack((B, D)) C = np.hstack((np.zeros((1, order)), np.ones((1, 1)))) D = np.array([0, 0]).reshape(1, 2) ABCD = np.vstack((np.hstack((A, B)), np.hstack((C, D)))) elif 'PFF' == form: # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b.T odd_1 = odd # !! Bold assumption !! odd_2 = 0. # !! Bold assumption !! gc = g.dot(c[odd_1::2]) theta = np.arccos(1 - gc/2.) if odd_1: theta0 = 0. else: theta0 = theta[0] order_2 = 2*max(np.flatnonzero(np.abs(theta - theta0) > 0.5).shape) order_1 = order - order_2 # B2 = -(c_1 0...0 c_n 0...0) ABCD[0, order + 1] = -c[0] ABCD[order_1, order + 1] = -c[order_1] # number of elements = order diagonal = np.arange(0, order*(order + 1), order + 2) ABCD[diagonal] = np.ones((1, order)) i = np.hstack((np.arange(0, order_1 - 1, 2), np.arange(order - order_2, order, 2))) subdiag = diagonal[i] + 1 ABCD[subdiag] = c[i + 1] if odd_1: if order_1 >= 3: supdiag = diagonal[2:order_1:2] - 1 ABCD[supdiag] = -g[:(order_1 - 1)/2.] else: # rows to have g*(following row) subtracted. multg = np.arange(0, order_1, 2) ABCD[multg, :] = ABCD[multg, :] - np.diag(g[:order_1/2.]) * \ ABCD[multg + 1, :] if odd_2: if order_2 >= 3: supdiag = diagonal[order_1 + 1:order:2] - 1 ABCD[supdiag] = -g[:(order_1 - 1)/2.] else: # rows to have g*(following row) subtracted. multg = np.arange(order_1, order, 2) gg = g[(order_1 - odd_1)/2:] ABCD[multg, :] = ABCD[multg, :] - np.diag(gg)*ABCD[multg + 1, :] # Rows to have c*(preceding row) added. multc = np.hstack((np.arange(2, order_1, 2), np.arange(order_1 + 2, order, 2) )).reshape(1, -1) ABCD[multc, :] = ABCD[multc, :] + np.diag(c[multc])*ABCD[multc - 1, :] # C portion of ABCD i = np.hstack(( np.arange(0, order_1, 2), np.arange(order_1, order, 2) )).reshape(1,-1) ABCD[order, i] = a[i] for i in np.hstack((np.arange(1, order_1, 2), np.arange(order_1 + 1, order, 2) )): ABCD[order, :] = ABCD[order, :] + a[i] * ABCD[i, :] elif 'Stratos' == form: # code copied from case 'CIFF': # # B1 = (b_1 b_2... b_n), D=(b_(n+1) 0) ABCD[:, order] = b # # B2 = -(c_1 0... 0) ABCD[0, order + 1] = -c[0, 0] diagonal = diagonal_indices(ABCD[:order, :order]) ABCD[diagonal] = np.ones((order, )) subdiag = diagonal_indices(ABCD[:order, :order], -1) ABCD[subdiag] = c[0, 1:] # code based on case 'CRFF': # # rows to have g*(following row) subtracted. multg = np.arange(odd, order - 1, 2) ABCD[multg, :] = ABCD[multg, :] - np.dot(np.diag(g.reshape((-1,))), \ ABCD[multg + 1, :]) # code copied from case 'CIFF': # # C = (a_1 a_2... a_n) ABCD[order, :order] = a[0, :order] else: raise ValueError('Form %s is not yet supported.', form) return ABCD