# -*- coding: utf-8 -*-
# _simulateQDSM.py
# Module providing the simulateQDSM function
# Copyright 2015 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 simulateQDSM() function
"""
from __future__ import division, print_function
import copy
import numpy as np
import scipy
from scipy.linalg import lstsq
from scipy.signal import freqz, tf2zpk
from ._config import _debug, setup_args
from ._ds_quantize import ds_quantize
from ._evalTF import evalTF
from ._partitionABCD import partitionABCD
from ._utils import carray, diagonal_indices, _is_zpk, _is_A_B_C_D, _is_num_den
# try to import and, if necessary, compile the Cython-optimized
# core of the simulation code`
try:
import pyximport
pyximport.install(setup_args=setup_args, inplace=True)
from ._simulateQDSM_core import simulateQDSM_core
except ImportError as e:
if _debug:
print(str(e))
# we'll just fall back to the Python version
pass
[docs]def simulateQDSM(u, arg2, nlev=2, x0=None):
"""Simulate a quadrature delta-sigma modulator.
This function computes the output of a quadrature delta-sigma modulator
corresponding to an input :math:`u`, with a description of the modulator, an
initial state :math:`x_0` (default all zeros) and a quantizer whose number
of levels is specified by :math:`n_{lev}`.
For multiple quantizers, make :math:`n_{lev}` a 1D vector, for complex
quantization to a diamond lattice, multiply :math:`n_{lev}` by :math:`j`.
Regarding the description of the modulator, it may be provided through an
ABCD matrix.
In this case, the shapes of the input parameters are:
* ``u.shape = (nu, N)``,
* ``nlev.shape = (nqi,)``,
* ``ABCD.shape = (order+nq, order+nq+nu)``.
Alternatively, the modulator may be described by a supported TF
representation, in particular it is recommended to use a zpk object. In this
case, the STF is assumed to be 1.
**Parameters:**
u : ndarray
The input signal to the modulator.
arg2 : ndarray or a supported LTI representation
A description of the modulator to simulate.
An ndarray instance is interpreted as an ABCD description. Equivalently,
the ABCD matrix may be supplied in ``(A, B, C, D)`` tuple form. All
other supported modulator specifications result in a conversion to a zpk
representation.
nlev : int or sequence-like, optional
The number of levels in the quantizer. If set to a sequence, each of the
elements is assumed to be the number of levels associated with a
quantizer. Defaults to ``2``.
x0 : float or sequence-like, optional
The initial states of the modulator. If it is set to a float, all states
are assumed to have the same value, ``x0``. If it is set to a
sequence-like object (list, tuple, 1D ndarray and similar), each entry is
assumed to be the value of one of the modulator states, in ascending
order. Defaults to ``0``.
**Returns:**
v : ndarray
The quantizer output.
xn : ndarray
The modulator states.
xmax : ndarray
The maximum value that each state reached during simulation.
y : ndarray
The quantizer input (ie the modulator output).
"""
if len(u.shape) == 1:
u = u.reshape((1, -1))
nu = u.shape[0]
if hasattr(nlev, '__len__'):
nlev = np.atleast_1d(nlev)
nq = max(nlev.shape)
else:
nq = 1
if isinstance(arg2, scipy.signal.lti):
k = arg2.k
zeros = np.asarray(arg2.z)
poles = np.asarray(arg2.p)
form = 2
order = max(zeros.shape)
elif _is_zpk(arg2):
zeros, poles, k = copy.deepcopy(arg2)
zeros = np.asarray(zeros)
poles = np.asarray(poles)
form = 2
order = max(zeros.shape)
elif isinstance(arg2, np.ndarray): # ABCD
if arg2.shape[1] > 2 and arg2.shape[1] == nu + arg2.shape[0]:
# ABCD dimesions OK
form = 1
ABCD = arg2
order = ABCD.shape[0] - nq
else:
raise ValueError('The ABCD argument does not have proper ' +
'dimensions.')
elif _is_A_B_C_D(arg2):
ABCD = np.vstack((np.hstack((np.atleast_2d(arg2[0]),
np.atleast_2d(arg2[1]))),
np.hstack((np.atleast_2d(arg2[2]),
np.atleast_2d(arg2[3])))))
form = 1
order = ABCD.shape[0] - nq
elif _is_num_den(arg2):
zeros, poles, k = tf2zpk(*arg2)
form = 2
order = max(zeros.shape)
else:
raise TypeError('The second argument is neither an ABCD matrix nor ' +
'an NTF.')
if x0 is None:
x0 = np.zeros(shape=(order, 1), dtype='complex128')
else:
x0 = carray(x0)
x0 = np.atleast_2d(x0).astype('complex128')
if form == 1:
A, B, C, D = partitionABCD(ABCD, nq + nu)
A = A.astype('complex128')
B = B.astype('complex128')
C = C.astype('complex128')
D = D.astype('complex128')
D1 = D[:, :nu].reshape((-1, nu))
else:
# Create a FF realization of 1-1/H.
# Note that MATLAB's zp2ss and canon functions don't work for complex
# TFs.
A = np.zeros(shape=(order, order), dtype='complex128')
B2 = np.vstack((np.atleast_2d(1), np.zeros(shape=(order-1, 1),
dtype='complex128')))
diag = diagonal_indices(A, 0)
A[diag] = zeros
subdiag = diagonal_indices(A, -1)
A[subdiag] = 1.
# Compute C st C*inv(zI-A)*B = 1-1/H(z);
w = 2*np.pi*np.random.rand(2*order)
desired = 1 - 1.0/evalTF((zeros, poles, k), np.exp(1j*w))
desired.reshape((1, -1))
# suppress warnings about complex TFs ???
sysresp = np.zeros((order, w.shape[0]), dtype='complex128')
for i in range(order):
Ctemp = np.zeros((1, order))
Ctemp[0, i] = 1
sys = (A, B2, Ctemp, np.zeros((1, 1)))
n, d = scipy.signal.ss2tf(*sys)
sysresp[i, :] = freqz(n[0, :], d, w)[1]
C = lstsq(sysresp.T, desired.T)[0].reshape((1, -1))
# !!!! Assume stf=1
B1 = -B2
B = np.hstack((B1, B2))
D1 = np.ones((1, 1), dtype='complex128')
v, xn, xmax, y = simulateQDSM_core(u, A, B, C, D1, order, nlev, nq, x0)
return v.squeeze(), xn.squeeze(), xmax, y.squeeze()