Source code for deltasigma._pulse

# -*- coding: utf-8 -*-
# This module provides the pulse 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.

"""This module provides the pulse() function, which calculates the sampled
pulse response of a CT system.

from __future__ import division

import collections

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

from ._utils import _is_A_B_C_D, _is_num_den, _is_zpk, lcm, rat

[docs]def pulse(S, tp=(0., 1.), dt=1., tfinal=10., nosum=False): """Calculate the sampled pulse response of a CT system. ``tp`` may be an array of pulse timings, one for each input, or even a simple 2-elements tuple. **Parameters:** S : sequence A sequence of LTI objects specifying the system. The sequence S should be assembled so that ``S[i][j]`` returns the LTI system description from input ``i`` to the output ``j``. In the case of a MISO system, a unidimensional sequence ``S[i]`` is also acceptable. tp : array-like An (n, 2) array of pulse timings dt : scalar The time increment tfinal : scalar The time of the last desired sample nosum : bool A flag indicating that the responses are not to be summed **Returns:** y : ndarray The pulse response """ tp = np.asarray(tp) if len(tp.shape) == 1: if not tp.shape[0] == 2: raise ValueError("tp is not (n, 2)-shaped") tp = tp.reshape((1, tp.shape[0])) if len(tp.shape) == 2: if not tp.shape[1] == 2: raise ValueError("tp is not (n, 2)-shaped") # Compute the time increment dd = 1 for tpi in np.nditer(tp.T.copy(order='C')): _, di = rat(tpi, 1e-3) dd = lcm(di, dd) _, ddt = rat(dt, 1e-3) _, df = rat(tfinal, 1e-3) delta_t = 1./lcm(dd, lcm(ddt, df)) delta_t = max(1e-3, delta_t) # Put a lower limit on delta_t if (isinstance(S, collections.Iterable) and len(S)) \ and (isinstance(S[0], collections.Iterable) and len(S[0])) \ and (isinstance(S[0][0], lti) or _is_zpk(S[0][0]) or _is_num_den(S[0][0]) \ or _is_A_B_C_D(S[0][0])): pass else: S = list(zip(S)) #S[input][output] y1 = None for Si in S: y2 = None for So in Si: _, y2i = step2(So, T=np.arange(0., tfinal + delta_t, delta_t)) if y2 is None: y2 = y2i.reshape((y2i.shape[0], 1, 1)) else: y2 = np.concatenate((y2, y2i.reshape((y2i.shape[0], 1, 1))), axis=1) if y1 is None: y1 = y2 else: y1 = np.concatenate((y1, y2), axis=2) nd = int(np.round(dt/delta_t, 0)) nf = int(np.round(tfinal/delta_t, 0)) ndac = tp.shape[0] ni = len(S) # number of inputs if ni % ndac != 0: raise ValueError('The number of inputs must be divisible by the number of dac timings.') # Original comment from the MATLAB sources: # This requirement comes from the complex case, where the number of inputs # is 2 times the number of dac timings. I think this could be tidied up. # nis: Number of inputs grouped together with a common DAC timing # (2 for the complex case) nis = int(ni/ndac) # notice len(S[0]) is the number of outputs for us if not nosum: # Sum the responses due to each input set y = np.zeros((np.ceil(tfinal/float(dt)) + 1, len(S[0]), nis)) else: y = np.zeros((np.ceil(tfinal/float(dt)) + 1, len(S[0]), ni)) for i in range(ndac): n1 = int(np.round(tp[i, 0]/delta_t, 0)) n2 = int(np.round(tp[i, 1]/delta_t, 0)) z1 = (n1, y1.shape[1], nis) z2 = (n2, y1.shape[1], nis) yy = + np.concatenate((np.zeros(z1), y1[:nf-n1+1, :, i*nis:(i + 1)*nis]), axis=0) \ - np.concatenate((np.zeros(z2), y1[:nf-n2+1, :, i*nis:(i + 1)*nis]), axis=0) yy = yy[::nd, :, :] if not nosum: # Sum the responses due to each input set y = y + yy else: y[:, :, i] = yy.reshape(yy.shape[0:2]) return y