Coverage for /home/martinb/.local/share/virtualenvs/camcops/lib/python3.6/site-packages/scipy/signal/waveforms.py : 9%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1# Author: Travis Oliphant
2# 2003
3#
4# Feb. 2010: Updated by Warren Weckesser:
5# Rewrote much of chirp()
6# Added sweep_poly()
7import numpy as np
8from numpy import asarray, zeros, place, nan, mod, pi, extract, log, sqrt, \
9 exp, cos, sin, polyval, polyint
12__all__ = ['sawtooth', 'square', 'gausspulse', 'chirp', 'sweep_poly',
13 'unit_impulse']
16def sawtooth(t, width=1):
17 """
18 Return a periodic sawtooth or triangle waveform.
20 The sawtooth waveform has a period ``2*pi``, rises from -1 to 1 on the
21 interval 0 to ``width*2*pi``, then drops from 1 to -1 on the interval
22 ``width*2*pi`` to ``2*pi``. `width` must be in the interval [0, 1].
24 Note that this is not band-limited. It produces an infinite number
25 of harmonics, which are aliased back and forth across the frequency
26 spectrum.
28 Parameters
29 ----------
30 t : array_like
31 Time.
32 width : array_like, optional
33 Width of the rising ramp as a proportion of the total cycle.
34 Default is 1, producing a rising ramp, while 0 produces a falling
35 ramp. `width` = 0.5 produces a triangle wave.
36 If an array, causes wave shape to change over time, and must be the
37 same length as t.
39 Returns
40 -------
41 y : ndarray
42 Output array containing the sawtooth waveform.
44 Examples
45 --------
46 A 5 Hz waveform sampled at 500 Hz for 1 second:
48 >>> from scipy import signal
49 >>> import matplotlib.pyplot as plt
50 >>> t = np.linspace(0, 1, 500)
51 >>> plt.plot(t, signal.sawtooth(2 * np.pi * 5 * t))
53 """
54 t, w = asarray(t), asarray(width)
55 w = asarray(w + (t - t))
56 t = asarray(t + (w - w))
57 if t.dtype.char in ['fFdD']:
58 ytype = t.dtype.char
59 else:
60 ytype = 'd'
61 y = zeros(t.shape, ytype)
63 # width must be between 0 and 1 inclusive
64 mask1 = (w > 1) | (w < 0)
65 place(y, mask1, nan)
67 # take t modulo 2*pi
68 tmod = mod(t, 2 * pi)
70 # on the interval 0 to width*2*pi function is
71 # tmod / (pi*w) - 1
72 mask2 = (1 - mask1) & (tmod < w * 2 * pi)
73 tsub = extract(mask2, tmod)
74 wsub = extract(mask2, w)
75 place(y, mask2, tsub / (pi * wsub) - 1)
77 # on the interval width*2*pi to 2*pi function is
78 # (pi*(w+1)-tmod) / (pi*(1-w))
80 mask3 = (1 - mask1) & (1 - mask2)
81 tsub = extract(mask3, tmod)
82 wsub = extract(mask3, w)
83 place(y, mask3, (pi * (wsub + 1) - tsub) / (pi * (1 - wsub)))
84 return y
87def square(t, duty=0.5):
88 """
89 Return a periodic square-wave waveform.
91 The square wave has a period ``2*pi``, has value +1 from 0 to
92 ``2*pi*duty`` and -1 from ``2*pi*duty`` to ``2*pi``. `duty` must be in
93 the interval [0,1].
95 Note that this is not band-limited. It produces an infinite number
96 of harmonics, which are aliased back and forth across the frequency
97 spectrum.
99 Parameters
100 ----------
101 t : array_like
102 The input time array.
103 duty : array_like, optional
104 Duty cycle. Default is 0.5 (50% duty cycle).
105 If an array, causes wave shape to change over time, and must be the
106 same length as t.
108 Returns
109 -------
110 y : ndarray
111 Output array containing the square waveform.
113 Examples
114 --------
115 A 5 Hz waveform sampled at 500 Hz for 1 second:
117 >>> from scipy import signal
118 >>> import matplotlib.pyplot as plt
119 >>> t = np.linspace(0, 1, 500, endpoint=False)
120 >>> plt.plot(t, signal.square(2 * np.pi * 5 * t))
121 >>> plt.ylim(-2, 2)
123 A pulse-width modulated sine wave:
125 >>> plt.figure()
126 >>> sig = np.sin(2 * np.pi * t)
127 >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
128 >>> plt.subplot(2, 1, 1)
129 >>> plt.plot(t, sig)
130 >>> plt.subplot(2, 1, 2)
131 >>> plt.plot(t, pwm)
132 >>> plt.ylim(-1.5, 1.5)
134 """
135 t, w = asarray(t), asarray(duty)
136 w = asarray(w + (t - t))
137 t = asarray(t + (w - w))
138 if t.dtype.char in ['fFdD']:
139 ytype = t.dtype.char
140 else:
141 ytype = 'd'
143 y = zeros(t.shape, ytype)
145 # width must be between 0 and 1 inclusive
146 mask1 = (w > 1) | (w < 0)
147 place(y, mask1, nan)
149 # on the interval 0 to duty*2*pi function is 1
150 tmod = mod(t, 2 * pi)
151 mask2 = (1 - mask1) & (tmod < w * 2 * pi)
152 place(y, mask2, 1)
154 # on the interval duty*2*pi to 2*pi function is
155 # (pi*(w+1)-tmod) / (pi*(1-w))
156 mask3 = (1 - mask1) & (1 - mask2)
157 place(y, mask3, -1)
158 return y
161def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
162 retenv=False):
163 """
164 Return a Gaussian modulated sinusoid:
166 ``exp(-a t^2) exp(1j*2*pi*fc*t).``
168 If `retquad` is True, then return the real and imaginary parts
169 (in-phase and quadrature).
170 If `retenv` is True, then return the envelope (unmodulated signal).
171 Otherwise, return the real part of the modulated sinusoid.
173 Parameters
174 ----------
175 t : ndarray or the string 'cutoff'
176 Input array.
177 fc : float, optional
178 Center frequency (e.g. Hz). Default is 1000.
179 bw : float, optional
180 Fractional bandwidth in frequency domain of pulse (e.g. Hz).
181 Default is 0.5.
182 bwr : float, optional
183 Reference level at which fractional bandwidth is calculated (dB).
184 Default is -6.
185 tpr : float, optional
186 If `t` is 'cutoff', then the function returns the cutoff
187 time for when the pulse amplitude falls below `tpr` (in dB).
188 Default is -60.
189 retquad : bool, optional
190 If True, return the quadrature (imaginary) as well as the real part
191 of the signal. Default is False.
192 retenv : bool, optional
193 If True, return the envelope of the signal. Default is False.
195 Returns
196 -------
197 yI : ndarray
198 Real part of signal. Always returned.
199 yQ : ndarray
200 Imaginary part of signal. Only returned if `retquad` is True.
201 yenv : ndarray
202 Envelope of signal. Only returned if `retenv` is True.
204 See Also
205 --------
206 scipy.signal.morlet
208 Examples
209 --------
210 Plot real component, imaginary component, and envelope for a 5 Hz pulse,
211 sampled at 100 Hz for 2 seconds:
213 >>> from scipy import signal
214 >>> import matplotlib.pyplot as plt
215 >>> t = np.linspace(-1, 1, 2 * 100, endpoint=False)
216 >>> i, q, e = signal.gausspulse(t, fc=5, retquad=True, retenv=True)
217 >>> plt.plot(t, i, t, q, t, e, '--')
219 """
220 if fc < 0:
221 raise ValueError("Center frequency (fc=%.2f) must be >=0." % fc)
222 if bw <= 0:
223 raise ValueError("Fractional bandwidth (bw=%.2f) must be > 0." % bw)
224 if bwr >= 0:
225 raise ValueError("Reference level for bandwidth (bwr=%.2f) must "
226 "be < 0 dB" % bwr)
228 # exp(-a t^2) <-> sqrt(pi/a) exp(-pi^2/a * f^2) = g(f)
230 ref = pow(10.0, bwr / 20.0)
231 # fdel = fc*bw/2: g(fdel) = ref --- solve this for a
232 #
233 # pi^2/a * fc^2 * bw^2 /4=-log(ref)
234 a = -(pi * fc * bw) ** 2 / (4.0 * log(ref))
236 if isinstance(t, str):
237 if t == 'cutoff': # compute cut_off point
238 # Solve exp(-a tc**2) = tref for tc
239 # tc = sqrt(-log(tref) / a) where tref = 10^(tpr/20)
240 if tpr >= 0:
241 raise ValueError("Reference level for time cutoff must "
242 "be < 0 dB")
243 tref = pow(10.0, tpr / 20.0)
244 return sqrt(-log(tref) / a)
245 else:
246 raise ValueError("If `t` is a string, it must be 'cutoff'")
248 yenv = exp(-a * t * t)
249 yI = yenv * cos(2 * pi * fc * t)
250 yQ = yenv * sin(2 * pi * fc * t)
251 if not retquad and not retenv:
252 return yI
253 if not retquad and retenv:
254 return yI, yenv
255 if retquad and not retenv:
256 return yI, yQ
257 if retquad and retenv:
258 return yI, yQ, yenv
261def chirp(t, f0, t1, f1, method='linear', phi=0, vertex_zero=True):
262 """Frequency-swept cosine generator.
264 In the following, 'Hz' should be interpreted as 'cycles per unit';
265 there is no requirement here that the unit is one second. The
266 important distinction is that the units of rotation are cycles, not
267 radians. Likewise, `t` could be a measurement of space instead of time.
269 Parameters
270 ----------
271 t : array_like
272 Times at which to evaluate the waveform.
273 f0 : float
274 Frequency (e.g. Hz) at time t=0.
275 t1 : float
276 Time at which `f1` is specified.
277 f1 : float
278 Frequency (e.g. Hz) of the waveform at time `t1`.
279 method : {'linear', 'quadratic', 'logarithmic', 'hyperbolic'}, optional
280 Kind of frequency sweep. If not given, `linear` is assumed. See
281 Notes below for more details.
282 phi : float, optional
283 Phase offset, in degrees. Default is 0.
284 vertex_zero : bool, optional
285 This parameter is only used when `method` is 'quadratic'.
286 It determines whether the vertex of the parabola that is the graph
287 of the frequency is at t=0 or t=t1.
289 Returns
290 -------
291 y : ndarray
292 A numpy array containing the signal evaluated at `t` with the
293 requested time-varying frequency. More precisely, the function
294 returns ``cos(phase + (pi/180)*phi)`` where `phase` is the integral
295 (from 0 to `t`) of ``2*pi*f(t)``. ``f(t)`` is defined below.
297 See Also
298 --------
299 sweep_poly
301 Notes
302 -----
303 There are four options for the `method`. The following formulas give
304 the instantaneous frequency (in Hz) of the signal generated by
305 `chirp()`. For convenience, the shorter names shown below may also be
306 used.
308 linear, lin, li:
310 ``f(t) = f0 + (f1 - f0) * t / t1``
312 quadratic, quad, q:
314 The graph of the frequency f(t) is a parabola through (0, f0) and
315 (t1, f1). By default, the vertex of the parabola is at (0, f0).
316 If `vertex_zero` is False, then the vertex is at (t1, f1). The
317 formula is:
319 if vertex_zero is True:
321 ``f(t) = f0 + (f1 - f0) * t**2 / t1**2``
323 else:
325 ``f(t) = f1 - (f1 - f0) * (t1 - t)**2 / t1**2``
327 To use a more general quadratic function, or an arbitrary
328 polynomial, use the function `scipy.signal.sweep_poly`.
330 logarithmic, log, lo:
332 ``f(t) = f0 * (f1/f0)**(t/t1)``
334 f0 and f1 must be nonzero and have the same sign.
336 This signal is also known as a geometric or exponential chirp.
338 hyperbolic, hyp:
340 ``f(t) = f0*f1*t1 / ((f0 - f1)*t + f1*t1)``
342 f0 and f1 must be nonzero.
344 Examples
345 --------
346 The following will be used in the examples:
348 >>> from scipy.signal import chirp, spectrogram
349 >>> import matplotlib.pyplot as plt
351 For the first example, we'll plot the waveform for a linear chirp
352 from 6 Hz to 1 Hz over 10 seconds:
354 >>> t = np.linspace(0, 10, 1500)
355 >>> w = chirp(t, f0=6, f1=1, t1=10, method='linear')
356 >>> plt.plot(t, w)
357 >>> plt.title("Linear Chirp, f(0)=6, f(10)=1")
358 >>> plt.xlabel('t (sec)')
359 >>> plt.show()
361 For the remaining examples, we'll use higher frequency ranges,
362 and demonstrate the result using `scipy.signal.spectrogram`.
363 We'll use a 4 second interval sampled at 7200 Hz.
365 >>> fs = 7200
366 >>> T = 4
367 >>> t = np.arange(0, int(T*fs)) / fs
369 We'll use this function to plot the spectrogram in each example.
371 >>> def plot_spectrogram(title, w, fs):
372 ... ff, tt, Sxx = spectrogram(w, fs=fs, nperseg=256, nfft=576)
373 ... plt.pcolormesh(tt, ff[:145], Sxx[:145], cmap='gray_r', shading='gouraud')
374 ... plt.title(title)
375 ... plt.xlabel('t (sec)')
376 ... plt.ylabel('Frequency (Hz)')
377 ... plt.grid()
378 ...
380 Quadratic chirp from 1500 Hz to 250 Hz
381 (vertex of the parabolic curve of the frequency is at t=0):
383 >>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic')
384 >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250', w, fs)
385 >>> plt.show()
387 Quadratic chirp from 1500 Hz to 250 Hz
388 (vertex of the parabolic curve of the frequency is at t=T):
390 >>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic',
391 ... vertex_zero=False)
392 >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f({T})=250\\n' +
393 ... '(vertex_zero=False)', w, fs)
394 >>> plt.show()
396 Logarithmic chirp from 1500 Hz to 250 Hz:
398 >>> w = chirp(t, f0=1500, f1=250, t1=T, method='logarithmic')
399 >>> plot_spectrogram(f'Logarithmic Chirp, f(0)=1500, f({T})=250', w, fs)
400 >>> plt.show()
402 Hyperbolic chirp from 1500 Hz to 250 Hz:
404 >>> w = chirp(t, f0=1500, f1=250, t1=T, method='hyperbolic')
405 >>> plot_spectrogram(f'Hyperbolic Chirp, f(0)=1500, f({T})=250', w, fs)
406 >>> plt.show()
408 """
409 # 'phase' is computed in _chirp_phase, to make testing easier.
410 phase = _chirp_phase(t, f0, t1, f1, method, vertex_zero)
411 # Convert phi to radians.
412 phi *= pi / 180
413 return cos(phase + phi)
416def _chirp_phase(t, f0, t1, f1, method='linear', vertex_zero=True):
417 """
418 Calculate the phase used by `chirp` to generate its output.
420 See `chirp` for a description of the arguments.
422 """
423 t = asarray(t)
424 f0 = float(f0)
425 t1 = float(t1)
426 f1 = float(f1)
427 if method in ['linear', 'lin', 'li']:
428 beta = (f1 - f0) / t1
429 phase = 2 * pi * (f0 * t + 0.5 * beta * t * t)
431 elif method in ['quadratic', 'quad', 'q']:
432 beta = (f1 - f0) / (t1 ** 2)
433 if vertex_zero:
434 phase = 2 * pi * (f0 * t + beta * t ** 3 / 3)
435 else:
436 phase = 2 * pi * (f1 * t + beta * ((t1 - t) ** 3 - t1 ** 3) / 3)
438 elif method in ['logarithmic', 'log', 'lo']:
439 if f0 * f1 <= 0.0:
440 raise ValueError("For a logarithmic chirp, f0 and f1 must be "
441 "nonzero and have the same sign.")
442 if f0 == f1:
443 phase = 2 * pi * f0 * t
444 else:
445 beta = t1 / log(f1 / f0)
446 phase = 2 * pi * beta * f0 * (pow(f1 / f0, t / t1) - 1.0)
448 elif method in ['hyperbolic', 'hyp']:
449 if f0 == 0 or f1 == 0:
450 raise ValueError("For a hyperbolic chirp, f0 and f1 must be "
451 "nonzero.")
452 if f0 == f1:
453 # Degenerate case: constant frequency.
454 phase = 2 * pi * f0 * t
455 else:
456 # Singular point: the instantaneous frequency blows up
457 # when t == sing.
458 sing = -f1 * t1 / (f0 - f1)
459 phase = 2 * pi * (-sing * f0) * log(np.abs(1 - t/sing))
461 else:
462 raise ValueError("method must be 'linear', 'quadratic', 'logarithmic',"
463 " or 'hyperbolic', but a value of %r was given."
464 % method)
466 return phase
469def sweep_poly(t, poly, phi=0):
470 """
471 Frequency-swept cosine generator, with a time-dependent frequency.
473 This function generates a sinusoidal function whose instantaneous
474 frequency varies with time. The frequency at time `t` is given by
475 the polynomial `poly`.
477 Parameters
478 ----------
479 t : ndarray
480 Times at which to evaluate the waveform.
481 poly : 1-D array_like or instance of numpy.poly1d
482 The desired frequency expressed as a polynomial. If `poly` is
483 a list or ndarray of length n, then the elements of `poly` are
484 the coefficients of the polynomial, and the instantaneous
485 frequency is
487 ``f(t) = poly[0]*t**(n-1) + poly[1]*t**(n-2) + ... + poly[n-1]``
489 If `poly` is an instance of numpy.poly1d, then the
490 instantaneous frequency is
492 ``f(t) = poly(t)``
494 phi : float, optional
495 Phase offset, in degrees, Default: 0.
497 Returns
498 -------
499 sweep_poly : ndarray
500 A numpy array containing the signal evaluated at `t` with the
501 requested time-varying frequency. More precisely, the function
502 returns ``cos(phase + (pi/180)*phi)``, where `phase` is the integral
503 (from 0 to t) of ``2 * pi * f(t)``; ``f(t)`` is defined above.
505 See Also
506 --------
507 chirp
509 Notes
510 -----
511 .. versionadded:: 0.8.0
513 If `poly` is a list or ndarray of length `n`, then the elements of
514 `poly` are the coefficients of the polynomial, and the instantaneous
515 frequency is:
517 ``f(t) = poly[0]*t**(n-1) + poly[1]*t**(n-2) + ... + poly[n-1]``
519 If `poly` is an instance of `numpy.poly1d`, then the instantaneous
520 frequency is:
522 ``f(t) = poly(t)``
524 Finally, the output `s` is:
526 ``cos(phase + (pi/180)*phi)``
528 where `phase` is the integral from 0 to `t` of ``2 * pi * f(t)``,
529 ``f(t)`` as defined above.
531 Examples
532 --------
533 Compute the waveform with instantaneous frequency::
535 f(t) = 0.025*t**3 - 0.36*t**2 + 1.25*t + 2
537 over the interval 0 <= t <= 10.
539 >>> from scipy.signal import sweep_poly
540 >>> p = np.poly1d([0.025, -0.36, 1.25, 2.0])
541 >>> t = np.linspace(0, 10, 5001)
542 >>> w = sweep_poly(t, p)
544 Plot it:
546 >>> import matplotlib.pyplot as plt
547 >>> plt.subplot(2, 1, 1)
548 >>> plt.plot(t, w)
549 >>> plt.title("Sweep Poly\\nwith frequency " +
550 ... "$f(t) = 0.025t^3 - 0.36t^2 + 1.25t + 2$")
551 >>> plt.subplot(2, 1, 2)
552 >>> plt.plot(t, p(t), 'r', label='f(t)')
553 >>> plt.legend()
554 >>> plt.xlabel('t')
555 >>> plt.tight_layout()
556 >>> plt.show()
558 """
559 # 'phase' is computed in _sweep_poly_phase, to make testing easier.
560 phase = _sweep_poly_phase(t, poly)
561 # Convert to radians.
562 phi *= pi / 180
563 return cos(phase + phi)
566def _sweep_poly_phase(t, poly):
567 """
568 Calculate the phase used by sweep_poly to generate its output.
570 See `sweep_poly` for a description of the arguments.
572 """
573 # polyint handles lists, ndarrays and instances of poly1d automatically.
574 intpoly = polyint(poly)
575 phase = 2 * pi * polyval(intpoly, t)
576 return phase
579def unit_impulse(shape, idx=None, dtype=float):
580 """
581 Unit impulse signal (discrete delta function) or unit basis vector.
583 Parameters
584 ----------
585 shape : int or tuple of int
586 Number of samples in the output (1-D), or a tuple that represents the
587 shape of the output (N-D).
588 idx : None or int or tuple of int or 'mid', optional
589 Index at which the value is 1. If None, defaults to the 0th element.
590 If ``idx='mid'``, the impulse will be centered at ``shape // 2`` in
591 all dimensions. If an int, the impulse will be at `idx` in all
592 dimensions.
593 dtype : data-type, optional
594 The desired data-type for the array, e.g., ``numpy.int8``. Default is
595 ``numpy.float64``.
597 Returns
598 -------
599 y : ndarray
600 Output array containing an impulse signal.
602 Notes
603 -----
604 The 1D case is also known as the Kronecker delta.
606 .. versionadded:: 0.19.0
608 Examples
609 --------
610 An impulse at the 0th element (:math:`\\delta[n]`):
612 >>> from scipy import signal
613 >>> signal.unit_impulse(8)
614 array([ 1., 0., 0., 0., 0., 0., 0., 0.])
616 Impulse offset by 2 samples (:math:`\\delta[n-2]`):
618 >>> signal.unit_impulse(7, 2)
619 array([ 0., 0., 1., 0., 0., 0., 0.])
621 2-dimensional impulse, centered:
623 >>> signal.unit_impulse((3, 3), 'mid')
624 array([[ 0., 0., 0.],
625 [ 0., 1., 0.],
626 [ 0., 0., 0.]])
628 Impulse at (2, 2), using broadcasting:
630 >>> signal.unit_impulse((4, 4), 2)
631 array([[ 0., 0., 0., 0.],
632 [ 0., 0., 0., 0.],
633 [ 0., 0., 1., 0.],
634 [ 0., 0., 0., 0.]])
636 Plot the impulse response of a 4th-order Butterworth lowpass filter:
638 >>> imp = signal.unit_impulse(100, 'mid')
639 >>> b, a = signal.butter(4, 0.2)
640 >>> response = signal.lfilter(b, a, imp)
642 >>> import matplotlib.pyplot as plt
643 >>> plt.plot(np.arange(-50, 50), imp)
644 >>> plt.plot(np.arange(-50, 50), response)
645 >>> plt.margins(0.1, 0.1)
646 >>> plt.xlabel('Time [samples]')
647 >>> plt.ylabel('Amplitude')
648 >>> plt.grid(True)
649 >>> plt.show()
651 """
652 out = zeros(shape, dtype)
654 shape = np.atleast_1d(shape)
656 if idx is None:
657 idx = (0,) * len(shape)
658 elif idx == 'mid':
659 idx = tuple(shape // 2)
660 elif not hasattr(idx, "__iter__"):
661 idx = (idx,) * len(shape)
663 out[idx] = 1
664 return out