Source code for deltasigma._realizeNTF

# -*- coding: utf-8 -*-
# _realizeNTF.py
# Module providing realizeNTF
# 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 realizeNTF()
"""

from __future__ import division, print_function

from warnings import warn

import numpy as np

from ._calculateTF import calculateTF
from ._evalRPoly import evalRPoly
from ._evalTF import evalTF
from ._stuffABCD import stuffABCD
from ._utils import _get_zpk, carray, cplxpair


[docs]def realizeNTF(ntf, form='CRFB', stf=None): """Convert an NTF into coefficients for the desired structure. **Parameters:** ntf : a zpk transfer function The modulator noise transfer function (NTF) form : string A structure identifier. Supported structures are: * *"CRFB"*: Cascade of resonators, feedback form. * *"CRFF"*: Cascade of resonators, feedforward form. * *"CIFB"*: Cascade of integrators, feedback form. * *"CIFF"*: Cascade of integrators, feedforward form. * *"CRFBD"*: CRFB with delaying quantizer. * *"CRFFD"*: CRFF with delaying quantizer. * *"PFF"*: Parallel feed-forward. * *"Stratos"*: A CIFF-like structure with non-delaying resonator feedbacks, contributed to the MATLAB Delta Sigma toolbox in 2007 by Jeff Gealow. See the accompanying documentation (:ref:`topologies-diagrams`) for block diagrams of each structure. .. note: The order of the NTF zeros must be (real, complex conj. pairs). The order of the zeros is used when mapping the NTF onto the chosen topology. stf : a zpk transfer function, optional the Signal Transfer Function **Returns:** a, g, b, c : tuple of ndarrays the coefficients for the desired structure **Example:** Determine the coefficients for a 5th-order modulator with the cascade-of-resonators structure, feedback (CRFB) form.:: from deltasigma import synthesizeNTF, realizeNTF H = synthesizeNTF(5, 32, 1) a, g, b, c = realizeNTF(H,'CRFB') Returns the values:: 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 """ # The basic idea is to equate the loop filter at a set of # points in the z-plane to L1 = 1-1/ntf at those points. # Code common to all functions\ if (hasattr(ntf, 'inputs') and not ntf.inputs == 1) or \ (hasattr(ntf, 'outputs') and not ntf.outputs == 1): raise TypeError("Only SISO transfer functions can be evaluated.") ntf_z, ntf_p, _ = _get_zpk(ntf) ntf_z = carray(ntf_z) ntf_p = carray(ntf_p) order = max(ntf_p.shape) order2 = int(np.floor(order / 2)) odd = (order - 2 * order2 == 1) a = np.zeros((1, order)) g = np.zeros((1, order2)) b = np.zeros((1, order + 1)) c = np.ones((1, order)) T = np.zeros((order, order * 2), dtype='complex') # instead of asuming enforce roots order ntf_z = cplxpair(ntf_z)[::-1] N = 200 min_distance = 0.09 zSet = np.zeros((N, ), dtype='complex') j = 0 for i in range(1, N + 1): z = 1.1 * np.exp(2j * np.pi * i / N) if np.all(np.abs(ntf_z - z) > min_distance): zSet[j] = z j = j + 1 zSet = zSet[:j] if form == 'CRFB': # Find g # Assume the roots are ordered, real first, then cx conj. pairs for i in range(0, order2): g[0, i] = 2 * (1 - np.real(ntf_z[2 * i + 1 * odd])) L1 = np.zeros((1, order * 2), dtype='complex') # Form the linear matrix equation a*T* = L1 for i in range(0, 2 * order): z = zSet[i] L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) Dfactor = (z - 1) / z product = 1 for j in range(order, 0 + 1 * odd, -2): product = z / evalRPoly(ntf_z[j - 2:j], z) * product T[j - 1, i] = product * Dfactor T[j - 2, i] = product if odd: T[0, i] = product / (z - 1) a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b[0,:order] = a[0, :] b[0, order] = 1 elif form == 'CRFF': # Find g # Assume the roots are ordered, real first, then cx conj. pairs for i in range(0, order2): g[0, i] = 2 * (1 - np.real(ntf_z[2 * i + 1 * odd])) L1 = np.zeros((1, order * 2), dtype='complex') # Form the linear matrix equation a*T*=L1 for i in range(0, 2 * order): z = zSet[i] # L1(z) = 1-1/H(z) L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) if odd: Dfactor = z - 1 product = 1 / Dfactor T[0, i] = product else: Dfactor = (z - 1) / z product = 1. for j in range(0 + 1 * odd, order, 2): product = z / evalRPoly(ntf_z[j:j + 2], z) * product T[j, i] = product * Dfactor T[j + 1, i] = product a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b = np.hstack((np.atleast_2d(1), np.zeros((1, order - 1)), np.atleast_2d(1) )) elif form == 'CIFB': if any(abs(np.real(ntf_z) - 1) > 0.001): warn("The ntf's zeros have had their real parts set to one.") ntf_z = 1 + 1j * np.imag(ntf_z) for i in range(order2): g[0, i] = np.imag(ntf_z[2 * i + 1 * odd]) ** 2 L1 = np.zeros((1, order * 2), dtype='complex') for i in range(order * 2): z = zSet[i] L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) Dfactor = z - 1. product = 1. for j in range(order - 1, 1 * odd, -2): product = product / evalRPoly(ntf_z[j - 1:j + 1], z) T[j, i] = product * Dfactor T[j - 1, i] = product if odd: T[0, i] = product / (z - 1.) a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b[0,:order] = a b[0, order] = 1. elif form == 'CIFF': if (abs(np.real(ntf_z) - 1) > 0.001).any(): warn("The ntf's zeros have had their real parts set to one.") ntf_z = 1. + 1j * np.imag(ntf_z) for i in range(order2): g[0, i] = np.imag(ntf_z[2 * i + 1 * odd]) ** 2 L1 = np.zeros((1, order * 2), dtype='complex') for i in range(order * 2): z = zSet[i] L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) Dfactor = (z - 1) if odd: product = 1. / (z - 1.) T[0, i] = product else: product = 1 for j in range(0 + 1 * odd, order - 1, 2): product = product / evalRPoly(ntf_z[j:j + 2], z) T[j, i] = product * Dfactor T[j + 1, i] = product a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b = np.zeros((1, order + 1)) b[0, 0] = 1. b[0, -1] = 1. elif form == 'CRFBD': for i in range(order2): g[0, i] = 2 * (1. - np.real(ntf_z[2 * i + 1 * odd])) L1 = np.zeros((1, order * 2), dtype='complex') for i in range(order * 2): z = zSet[i] L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) Dfactor = (z - 1.) product = 1. / z for j in range(order - 1, 0 + 1 * odd, -2): product = z / evalRPoly(ntf_z[j - 1:j + 1], z) * product T[j, i] = product * Dfactor T[j - 1, i] = product if odd: T[0, i] = product * z / (z - 1.) a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b = np.hstack((a, np.zeros((1, 1)))) b[0, order] = 1. elif form == 'CRFFD': for i in range(order2): g[0, i] = 2 * (1. - np.real(ntf_z[2 * i + 1 * odd])) zL1 = np.zeros((1, order*2), dtype='complex') # Form the linear matrix equation a*T*=zL1 for i in range(order * 2): z = zSet[i] # zL1 = z*(1-1/H(z)) zL1[0, i] = z*(1. - 1.0/evalTF(ntf, z)) if odd: Dfactor = (z - 1.) / z product = 1. / Dfactor T[0, i] = product else: Dfactor = z - 1 product = 1 for j in range(0 + 1 * odd, order, 2): product = z / evalRPoly(ntf_z[j:j + 2], z) * product T[j, i] = product * Dfactor T[j + 1, i] = product a = -np.real(np.linalg.lstsq(T.T, zL1.T)[0]).T if stf is None: b = np.hstack((np.atleast_2d(1), np.zeros((1, order - 1)), np.atleast_2d(1)) ) elif form == 'PFF': warn("Untested code accessed: realizeNTF('PFF')") # Find g # Assume the roots are ordered, real first, then cx conj. pairs # with the secondary zeros after the primary zeros for i in range(order2): g[0, i] = 2 * (1. - np.real(ntf_z[2 * i + 1 * odd])) # Find the dividing line between the zeros theta0 = np.abs(np.angle(ntf_z[0])) # !! 0.5 radians is an arbitrary separation !! i = np.flatnonzero(np.abs(np.abs(np.angle(ntf_z)) - theta0) > 0.5) if i.shape[0]: order_1 = i[0] - 1 else: order_1 = 0 order_2 = order - order_1 if i.shape[0] != order_2: raise ValueError('For the PFF form, the NTF zeros must ' + 'be sorted into primary and secondary zeros') odd_1 = order_1 % 2 odd_2 = order_2 % 2 L1 = np.zeros((1, order * 2), dtype='complex') # Form the linear matrix equation a*T*=L1 for i in range(order * 2): z = zSet[i] # L1(z) = 1-1/H(z) L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) if odd_1: Dfactor = z - 1 product = 1. / Dfactor T[0, i] = product else: Dfactor = (z - 1) / z product = 1 for j in range(odd_1, order_1, 2): product = z / evalRPoly(ntf_z[j:j + 2], z) * product T[j, i] = product * Dfactor T[j + 1, i] = product if odd_2: Dfactor = z - 1 product = 1. / Dfactor T[order_1, i] = product else: Dfactor = (z - 1) / z product = 1. for j in range(order_1 + odd_2, order, 2): product = z / evalRPoly(ntf_z[j:j + 2], z) * product T[j, i] = product * Dfactor T[j + 1, i] = product a = -np.real(np.linalg.lstsq(T.T, zL1.T)[0]).T if stf is None: b = np.hstack((np.atleast_2d(1), np.zeros((1, order_1 - 1)), np.atleast_2d(1), np.zeros(1, order_2 - 1), np.atleast_2d(1) )) elif form == 'Stratos': # code copied from case 'CRFF': # Find g # Assume the roots are ordered, real first, then cx conj. pairs for i in range(order2): g[0, i] = 2 * (1 - np.real(ntf_z[2 * i + 1 * odd])) # code copied from case 'CIFF': L1 = np.zeros((1, order * 2), dtype='complex') # Form the linear matrix equation a*T*=L1 for i in range(order * 2): z = zSet[i] # L1(z) = 1-1/H(z) L1[0, i] = 1. - evalRPoly(ntf_p, z) / evalRPoly(ntf_z, z) Dfactor = (z - 1.) if odd: product = 1. / (z - 1.) T[0, i] = product else: product = 1. for j in range(0 + 1 * odd, order - 1, 2): product = product / evalRPoly(ntf_z[j:j + 2], z) T[j, i] = product * Dfactor T[j + 1, i] = product a = -np.real(np.linalg.lstsq(T.T, L1.T)[0]).T if stf is None: b = np.hstack(( np.atleast_2d(1), np.zeros((1, order - 1)), np.atleast_2d(1) )) if stf is not None: # Compute the TF from each feed-in to the output # and solve for coefficients which yield the best match # THIS CODE IS NOT OPTIMAL, in terms of computational efficiency. stfList = [] for i in range(0, order + 1): bi = np.zeros((1, order + 1)) bi[i] = 1 ABCD = stuffABCD(a, g, bi, c, form) if form[2:4] == 'FF': ABCD[0, order + 1] = -1 _, stfListi = calculateTF(ABCD) stfList.append(stfListi) # Build the matrix equation b A = x and solve it. A = np.zeros((order + 1, max(zSet.shape))) for i in range(0, order + 1): A[i, :] = evalTF(stfList[i], zSet) x = evalTF(stf, zSet) x = x[:].T b = np.real(x / A) return a.squeeze(), g.squeeze(), b.squeeze(), c.squeeze()