Source code for deltasigma._calculateTF

# -*- coding: utf-8 -*-
# Module providing the calculateTF 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
# LICENSE file for the licensing terms.

"""Module providing the calculateTF() function

from warnings import warn

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

from ._constants import eps
from ._partitionABCD import partitionABCD
from ._utils import carray, minreal

[docs]def calculateTF(ABCD, k=1.): """Calculate the NTF and STF of a delta-sigma modulator. The calculation is performed for a given loop filter ABCD matrix, assuming a quantizer gain of ``k``. **Parameters:** ABCD : array_like, The ABCD matrix that describes the system. k : float or ndarray-like, optional The quantizer gains. If only one quantizer is present, it may be set to a float, corresponding to the quantizer gain. If multiple quantizers are present, a list should be used, with quantizer gains ordered according to the order in which the quantizer inputs appear in the ``C`` and ``D`` submatrices. If not specified, a default of one quantizer with gain ``1.`` is assumed. **Returns:** (NTF, STF) : a tuple of two LTI objects (or of two lists of LTI objects). If a version of the ``scipy`` library equal to 0.16.x or greater is in use, the objects will be ``ZeroPolesGain`` objects, a subclass of ``scipy.signal.lti``. If the system has multiple quantizers, multiple STFs and NTFs will be returned. In that case: * ``STF[i]`` is the STF from ``u`` to output number ``i``. * ``NTF[i, j]`` is the NTF from the quantization noise of the quantizer number ``j`` to output number ``i``. **Note:** Setting ``k`` to a list is unsupported in the MATLAB code (last checked Nov. 2014). **Example:** Realize a fifth-order modulator with the cascade-of-resonators structure, feedback form. Calculate the ABCD matrix of the loop filter and verify that the NTF and STF are correct. .. code-block:: python from deltasigma import * H = synthesizeNTF(5, 32, 1) a, g, b, c = realizeNTF(H) ABCD = stuffABCD(a,g,b,c) ntf, stf = calculateTF(ABCD) From which we get: ``H``:: (z -1) (z^2 -1.997z +1) (z^2 -1.992z +0.9999) -------------------------------------------------------- (z -0.7778) (z^2 -1.796z +0.8549) (z^2 -1.613z +0.665) coefficients:: a: 0.0007, 0.0084, 0.055, 0.2443, 0.5579 g: 0.0028, 0.0079 b: 0.0007, 0.0084, 0.055, 0.2443, 0.5579, 1.0 c: 1.0, 1.0, 1.0, 1.0, 1.0 ABCD matrix:: [[ 1.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 6.75559806e-04 -6.75559806e-04] [ 1.00000000e+00 1.00000000e+00 -2.79396240e-03 0.00000000e+00 0.00000000e+00 8.37752565e-03 -8.37752565e-03] [ 1.00000000e+00 1.00000000e+00 9.97206038e-01 0.00000000e+00 0.00000000e+00 6.33294166e-02 -6.33294166e-02] [ 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 -7.90937431e-03 2.44344030e-01 -2.44344030e-01] [ 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 9.92090626e-01 8.02273699e-01 -8.02273699e-01] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 0.00000000e+00]] NTF:: (z -1) (z^2 -1.997z +1) (z^2 -1.992z +0.9999) -------------------------------------------------------- (z -0.7778) (z^2 -1.796z +0.8549) (z^2 -1.613z +0.665) STF:: 1 """ nq = len(k) if type(k) in (tuple, list) else 1 A, B, C, D = partitionABCD(ABCD, m=nq+1, r=nq) k = carray(k) diagk = np.atleast_2d(np.diag(k)) B1 = B[:, 0] B2 = B[:, 1:] B1 = B1.reshape((B1.shape[0], 1)) if len(B1.shape) == 1 else B1 B2 = B2.reshape((B2.shape[0], 1)) if len(B2.shape) == 1 else B2 # In a single-quantizer system, D2 should be all zeros, as any # non-zero term in D2 would imply we got a delay-free loop. # The original MATLAB code implies so. # The trouble arises when we adapt and extend the code to calculate # the transfer functions for multiple-quantizer systems. There the # off-diagonal terms of D2 may very well be non-zero, therefore # in the following we consider D2. # this means that if you supply an ABCD matrix, with one quantizer # and an (erroneously) zero-delay loop, the MATLAB toolbox will # disregard D2 and pretend it's zero and the loop is correctly # delayed, spitting out the corresponding TF. # Instead here we print out a warning and process the ABCD matrix # with the user-supplied, non-zero D2 matrix. # The resulting TFs will obviously be different. D1 = D[:, 0] D2 = D[:, 1:] D1 = D1.reshape((D1.shape[0], 1)) if len(D1.shape) == 1 else D1 D2 = D2.reshape((D2.shape[0], 1)) if len(D2.shape) == 1 else D2 # WARN DELAY FREE LOOPS if np.diag(D2).any(): warn("Delay free loop detected! D2 diag: %s", str(np.diag(D2))) # Find the noise transfer function by forming the closed-loop # system (sys_cl) in state-space form. Ct = np.linalg.inv(np.eye(nq) - D2*diagk) Acl = A +,,, C))) Bcl = np.hstack((B1 +,,, D1))),, Ct))) Ccl =,, C)) Dcl =, np.hstack((, D1), np.eye(nq)))) tol = min(1e-3, max(1e-6, eps**(1/ABCD.shape[0]))) ntfs = np.empty((nq, Dcl.shape[0]), dtype=np.object_) stfs = np.empty((Dcl.shape[0],), dtype=np.object_) # sweep the outputs 'cause scipy is silly but we love it anyhow. for i in range(Dcl.shape[0]): # input #0 is the signal # inputs #1,... are quantization noise stf_z, stf_p, stf_k = ss2zpk(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=0) stf = lti(stf_z, stf_p, stf_k) for j in range(nq): ntf_z, ntf_p, ntf_k = ss2zpk(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=j+1) ntf = lti(ntf_z, ntf_p, ntf_k) stf_min, ntf_min = minreal((stf, ntf), tol) ntfs[i, j] = ntf_min stfs[i] = stf_min # if we have one stf and one ntf, then just return those in a list if ntfs.shape == (1, 1): return [ntfs[0, 0], stfs[0]] return ntfs, stfs