Merge branch 'v02' of https://github.com/ggventurini/python-deltasigma

deltasigma/_DocumentNTF.py

@@ -102,7 +102,7 @@ def DocumentNTF(arg1, osr=64, f0=0, quadrature=False):
         plt.plot(-np.array([f1, f2]), ING0*np.array([1, 1]), 'k', linewidth=3)
         plt.text(-f0, ING0 + 1, '%.0fdB' % ING0, horizontalalignment = 'center',
                  verticalalignment = 'bottom')
-        f_left = - 0.5
+        f_left = -0.5
     else:
         f_left = 0
     # variable 'f_left' used before assignment in DocumentNTF.m #REP
@@ -112,7 +112,8 @@ def DocumentNTF(arg1, osr=64, f0=0, quadrature=False):
         plt.text(f_left, 0, msg, horizontalalignment = 'left', verticalalignment = 'top')
     plt.grid(True)
     y_bot = min(-80, np.round(NG0*1.1, -1))
-    figureMagic(xRange=(f_left, 0.5), dx=1/20., yRange=(y_bot, 15), dy=10)
+    dx = 1./10 if quadrature else 1./20
+    figureMagic(xRange=(f_left, 0.5), dx=dx, yRange=(y_bot, 15), dy=10)
     plt.ylabel('|H(f)| dB')
     plt.xlabel('Frequency ($1 \\rightarrow f_{s}$)')
     plt.title('Frequency Response')

deltasigma/_PlotExampleSpectrum.py

@@ -32,6 +32,7 @@
 from ._evalTF import evalTF
 from ._figureMagic import figureMagic
 from ._simulateDSM import simulateDSM
+from ._simulateQDSM import simulateQDSM
 from ._undbv import undbv
 
 
@@ -63,13 +64,6 @@ def PlotExampleSpectrum(ntf, M=1, osr=64, f0=0, quadrature=False):
                  Whether the delta sigma modulator is a quadrature
                  modulator or not. Defaults to ``False``.
 
-    .. note::
-
-        Quadrature modulators require the following currently unimplemented
-        functions: :func:`simulateQDSM` (expected v. 0.2).
-        Setting ``quadrature`` to ``True`` results in a
-        ``NotImplementedError`` being raised.
-
     .. plot::
 
         import pylab as plt
@@ -98,11 +92,9 @@ def PlotExampleSpectrum(ntf, M=1, osr=64, f0=0, quadrature=False):
         u = Amp*M*np.cos((2*np.pi/N)*fin*t)
         v, _, xmax, y = simulateDSM(u, ntf, M+1)
     else:
-        raise NotImplementedError("The required simulateQDSM function " + \
-                                  "is not available yet.")
         t = np.arange(0, N).reshape((1, -1))
         u = Amp*M*np.exp((2j*np.pi/N)*fin*t)
-        v, xn, xmax, y = simulateQDSM(u, ntf, M + 1)
+        v, _, xmax, y = simulateQDSM(u, ntf, M+1)
     window = ds_hann(N)
     NBW = 1.5/N
     spec0 = fft(v * window)/(M*N/4)
@@ -131,13 +123,13 @@ def PlotExampleSpectrum(ntf, M=1, osr=64, f0=0, quadrature=False):
         freq = np.linspace(-0.5, 0.5, N + 1)
         freq = freq[:-1]
         plt.plot(freq, dbv(spec0), 'c', linewidth=1)
-        plt.hold('on')
-        spec_smoothed = circ_smooth(abs(spec0) ** 2, 16)
+        plt.hold(True)
+        spec_smoothed = circ_smooth(abs(spec0)**2, 16)
         plt.plot(freq, dbp(spec_smoothed), 'b', linewidth=3)
         Snn = abs(evalTF(ntf, np.exp(2j * np.pi * freq))) ** 2 * 2 / 12 * (delta / M) ** 2
-        plt.plot(freq, dbp(Snn * NBW), 'm', linewidth=1)
-        snr = calculateSNR(spec0[N/2 + np.arange(f1_bin, f2_bin + 1)], fin - f1_bin)
-        msg = 'SQNR  =  %.1fdB\\n @ A = %.1fdBFS & osr = %.0f\\n' % \
+        plt.plot(freq, dbp(Snn*NBW), 'm', linewidth=1)
+        snr = calculateSNR(spec0[N/2 + f1_bin:N/2 + f2_bin + 1], fin - f1_bin)
+        msg = 'SQNR  =  %.1fdB\n @ A = %.1fdBFS & osr = %.0f' % \
               (snr, dbv(spec0[N/2 + fin]), osr)
         if f0 >=  0:
             plt.text(f0 - 0.05, - 15, msg, horizontalalignment='right',

deltasigma/_calculateQTF.py

@@ -0,0 +1,102 @@
+# -*- coding: utf-8 -*-
+# _calculateQTF.py
+# Module providing the calculateQTF 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 calculateQTF() function
+"""
+
+from __future__ import division
+
+import numpy as np
+
+from scipy.signal import ss2tf, tf2zpk
+
+from ._cancelPZ import cancelPZ
+from ._partitionABCD import partitionABCD
+
+
+def calculateQTF(ABCDr):
+    """Calculate noise and signal transfer functions for a quadrature modulator
+
+    **Parameters:**
+
+    ABCDr : ndarray
+        The ABCD matrix, in real form. You may call :func:`mapQtoR` to convert
+        an imaginary (quadrature) ABCD matrix to a real one.
+
+    **Returns:**
+
+    ntf, stf, intf, istf : tuple of zpk tuples
+        The quadrature noise and signal transfer functions.
+
+    :raises RuntimeError: if the supplied ABCD matrix results in denominator mismatches.
+    """
+    A, B, C, D = partitionABCD(ABCDr, 4)
+
+    #Construct an ABCD description of the closed-loop system
+    # sys is a tuple in A, B, C, D form
+    Acl = A + np.dot(B[:, 2:4], C)
+    Bcl = B
+    Ccl = C
+    Dcl = np.hstack((D[:, 0:2], np.eye(2)))
+    #sys = (A + np.dot(B[:, 2:4], C), B, C,
+    #       np.hstack((D[:, 0:2], np.eye(2))))
+
+    #Calculate the 2x4 matrix of transfer functions
+    tfs = np.empty((2, 4), dtype=object)
+    # Each tf is a tuple in num, den form
+    # tf[i, j] corresponds to the TF from input j to output i
+    for i in range(2):
+        for j in range(4):
+            tfs[i, j] = ss2tf(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=j)
+
+    #Reduce these to NTF, STF, INTF and ISTF
+    if any(tfs[0, 2][1] != tfs[1, 3][1]):
+        raise RuntimeError('TF Denominator mismatch. Location 1')
+
+    ntf_x = (0.5 * (tfs[0, 2][0] + tfs[1, 3][0]), tfs[0, 2][1])
+    intf_x = (0.5 * (tfs[0, 2][0] - tfs[1, 3][0]), tfs[0, 2][1])
+
+    if any(tfs[0, 3][1] != tfs[1, 2][1]):
+        raise RuntimeError('TF Denominator mismatch. Location 2')
+
+    ntf_y = (0.5 * (tfs[1, 2][0] - tfs[0, 3][0]), tfs[1, 2][1])
+    intf_y = (0.5 * (tfs[1, 2][0] + tfs[0, 3][0]), tfs[1, 2][1])
+
+    if any(ntf_x[1] != ntf_y[1]):
+        raise RuntimeError('TF Denominator mismatch. Location 3')
+    if any(tfs[0, 0][1] != tfs[1, 1][1]):
+        raise RuntimeError('TF Denominator mismatch. Location 4')
+
+    stf_x = (0.5 * (tfs[0, 0][0] + tfs[1, 1][0]), tfs[0, 0][1])
+    istf_x = (0.5 * (tfs[0, 0][0] - tfs[1, 1][0]), tfs[0, 0][1])
+
+    if any(tfs[0, 1][1] != tfs[1, 0][1]):
+        raise RuntimeError('TF Denominator mismatch. Location 5')
+
+    stf_y = (0.5 * (tfs[1, 0][0] - tfs[0, 1][0]), tfs[1, 0][1])
+    istf_y = (0.5 * (tfs[1, 0][0] + tfs[0, 1][0]), tfs[1, 0][1])
+
+    if any(stf_x[1] != stf_y[1]):
+        raise RuntimeError('TF Denominator mismatch. Location 6')
+
+    # suppress warnings about complex TFs
+    #warning('off')
+    ntf = cancelPZ(tf2zpk(ntf_x[0] + 1j*ntf_y[0], ntf_x[1]))
+    intf = cancelPZ(tf2zpk(intf_x[0] + 1j*intf_y[0], intf_x[1]))
+    stf = cancelPZ(tf2zpk(stf_x[0] + 1j*stf_y[0], ntf_x[1]))
+    istf = cancelPZ(tf2zpk(istf_x[0] + 1j*istf_y[0], intf_x[1]))
+    #warning('on')
+
+    return ntf, stf, intf, istf

deltasigma/_calculateTF.py

@@ -169,8 +169,8 @@ def calculateTF(ABCD, k=1.):
     for i in range(Dcl.shape[0]):
         # input #0 is the signal
         # inputs #1,... are quantization noise
-        stf_p, stf_z, stf_k  = ss2zpk(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=0)
-        stf = lti(stf_p, stf_z, stf_k)
+        stf_z, stf_p, stf_k  = ss2zpk(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=0)
+        stf = lti(stf_z, stf_p, stf_k)
         for j in range(nq):
             ntf_z, ntf_p, ntf_k = ss2zpk(Acl, Bcl, Ccl[i, :], Dcl[i, :], input=j+1)
             ntf = lti(ntf_z, ntf_p, ntf_k)

deltasigma/_partitionABCD.py

@@ -91,6 +91,10 @@ def partitionABCD(ABCD, m=None, r=None):
     # remember the ABCD matrix is assembled like this:
     # [[A, B],
     #  [C, D]]
+    #
+    # n is the number of states,
+    # m is the number of inputs
+    # r is the number of outputs
     if m is None:
         # we guess: the system has either one input or one output. Most
         # likely one output, but we keep this for consistency with the

deltasigma/_realizeNTF.py

@@ -363,7 +363,7 @@ def realizeNTF(ntf, form='CRFB', stf=None):
         # 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 = np.zeros((1, order + 1))
             bi[i] = 1
             ABCD = stuffABCD(a, g, bi, c, form)
             if form[2:4] == 'FF':