# -*- coding: utf-8 -*-
# _mapCtoD.py
# Module providing the mapCtoD 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 mapCtoD() function
"""
from __future__ import division, print_function
import collections
from warnings import warn
import numpy as np
from numpy.linalg import cond
from numpy.random import rand
from scipy.linalg import expm, inv, norm
from scipy.signal import cont2discrete, lti, ss2zpk
from ._constants import eps
from ._evalMixedTF import evalMixedTF
from ._padb import padb
from ._padr import padr
from ._utils import _getABCD
[docs]def mapCtoD(sys_c, t=(0, 1), f0=0.):
"""Map a MIMO continuous-time to an equiv. SIMO discrete-time system.
The criterion for equivalence is that the sampled pulse response
of the CT system must be identical to the impulse response of the DT system.
i.e. If ``yc`` is the output of the CT system with an input ``vc`` taken
from a set of DACs fed with a single DT input ``v``, then ``y``, the output
of the equivalent DT system with input ``v`` satisfies:
``y(n) = yc(n-)`` for integer ``n``. The DACs are characterized by
rectangular impulse responses with edge times specified in the t list.
**Input:**
sys_c : object
the LTI description of the CT system, which can be:
* the ABCD matrix,
* a list-like containing the A, B, C, D matrices,
* a list of zpk tuples (internally converted to SS representation),
* a list of LTI objects.
t : array_like
The edge times of the DAC pulse used to make CT waveforms
from DT inputs. Each row corresponds to one of the system
inputs; [-1 -1] denotes a CT input. The default is [0 1],
for all inputs except the first.
f0 : float
The (normalized) frequency at which the Gp filters' gains are
to be set to unity. Default 0 (DC).
**Output:**
sys : tuple
the LTI description for the DT equivalent, in A, B, C, D
representation.
Gp : list of lists
the mixed CT/DT prefilters which form the samples
fed to each state for the CT inputs.
**Example:**
Map the standard second order CT modulator shown below to its CT
equivalent and verify that its NTF is :math:`(1-z^{-1})^2`.
.. image:: ../doc/_static/mapCtoD.png
:align: center
:alt: mapCtoD block diagram
It can be done as follows::
from __future__ import print_function
import numpy as np
from scipy.signal import lti
from deltasigma import *
LFc = lti([[0, 0], [1, 0]], [[1, -1], [0, -1.5]], [[0, 1]], [[0, 0]])
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))
))
NTF, STF = calculateTF(ABCD)
print("NTF:") # after rounding to a 1e-6 resolution
print("Zeros:", np.real_if_close(np.round(NTF.zeros, 6)))
print("Poles:", np.real_if_close(np.round(NTF.poles, 6)))
Prints::
Zeros: [ 1. 1.]
Poles: [ 0. 0.]
Equivalent to::
(z -1)^2
NTF = ----------
z^2
.. seealso:: R. Schreier and B. Zhang, "Delta-sigma modulators employing \
continuous-time circuitry," IEEE Transactions on Circuits and Systems I, \
vol. 43, no. 4, pp. 324-332, April 1996.
"""
# You need to have A, B, C, D specification of the system
Ac, Bc, Cc, Dc = _getABCD(sys_c)
ni = Bc.shape[1]
# Sanitize t
if hasattr(t, 'tolist'):
t = t.tolist()
if (type(t) == tuple or type(t) == list) and np.isscalar(t[0]):
t = [t] # we got a simple list, like the default value
if not (type(t) == tuple or type(t) == list) and \
not (type(t[0]) == tuple or type(t[0]) == list):
raise ValueError("The t argument has an unrecognized shape")
# back to business
t = np.array(t)
if t.shape == (1, 2) and ni > 1:
t = np.vstack((np.array([[-1, -1]]), np.dot(np.ones((ni - 1, 1)), t)))
if t.shape != (ni, 2):
raise ValueError('The t argument has the wrong dimensions.')
di = np.ones(ni).astype(bool)
for i in range(ni):
if t[i, 0] == -1 and t[i, 1] == -1:
di[i] = False
# c2d assumes t1=0, t2=1.
# Also c2d often complains about poor scaling and can even produce
# incorrect results.
A, B, C, D, _ = cont2discrete((Ac, Bc, Cc, Dc), 1, method='zoh')
Bc1 = Bc[:, ~di]
# Examine the discrete-time inputs to see how big the
# augmented matrices need to be.
B1 = B[:, ~di]
D1 = D[:, ~di]
n = A.shape[0]
t2 = np.ceil(t[di, 1]).astype(np.int_)
esn = (t2 == t[di, 1]) and (D[0, di] != 0).T # extra states needed?
npp = n + np.max(t2 - 1 + 1*esn)
# Augment A to npp x npp, B to np x 1, C to 1 x np.
Ap = padb(padr(A, npp), npp)
for i in range(n + 1, npp):
Ap[i, i - 1] = 1
Bp = np.zeros((npp, 1))
if npp > n:
Bp[n, 0] = 1
Cp = padr(C, npp)
Dp = np.zeros((1, 1))
# Add in the contributions from each DAC
for i in np.flatnonzero(di):
t1 = t[i, 0]
t2 = t[i, 1]
B2 = B[:, i]
D2 = D[:, i]
if t1 == 0 and t2 == 1 and D2 == 0: # No fancy stuff necessary
Bp = Bp + padb(B2, npp)
else:
n1 = np.floor(t1)
n2 = np.ceil(t2) - n1 - 1
t1 = t1 - n1
t2 = t2 - n2 - n1
if t2 == 1 and D2 != 0:
n2 = n2 + 1
extraStateNeeded = 1
else:
extraStateNeeded = 0
nt = n + n1 + n2
if n2 > 0:
if t2 == 1:
Ap[:n, nt - n2:nt] = Ap[:n, nt - n2:nt] + np.tile(B2, (1, n2))
else:
Ap[:n, nt - n2:nt - 1] = Ap[:n, nt - n2:nt - 1] + np.tile(B2, (1, n2 - 1))
Ap[:n, (nt-1)] = Ap[:n, (nt-1)] + _B2formula(Ac, 0, t2, B2)
if n2 > 0: # pulse extends to the next period
Btmp = _B2formula(Ac, t1, 1, B2)
else: # pulse ends in this period
Btmp = _B2formula(Ac, t1, t2, B2)
if n1 > 0:
Ap[:n, n + n1 - 1] = Ap[:n, n + n1 - 1] + Btmp
else:
Bp = Bp + padb(Btmp, npp)
if n2 > 0:
Cp = Cp + padr(np.hstack((np.zeros((D2.shape[0], n + n1)), D2*np.ones((1, n2)))), npp)
sys = (Ap, Bp, Cp, Dp)
if np.any(~di):
# Compute the prefilters and add in the CT feed-ins.
# Gp = inv(sI - Ac)*(zI - A)/z*Bc1
n, m = Bc1.shape
Gp = np.empty_like(np.zeros((n, m)), dtype=object)
# !!Make this like stf: an array of zpk objects
ztf = np.empty_like(Bc1, dtype=object)
# Compute the z-domain portions of the filters
ABc1 = np.dot(A, Bc1)
for h in range(m):
for i in range(n):
if Bc1[i, h] == 0:
ztf[i, h] = (np.array([]),
np.array([0.]),
-ABc1[i, h]) # dt=1
else:
ztf[i, h] = (np.atleast_1d(ABc1[i, h]/Bc1[i, h]),
np.array([0.]),
Bc1[i, h]) # dt = 1
# Compute the s-domain portions of each of the filters
stf = np.empty_like(np.zeros((n, n)), dtype=object) # stf[out, in] = zpk
for oi in range(n):
for ii in range(n):
# Doesn't do pole-zero cancellation
stf[oi, ii] = ss2zpk(Ac, np.eye(n), np.eye(n)[oi, :],
np.zeros((1, n)), input=ii)
# scipy as of v 0.13 has no support for LTI MIMO systems
# only 'MISO', therefore you can't write:
# stf = ss2zpk(Ac, eye(n), eye(n), np.zeros(n, n)))
for h in range(m):
for i in range(n):
# k = 1 unneded, see below
for j in range(n):
# check the k values for a non-zero term
if stf[i, j][2] != 0 and ztf[j, h][2] != 0:
if Gp[i, h] is None:
Gp[i, h] = {}
Gp[i, h].update({'Hs':[list(stf[i, j])]})
Gp[i, h].update({'Hz':[list(ztf[j, h])]})
else:
Gp[i, h].update({'Hs':Gp[i, h]['Hs'] + [list(stf[i, j])]})
Gp[i, h].update({'Hz':Gp[i, h]['Hz'] + [list(ztf[j, h])]})
# the MATLAB-like cell code for the above statements would have
# been:
#Gp[i, h](k).Hs = stf[i, j]
#Gp[i, h](k).Hz = ztf[j, h]
#k = k + 1
if f0 != 0: # Need to correct the gain terms calculated by c2d
# B1 = gains of Gp @f0;
for h in range(m):
for i in range(n):
B1ih = np.real_if_close(evalMixedTF(Gp[i, h], f0))
# abs() used because ss() whines if B has complex entries...
# This is clearly incorrect.
# I've fudged the complex stuff by including a sign....
B1[i, h] = np.abs(B1ih) * np.sign(np.real(B1ih))
if np.abs(B1[i, h]) < 1e-09:
B1[i, h] = 1e-09 # This prevents NaN in "line 174" below
# Adjust the gains of the pre-filters
for h in range(m):
for i in range(n):
for j in range(max(len(Gp[i, h]['Hs']), len(Gp[i, h]['Hz']))):
# The next is "line 174"
Gp[i, h]['Hs'][j][2] = Gp[i, h]['Hs'][j][2]/B1[i, h]
sys = (sys[0], # Ap
np.hstack((padb(B1, npp), sys[1])), # new B
sys[2], # Cp
np.hstack((D1, sys[3]))) # new D
return sys, Gp
def _B2formula(Ac, t1, t2, B2):
if t1 == 0 and t2 == 0:
term = B2
return term
n = Ac.shape[0]
tmp = np.eye(n) - expm(-Ac)
if cond(tmp) < 1000000.0:
term = np.dot(((expm(-Ac*t1) - expm(-Ac*t2))*inv(tmp)), B2)
return term
# Numerical trouble. Perturb slightly and check the result
ntry = 0
k = np.sqrt(eps)
Ac0 = Ac
while ntry < 2:
Ac = Ac0 + k*rand(n,n)
tmp = np.eye(n) - expm(-Ac)
if cond(tmp) < 1/np.sqrt(eps):
ntry = ntry + 1
if ntry == 1:
term = np.dot(np.dot(expm(-Ac*t1) - expm(-Ac*t2), inv(tmp)), B2)
else:
term1 = np.dot(np.dot(expm(-Ac*t1) - expm(-Ac*t2), inv(tmp)), B2)
k = k*np.sqrt(2)
if norm(term1 - term) > 0.001:
warn('Inaccurate calculation in mapCtoD.')
return term