[Scipy-svn] r6241 - in trunk/scipy/signal: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Feb 19 06:01:32 EST 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