Source code for deltasigma._evalTFP

# -*- coding: utf-8 -*-
# _evalTFP.py
# Module providing the evalTFP 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.

from __future__ import division

import numpy as np

from ._constants import eps
from ._evalRPoly import evalRPoly
from ._evalTF import evalTF
from ._utils import carray, restore_input_form, save_input_form


[docs]def evalTFP(Hs, Hz, f): """Evaluate a CT-DT transfer function product. Compute the value of a transfer function product Hs*Hz at a frequency f, where Hs is a continuous-time TF and Hz is a discrete-time TF. Use this function to evaluate the signal transfer function of a continuous-time (CT) system. In this context ``Hs`` is the open-loop response of the loop filter from the ``u`` input and ``Hz`` is the closed-loop noise transfer function. .. note:: This function attempts to **cancel poles** in ``Hs`` with **zeros** in ``Hz``. **Parameters:** Hs : tuple Hs is a CT (SISO) TF in zpk-tuple form. Hz : tuple Hz is a DT (SISO) TF in zpk-tuple form. f : scalar or 1D array or sequence Frequency values. **Returns:** H : scalar, or ndarray or sequence The calculated values: .. math:: H(f) = H_s(j2\\pi f)\ H_z(e^{j2\\pi f}) ``H`` has the same form as ``f``. .. seealso:: :func:`evalMixedTF`, a more advanced version of this function which is used to evaluate the individual feed-in transfer functions of a CT modulator. **Example:** Plot the STF of the 2nd-order CT system depicted at :func:`mapCtoD`.:: import numpy as np import pylab as plt from scipy.signal import lti from deltasigma import * from deltasigma._utils import _get_zpk Ac = np.array([[0, 0], [1, 0]]) Bc = np.array([[1, -1], [0, -1.5]]) Cc = np.array([[0, 1]]) Dc = np.array([[0, 0]]) LFc = (Ac, Bc, Cc, Dc) L0c = _get_zpk((Ac, Bc[:, 0].reshape((-1, 1)), Cc, Dc[0, 0].reshape(1, 1))) tdac = [0, 1] LF, Gp = mapCtoD(LFc, tdac) LF = lti(*LF) ABCD = np.vstack((np.hstack((LF.A, LF.B)), np.hstack((LF.C, LF.D)) )) H, _ = calculateTF(ABCD) # Yields H=(1-z^-1)^2 f = np.linspace(0, 2, 300); STF = evalTFP(L0c, _get_zpk(H), f) plt.figure(figsize=(12, 5)) plt.plot(f, dbv(STF)) plt.ylabel("|STF| [dB]") plt.xlabel("frequency ($1 \\\\rightarrow f_s$)") figureMagic((0, 2), .5, None, (-60, 0), 20, None) .. plot:: import numpy as np import pylab as plt from scipy.signal import lti from deltasigma import * from deltasigma._utils import _get_zpk Ac = np.array([[0, 0], [1, 0]]) Bc = np.array([[1, -1], [0, -1.5]]) Cc = np.array([[0, 1]]) Dc = np.array([[0, 0]]) LFc = (Ac, Bc, Cc, Dc) L0c = _get_zpk((Ac, Bc[:, 0].reshape((-1, 1)), Cc, Dc[0, 0].reshape(1, 1))) tdac = [0, 1] LF, Gp = mapCtoD(LFc, tdac) LF = lti(*LF) ABCD = np.vstack((np.hstack((LF.A, LF.B)), np.hstack((LF.C, LF.D)) )) H, _ = calculateTF(ABCD) # Yields H=(1-z^-1)^2 f = np.linspace(0, 2, 300); STF = evalTFP(L0c, _get_zpk(H), f) plt.figure(figsize=(12, 5)) plt.plot(f, dbv(STF)) plt.ylabel("|STF| [dB]") plt.xlabel("frequency ($1 \\\\rightarrow f_s$)") figureMagic((0, 2), .5, None, (-60, 0), 20, None) """ szeros, spoles, sk = Hs zzeros, zpoles, zk = Hz # sanitize f form_f = save_input_form(f) f = carray(f) if not np.prod(f.shape) == max(f.shape): raise ValueError("The f array must have shape " + "(N,) or (N, 1) or (N, 1, 1) ...") f = f.reshape((-1,)) # sanitize poles and zeros szeros, spoles = np.asarray(szeros), np.asarray(spoles) zzeros, zzeros = np.asarray(zzeros), np.asarray(zzeros) # back to business slim = min(0.001, max(1e-05, eps**(1./(1 + spoles.size)))) zlim = min(0.001, max(1e-05, eps**(1./(1 + zzeros.size)))) H = np.zeros(f.shape, dtype=np.complex128) w = 2*np.pi*f s = 1j*w z = np.exp(s) for i in range(f.shape[0]): wi = w[i] si = s[i] zi = z[i] if spoles.size == 0: cancel = False else: cancel = (np.abs(si - spoles) < slim) if not np.any(cancel): # wi is far from a pole, so just use the product Hs*Hz H[i] = evalTF(Hs, si)*evalTF(Hz, zi) else: # cancel pole(s) of Hs with corresponding zero(s) of Hz cancelz = np.abs(zi - zzeros) < zlim if np.sum(cancelz) > np.sum(cancel): H[i] = 0. else: if np.sum(cancelz) < np.sum(cancel): H[i] = np.Inf else: H[i] = evalRPoly(szeros, si, sk)*zi**np.sum(cancel) \ * evalRPoly(zzeros[~cancelz], zi, zk) \ / (evalRPoly(spoles[~cancel], si, 1.) \ * evalRPoly(zpoles, zi, 1.)) # return H matching the shape of f H = restore_input_form(H, form_f) return H