[Scipy-svn] r6241 - in trunk/scipy/signal: . tests

scipy-svn@scip... scipy-svn@scip...
Fri Feb 19 05:01:32 CST 2010


Author: stefan
Date: 2010-02-19 05:01:31 -0600 (Fri, 19 Feb 2010)
New Revision: 6241

Modified:
   trunk/scipy/signal/tests/test_waveforms.py
   trunk/scipy/signal/waveforms.py
Log:
Rewrite of chirp by Warren Weckesser.

Modified: trunk/scipy/signal/tests/test_waveforms.py
===================================================================
--- trunk/scipy/signal/tests/test_waveforms.py	2010-02-17 03:25:50 UTC (rev 6240)
+++ trunk/scipy/signal/tests/test_waveforms.py	2010-02-19 11:01:31 UTC (rev 6241)
@@ -1,14 +1,246 @@
-#this program corresponds to special.py
 
 import numpy as np
-from numpy.testing import *
+from numpy.testing import TestCase, assert_almost_equal, assert_equal, assert_, \
+        assert_raises, run_module_suite
 
-import scipy.signal as signal
+import scipy.signal.waveforms as waveforms
 
+
+# These chirp_* functions are the instantaneous frequencies of the signals
+# returned by chirp().
+
+def chirp_linear(t, f0, f1, t1):
+    f = f0 + (f1 - f0) * t / t1
+    return f
+
+def chirp_quadratic(t, f0, f1, t1, vertex_zero=True):
+    if vertex_zero:
+        f = f0 + (f1 - f0) * t**2 / t1**2
+    else:
+        f = f1 - (f1 - f0) * (t1 - t)**2 / t1**2
+    return f
+
+def chirp_geometric(t, f0, f1, t1):
+    f = f0 * (f1/f0)**(t/t1)
+    return f
+
+def chirp_hyperbolic(t, f0, f1, t1):
+    f = f0*f1*t1 / ((f0 - f1)*t + f1*t1)
+    return f
+
+
+def compute_frequency(t, theta):
+    """Compute theta'(t)/(2*pi), where theta'(t) is the derivative of theta(t)."""
+    # Assume theta and t are 1D numpy arrays.
+    # Assume that t is uniformly spaced.
+    dt = t[1] - t[0]
+    f = np.diff(theta)/(2*np.pi) / dt
+    tf = 0.5*(t[1:] + t[:-1])
+    return tf, f
+
+
 class TestChirp(TestCase):
-    def test_log_chirp_at_zero(self):
-        assert_almost_equal(signal.waveforms.chirp(t=0, method='log'),
-                            1.0)
 
+    def test_linear_at_zero(self):
+        w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='linear')
+        assert_almost_equal(w, 1.0)
+
+    def test_linear_freq_01(self):
+        method = 'linear'
+        f0 = 1.0
+        f1 = 2.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 100)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_linear(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_linear_freq_02(self):
+        method = 'linear'
+        f0 = 200.0
+        f1 = 100.0
+        t1 = 10.0
+        t = np.linspace(0, t1, 100)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_linear(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_quadratic_at_zero(self):
+        w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='quadratic')
+        assert_almost_equal(w, 1.0)
+
+    def test_quadratic_at_zero2(self):
+        w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='quadratic',
+                                                                vertex_zero=False)
+        assert_almost_equal(w, 1.0)
+
+    def test_quadratic_freq_01(self):
+        method = 'quadratic'
+        f0 = 1.0
+        f1 = 2.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 2000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_quadratic(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_quadratic_freq_02(self):
+        method = 'quadratic'
+        f0 = 20.0
+        f1 = 10.0
+        t1 = 10.0
+        t = np.linspace(0, t1, 2000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_quadratic(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_logarithmic_at_zero(self):
+        w = waveforms.chirp(t=0, f0=1.0, f1=2.0, t1=1.0, method='logarithmic')
+        assert_almost_equal(w, 1.0)
+
+    def test_logarithmic_freq_01(self):
+        method = 'logarithmic'
+        f0 = 1.0
+        f1 = 2.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 10000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_geometric(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_logarithmic_freq_02(self):
+        method = 'logarithmic'
+        f0 = 200.0
+        f1 = 100.0
+        t1 = 10.0
+        t = np.linspace(0, t1, 10000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_geometric(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_logarithmic_freq_03(self):
+        method = 'logarithmic'
+        f0 = 100.0
+        f1 = 100.0
+        t1 = 10.0
+        t = np.linspace(0, t1, 10000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_geometric(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_hyperbolic_at_zero(self):
+        w = waveforms.chirp(t=0, f0=10.0, f1=1.0, t1=1.0, method='hyperbolic')
+        assert_almost_equal(w, 1.0)
+
+    def test_hyperbolic_freq_01(self):
+        method = 'hyperbolic'
+        f0 = 10.0
+        f1 = 1.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 10000)
+        phase = waveforms._chirp_phase(t, f0, t1, f1, method)
+        tf, f = compute_frequency(t, phase)
+        abserr = np.max(np.abs(f - chirp_hyperbolic(tf, f0, f1, t1)))
+        assert_(abserr < 1e-6)
+
+    def test_hyperbolic_freq_02(self):
+        method = 'hyperbolic'
+        f0 = 10.0
+        f1 = 100.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 10)
+        assert_raises(ValueError, waveforms.chirp, t, f0, t1, f1, method)
+
+    def test_hyperbolic_freq_03(self):
+        method = 'hyperbolic'
+        f0 = -10.0
+        f1 = 0.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 10)
+        assert_raises(ValueError, waveforms.chirp, t, f0, t1, f1, method)
+
+    def test_unknown_method(self):
+        method = "foo"
+        f0 = 10.0
+        f1 = 20.0
+        t1 = 1.0
+        t = np.linspace(0, t1, 10)        
+        assert_raises(ValueError, waveforms.chirp, t, f0, t1, f1, method)
+
+
+class TestSweepPoly(TestCase):
+
+    def test_sweep_poly_quad1(self):
+        p = np.poly1d([1.0, 0.0, 1.0])
+        t = np.linspace(0, 3.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = p(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_const(self):
+        p = np.poly1d(2.0)
+        t = np.linspace(0, 3.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = p(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_linear(self):
+        p = np.poly1d([-1.0, 10.0])
+        t = np.linspace(0, 3.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = p(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_quad2(self):
+        p = np.poly1d([1.0, 0.0, -2.0])
+        t = np.linspace(0, 3.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = p(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_cubic(self):
+        p = np.poly1d([2.0, 1.0, 0.0, -2.0])
+        t = np.linspace(0, 2.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = p(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_cubic2(self):
+        """Use an array of coefficients instead of a poly1d."""
+        p = np.array([2.0, 1.0, 0.0, -2.0])
+        t = np.linspace(0, 2.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = np.poly1d(p)(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
+    def test_sweep_poly_cubic3(self):
+        """Use a list of coefficients instead of a poly1d."""
+        p = [2.0, 1.0, 0.0, -2.0]
+        t = np.linspace(0, 2.0, 10000)
+        phase = waveforms._sweep_poly_phase(t, p)
+        tf, f = compute_frequency(t, phase)
+        expected = np.poly1d(p)(tf)
+        abserr = np.max(np.abs(f - expected))
+        assert_(abserr < 1e-6)
+
 if __name__ == "__main__":
     run_module_suite()

Modified: trunk/scipy/signal/waveforms.py
===================================================================
--- trunk/scipy/signal/waveforms.py	2010-02-17 03:25:50 UTC (rev 6240)
+++ trunk/scipy/signal/waveforms.py	2010-02-19 11:01:31 UTC (rev 6241)
@@ -1,8 +1,12 @@
 # Author: Travis Oliphant
 # 2003
+#
+# Feb. 2010: Updated by Warren Weckesser:
+#   Rewrote much of chirp()
+#   Added sweep_poly()
 
 from numpy import asarray, zeros, place, nan, mod, pi, extract, log, sqrt, \
-     exp, cos, sin, size, polyval, polyint, log10
+     exp, cos, sin, polyval, polyint
 
 def sawtooth(t,width=1):
     """Returns a periodic sawtooth waveform with period 2*pi
@@ -135,67 +139,174 @@
     if retquad and retenv:
         return yI, yQ, yenv
 
-def chirp(t, f0=0, t1=1, f1=100, method='linear', phi=0, qshape=None):
+
+def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
     """Frequency-swept cosine generator.
 
+    In the following, 'Hz' should be interpreted as 'cycles per time unit'; there is
+    no assumption here that the time unit is one second.  The important distinction
+    is that the units of rotation are cycles, not radians.
+
     Parameters
     ----------
     t : ndarray
         Times at which to evaluate the waveform.
-    f0 : float or ndarray, optional
-        Frequency (in Hz) of the waveform at time 0.  If `f0` is an
-        ndarray, it specifies the frequency change as a polynomial in
-        `t` (see Notes below).
-    t1 : float, optional
+    f0 : float
+        Frequency (in Hz) at time t=0. 
+    t1 : float
         Time at which `f1` is specified.
-    f1 : float, optional
+    f1 : float
         Frequency (in Hz) of the waveform at time `t1`.
-    method : {'linear', 'quadratic', 'logarithmic'}, optional
-        Kind of frequency sweep.
-    phi : float
-        Phase offset, in degrees.
-    qshape : {'convex', 'concave'}
-        If method is 'quadratic', `qshape` specifies its shape.
+    method : {'linear', 'quadratic', 'logarithmic', 'hyperbolic'}, optional
+        Kind of frequency sweep.  If not given, `linear` is assumed.  See Notes below
+        for more details.
+    phi : float, optional
+        Phase offset, in degrees. Default is 0.
+    vertex_zero : bool, optional
+        This parameter is only used when `method` is 'quadratic'.
+        It determines whether the vertex of the parabola that is the graph of the
+        frequency is at t=0 or t=t1.
 
+    Returns
+    -------
+    A numpy array containing the signal evaluated at 't' with the requested
+    time-varying frequency.  More precisely, the function returns
+        cos(phase + (pi/180)*phi)
+    where `phase` is the integral (from 0 to t) of 2*pi*f(t).  f(t) is defined below.
+
     Notes
     -----
-    If `f0` is an array, it forms the coefficients of a polynomial in
-    `t` (see `numpy.polval`). The polynomial determines the waveform
-    frequency change in time.  In this case, the values of `f1`, `t1`,
-    `method`, and `qshape` are ignored.
+    There are four options for the `method`.  The following formulas give the
+    instantaneous frequency (in Hz) of the signal generated by `chirp()`.
+    For convenience, the shorter names shown below may also be used.
+    
+    linear, lin, li:
 
+        f(t) = f0 + (f1 - f0) * t / t1
+        
+    quadratic, quad, q:
+
+        The graph of the frequency f(t) is a parabola through (0, f0) and (t1, f1).
+        By default, the vertex of the parabola is at (0, f0).  If `vertex_zero`
+        is False, then the vertex is at (t1, f1).  The formula is:
+
+        if vertex_zero is True::
+            f(t) = f0 + (f1 - f0) * t**2 / t1**2
+        else:
+            f(t) = f1 - (f1 - f0) * (t1 - t)**2 / t1**2
+
+        To use a more general quadratic function, or an arbitrary polynomial,
+        use the function `scipy.signal.waveforms.sweep_poly`.
+
+    logarithmic, log, lo:
+
+        f(t) = f0 * (f1/f0)**(t/t1)
+        
+        f0 and f1 must be nonzero and have the same sign.
+
+        This signal is also known as a geometric or exponential chirp.
+
+    hyperbolic, hyp:
+
+        f(t) = f0*f1*t1 / ((f0 - f1)*t + f1*t1)
+        
+        f1 must be positive, and f0 must be greater than f1.
+
+    See Also
+    --------
+    scipy.signal.waveforms.sweep_poly
     """
-    # Convert to radians.
+    # 'phase' is computed in _chirp_phase, to make testing easier.
+    phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero)
+    # Convert  phi to radians.
     phi *= pi / 180
-    if size(f0) > 1:
-        # We were given a polynomial.
-        return cos(2*pi*polyval(polyint(f0),t)+phi)
-    if method in ['linear','lin','li']:
-        beta = (f1-f0)/t1
-        phase_angle = 2*pi * (f0*t + 0.5*beta*t*t)
+    return cos(phase + phi)
+
+    
+def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True):
+    """
+    Calculate the phase used by chirp_phase to generate its output.  See
+    chirp_phase for a description of the arguments.
+    """
+    if method in ['linear', 'lin', 'li']:
+        beta = (f1 - f0) / t1
+        phase = 2*pi * (f0*t + 0.5*beta*t*t)
+
     elif method in ['quadratic','quad','q']:
-        if qshape == 'concave':
-            mxf = max(f0,f1)
-            mnf = min(f0,f1)
-            f1,f0 = mxf, mnf
-        elif qshape == 'convex':
-            mxf = max(f0,f1)
-            mnf = min(f0,f1)
-            f1,f0 = mnf, mxf
+        beta = (f1 - f0)/(t1**2)
+        if vertex_zero:
+            phase = 2*pi * (f0*t + beta * t**3/3)
         else:
-            raise ValueError("qshape must be either 'concave' or 'convex' but "
-                "a value of %r was given." % qshape)
-        beta = (f1-f0)/t1/t1
-        phase_angle = 2*pi * (f0*t + beta*t*t*t/3)
-    elif method in ['logarithmic','log','lo']:
-        if f1 <= f0:
-            raise ValueError(
-                "For a logarithmic sweep, f1=%f must be larger than f0=%f."
-                % (f1, f0))
-        beta = log10(f1-f0)/t1
-        phase_angle = 2*pi * (f0*t + (pow(10,beta*t)-1)/(beta*log(10)))
+            phase = 2*pi * (f1*t + beta * ((t1 - t)**3 - t1**3)/3)               
+
+    elif method in ['logarithmic', 'log', 'lo']:
+        if f0*f1 <= 0.0:
+            raise ValueError("For a geometric chirp, f0 and f1 must be nonzero " \
+                                "and have the same sign.")
+        if f0 == f1:
+            phase = 2*pi * f0 * t
+        else:
+            beta = t1 / log(f1/f0)
+            phase = 2*pi * beta * f0 * (pow(f1/f0, t/t1) - 1.0)
+
+    elif method in ['hyperbolic', 'hyp']:
+        if f1 <= 0.0 or f0 <= f1:
+            raise ValueError("hyperbolic chirp requires f0 > f1 > 0.0.")
+        c = f1*t1
+        df = f0 - f1
+        phase = 2*pi * (f0 * c / df) * log((df*t + c)/c)
+
     else:
-        raise ValueError("method must be 'linear', 'quadratic', or "
-            "'logarithmic' but a value of %r was given." % method)
+        raise ValueError("method must be 'linear', 'quadratic', 'logarithmic', "
+                "or 'hyperbolic', but a value of %r was given." % method)
 
-    return cos(phase_angle + phi)
+    return phase
+
+
+def sweep_poly(t, poly, phi=0):
+    """Frequency-swept cosine generator, with a time-dependent frequency specified
+    as a polynomial.
+    
+    This function generates a sinusoidal function whose instantaneous frequency
+    varies with time.  The frequency at time `t` is given by the polynomial `poly`.
+
+    Parameters
+    ----------
+    t : ndarray
+        Times at which to evaluate the waveform.
+    poly : 1D ndarray (or array-like), or instance of numpy.poly1d
+        The desired frequency expressed as a polynomial.  If `poly` is a list or
+        ndarray of length n, then the elements of `poly` are the coefficients of the
+        polynomial, and the instantaneous frequency is
+            f(t) = poly[0]*t**(n-1) + poly[1]*t**(n-2) + ... + poly[n-1]
+        If `poly` is an instance of numpy.poly1d, then the instantaneous frequency is
+            f(t) = poly(t)
+    phi : float, optional
+        Phase offset, in degrees. Default is 0.
+
+    Returns
+    -------
+    A numpy array containing the signal evaluated at 't' with the requested
+    time-varying frequency.  More precisely, the function returns
+        cos(phase + (pi/180)*phi)
+    where `phase` is the integral (from 0 to t) of 2*pi*f(t); f(t) is defined above. 
+    
+    See Also
+    --------
+    scipy.signal.waveforms.chirp
+    """
+    # 'phase' is computed in _sweep_poly_phase, to make testing easier.
+    phase = _sweep_poly_phase(t, poly)
+    # Convert to radians.
+    phi *= pi / 180
+    return cos(phase + phi)
+
+def _sweep_poly_phase(t, poly):
+    """
+    Calculate the phase used by sweep_poly to generate its output.  See
+    sweep_poly for a description of the arguments.
+    """
+    # polyint handles lists, ndarrays and instances of poly1d automatically.
+    intpoly = polyint(poly)
+    phase = 2*pi * polyval(intpoly, t)
+    return phase



More information about the Scipy-svn mailing list