package scipy

Get an attribute of this module as a Py.Object.t. This is useful to pass a Python function to another function.

Check state-space matrices and ensure they are 2-D.

If enough information on the system is provided, that is, enough properly-shaped arrays are passed to the function, the missing ones are built from this information, ensuring the correct number of rows and columns. Otherwise a ValueError is raised.

Parameters ---------- A, B, C, D : array_like, optional State-space matrices. All of them are None (missing) by default. See `ss2tf` for format.

Returns ------- A, B, C, D : array Properly shaped state-space matrices.

Raises ------ ValueError If not enough information on the system was provided.

Calculate the relative extrema of `data`.

Parameters ---------- data : ndarray Array in which to find the relative extrema. comparator : callable Function to use to compare two data points. Should take two arrays as arguments. axis : int, optional Axis over which to select from `data`. Default is 0. order : int, optional How many points on each side to use for the comparison to consider ``comparator(n, n+x)`` to be True. mode : str, optional How the edges of the vector are treated. 'wrap' (wrap around) or 'clip' (treat overflow as the same as the last (or first) element). Default is 'clip'. See `numpy.take`.

Returns ------- extrema : tuple of ndarrays Indices of the maxima in arrays of integers. ``extremak`` is the array of indices of axis `k` of `data`. Note that the return value is a tuple even when `data` is 1-D.

See Also -------- argrelmin, argrelmax

Notes -----

.. versionadded:: 0.11.0

Examples -------- >>> from scipy.signal import argrelextrema >>> x = np.array(2, 1, 2, 3, 2, 0, 1, 0) >>> argrelextrema(x, np.greater) (array(3, 6),) >>> y = np.array([1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]) ... >>> argrelextrema(y, np.less, axis=1) (array(0, 2), array(2, 1))

Calculate the relative maxima of `data`.

Parameters ---------- data : ndarray Array in which to find the relative maxima. axis : int, optional Axis over which to select from `data`. Default is 0. order : int, optional How many points on each side to use for the comparison to consider ``comparator(n, n+x)`` to be True. mode : str, optional How the edges of the vector are treated. Available options are 'wrap' (wrap around) or 'clip' (treat overflow as the same as the last (or first) element). Default 'clip'. See `numpy.take`.

Returns ------- extrema : tuple of ndarrays Indices of the maxima in arrays of integers. ``extremak`` is the array of indices of axis `k` of `data`. Note that the return value is a tuple even when `data` is 1-D.

See Also -------- argrelextrema, argrelmin, find_peaks

Notes ----- This function uses `argrelextrema` with np.greater as comparator. Therefore, it requires a strict inequality on both sides of a value to consider it a maximum. This means flat maxima (more than one sample wide) are not detected. In case of 1-D `data` `find_peaks` can be used to detect all local maxima, including flat ones.

.. versionadded:: 0.11.0

Examples -------- >>> from scipy.signal import argrelmax >>> x = np.array(2, 1, 2, 3, 2, 0, 1, 0) >>> argrelmax(x) (array(3, 6),) >>> y = np.array([1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]) ... >>> argrelmax(y, axis=1) (array(0), array(1))

Calculate the relative minima of `data`.

Parameters ---------- data : ndarray Array in which to find the relative minima. axis : int, optional Axis over which to select from `data`. Default is 0. order : int, optional How many points on each side to use for the comparison to consider ``comparator(n, n+x)`` to be True. mode : str, optional How the edges of the vector are treated. Available options are 'wrap' (wrap around) or 'clip' (treat overflow as the same as the last (or first) element). Default 'clip'. See numpy.take.

Returns ------- extrema : tuple of ndarrays Indices of the minima in arrays of integers. ``extremak`` is the array of indices of axis `k` of `data`. Note that the return value is a tuple even when `data` is 1-D.

See Also -------- argrelextrema, argrelmax, find_peaks

Notes ----- This function uses `argrelextrema` with np.less as comparator. Therefore, it requires a strict inequality on both sides of a value to consider it a minimum. This means flat minima (more than one sample wide) are not detected. In case of 1-D `data` `find_peaks` can be used to detect all local minima, including flat ones, by calling it with negated `data`.

.. versionadded:: 0.11.0

Examples -------- >>> from scipy.signal import argrelmin >>> x = np.array(2, 1, 2, 3, 2, 0, 1, 0) >>> argrelmin(x) (array(1, 5),) >>> y = np.array([1, 2, 1, 2], ... [2, 2, 0, 0], ... [5, 3, 4, 4]) ... >>> argrelmin(y, axis=1) (array(0, 2), array(2, 1))

Band Stop Objective Function for order minimization.

Returns the non-integer order for an analog band stop filter.

Parameters ---------- wp : scalar Edge of passband `passb`. ind : int,

, 1

Index specifying which `passb` edge to vary (0 or 1). passb : ndarray Two element sequence of fixed passband edges. stopb : ndarray Two element sequence of fixed stopband edges. gstop : float Amount of attenuation in stopband in dB. gpass : float Amount of ripple in the passband in dB. type : 'butter', 'cheby', 'ellip' Type of filter.

Returns ------- n : scalar Filter order (possibly non-integer).

Return a modified Bartlett-Hann window.

.. warning:: scipy.signal.barthann is deprecated, use scipy.signal.windows.barthann instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.barthann(51) >>> plt.plot(window) >>> plt.title('Bartlett-Hann window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Bartlett-Hann window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Return a Bartlett window.

The Bartlett window is very similar to a triangular window, except that the end points are at zero. It is often used in signal processing for tapering a signal, without generating too much ripple in the frequency domain.

.. warning:: scipy.signal.bartlett is deprecated, use scipy.signal.windows.bartlett instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The triangular window, with the first and last samples equal to zero and the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

See Also -------- triang : A triangular window that does not touch zero at the ends

Notes ----- The Bartlett window is defined as

.. math:: w(n) = \frac

M-1 \left( \fracM-1

• \left|n - \fracM-1

  \right| \right)

Most references to the Bartlett window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. Note that convolution with this window produces linear interpolation. It is also known as an apodization (which means'removing the foot', i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. The Fourier transform of the Bartlett is the product of two sinc functions. Note the excellent discussion in Kanasewich. 2_

References ---------- .. 1 M.S. Bartlett, 'Periodogram Analysis and Continuous Spectra', Biometrika 37, 1-16, 1950. .. 2 E.R. Kanasewich, 'Time Sequence Analysis in Geophysics', The University of Alberta Press, 1975, pp. 109-110. .. 3 A.V. Oppenheim and R.W. Schafer, 'Discrete-Time Signal Processing', Prentice-Hall, 1999, pp. 468-471. .. 4 Wikipedia, 'Window function', https://en.wikipedia.org/wiki/Window_function .. 5 W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, 'Numerical Recipes', Cambridge University Press, 1986, page 429.

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.bartlett(51) >>> plt.plot(window) >>> plt.title('Bartlett window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Bartlett window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Bessel/Thomson digital and analog filter design.

Design an Nth-order digital or analog Bessel filter and return the filter coefficients.

Parameters ---------- N : int The order of the filter. Wn : array_like A scalar or length-2 sequence giving the critical frequencies (defined by the `norm` parameter). For analog filters, `Wn` is an angular frequency (e.g., rad/s).

For digital filters, `Wn` are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`Wn` is thus in half-cycles / sample.) btype : 'lowpass', 'highpass', 'bandpass', 'bandstop', optional The type of filter. Default is 'lowpass'. analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. (See Notes.) output : 'ba', 'zpk', 'sos', optional Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or second-order sections ('sos'). Default is 'ba'. norm : 'phase', 'delay', 'mag', optional Critical frequency normalization:

``phase`` The filter is normalized such that the phase response reaches its midpoint at angular (e.g. rad/s) frequency `Wn`. This happens for both low-pass and high-pass filters, so this is the 'phase-matched' case.

The magnitude response asymptotes are the same as a Butterworth filter of the same order with a cutoff of `Wn`.

This is the default, and matches MATLAB's implementation.

``delay`` The filter is normalized such that the group delay in the passband is 1/`Wn` (e.g., seconds). This is the 'natural' type obtained by solving Bessel polynomials.

``mag`` The filter is normalized such that the gain magnitude is -3 dB at angular frequency `Wn`.

.. versionadded:: 0.18.0 fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- b, a : ndarray, ndarray Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. Only returned if ``output='ba'``. z, p, k : ndarray, ndarray, float Zeros, poles, and system gain of the IIR filter transfer function. Only returned if ``output='zpk'``. sos : ndarray Second-order sections representation of the IIR filter. Only returned if ``output=='sos'``.

Notes ----- Also known as a Thomson filter, the analog Bessel filter has maximally flat group delay and maximally linear phase response, with very little ringing in the step response. 1_

The Bessel is inherently an analog filter. This function generates digital Bessel filters using the bilinear transform, which does not preserve the phase response of the analog filter. As such, it is only approximately correct at frequencies below about fs/4. To get maximally-flat group delay at higher frequencies, the analog Bessel filter must be transformed using phase-preserving techniques.

See `besselap` for implementation details and references.

The ``'sos'`` output parameter was added in 0.16.0.

Examples -------- Plot the phase-normalized frequency response, showing the relationship to the Butterworth's cutoff frequency (green):

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> b, a = signal.butter(4, 100, 'low', analog=True) >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(np.abs(h)), color='silver', ls='dashed') >>> b, a = signal.bessel(4, 100, 'low', analog=True, norm='phase') >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(np.abs(h))) >>> plt.title('Bessel filter magnitude response (with Butterworth)') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.show()

and the phase midpoint:

>>> plt.figure() >>> plt.semilogx(w, np.unwrap(np.angle(h))) >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.axhline(-np.pi, color='red') # phase midpoint >>> plt.title('Bessel filter phase response') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Phase radians') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.show()

Plot the magnitude-normalized frequency response, showing the -3 dB cutoff:

>>> b, a = signal.bessel(3, 10, 'low', analog=True, norm='mag') >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(np.abs(h))) >>> plt.axhline(-3, color='red') # -3 dB magnitude >>> plt.axvline(10, color='green') # cutoff frequency >>> plt.title('Magnitude-normalized Bessel filter frequency response') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.show()

Plot the delay-normalized filter, showing the maximally-flat group delay at 0.1 seconds:

>>> b, a = signal.bessel(5, 1/0.1, 'low', analog=True, norm='delay') >>> w, h = signal.freqs(b, a) >>> plt.figure() >>> plt.semilogx(w1:, -np.diff(np.unwrap(np.angle(h)))/np.diff(w)) >>> plt.axhline(0.1, color='red') # 0.1 seconds group delay >>> plt.title('Bessel filter group delay') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Group delay seconds') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.show()

References ---------- .. 1 Thomson, W.E., 'Delay Networks having Maximally Flat Frequency Characteristics', Proceedings of the Institution of Electrical Engineers, Part III, November 1949, Vol. 96, No. 44, pp. 487-490.

Return (z,p,k) for analog prototype of an Nth-order Bessel filter.

Parameters ---------- N : int The order of the filter. norm : 'phase', 'delay', 'mag', optional Frequency normalization:

``phase`` The filter is normalized such that the phase response reaches its midpoint at an angular (e.g., rad/s) cutoff frequency of 1. This happens for both low-pass and high-pass filters, so this is the 'phase-matched' case. 6_

The magnitude response asymptotes are the same as a Butterworth filter of the same order with a cutoff of `Wn`.

This is the default, and matches MATLAB's implementation.

``delay`` The filter is normalized such that the group delay in the passband is 1 (e.g., 1 second). This is the 'natural' type obtained by solving Bessel polynomials

``mag`` The filter is normalized such that the gain magnitude is -3 dB at angular frequency 1. This is called 'frequency normalization' by Bond. 1_

.. versionadded:: 0.18.0

Returns ------- z : ndarray Zeros of the transfer function. Is always an empty array. p : ndarray Poles of the transfer function. k : scalar Gain of the transfer function. For phase-normalized, this is always 1.

See Also -------- bessel : Filter design function using this prototype

Notes ----- To find the pole locations, approximate starting points are generated 2_ for the zeros of the ordinary Bessel polynomial 3_, then the Aberth-Ehrlich method 4_ 5_ is used on the Kv(x) Bessel function to calculate more accurate zeros, and these locations are then inverted about the unit circle.

References ---------- .. 1 C.R. Bond, 'Bessel Filter Constants', http://www.crbond.com/papers/bsf.pdf .. 2 Campos and Calderon, 'Approximate closed-form formulas for the zeros of the Bessel Polynomials', :arXiv:`1105.0957`. .. 3 Thomson, W.E., 'Delay Networks having Maximally Flat Frequency Characteristics', Proceedings of the Institution of Electrical Engineers, Part III, November 1949, Vol. 96, No. 44, pp. 487-490. .. 4 Aberth, 'Iteration Methods for Finding all Zeros of a Polynomial Simultaneously', Mathematics of Computation, Vol. 27, No. 122, April 1973 .. 5 Ehrlich, 'A modified Newton method for polynomials', Communications of the ACM, Vol. 10, Issue 2, pp. 107-108, Feb. 1967, :DOI:`10.1145/363067.363115` .. 6 Miller and Bohn, 'A Bessel Filter Crossover, and Its Relation to Others', RaneNote 147, 1998, https://www.ranecommercial.com/legacy/note147.html

Return a digital IIR filter from an analog one using a bilinear transform.

Transform a set of poles and zeros from the analog s-plane to the digital z-plane using Tustin's method, which substitutes ``(z-1) / (z+1)`` for ``s``, maintaining the shape of the frequency response.

Parameters ---------- b : array_like Numerator of the analog filter transfer function. a : array_like Denominator of the analog filter transfer function. fs : float Sample rate, as ordinary frequency (e.g., hertz). No prewarping is done in this function.

Returns ------- z : ndarray Numerator of the transformed digital filter transfer function. p : ndarray Denominator of the transformed digital filter transfer function.

See Also -------- lp2lp, lp2hp, lp2bp, lp2bs bilinear_zpk

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> fs = 100 >>> bf = 2 * np.pi * np.array(7, 13) >>> filts = signal.lti( *signal.butter(4, bf, btype='bandpass', analog=True)) >>> filtz = signal.lti( *signal.bilinear(filts.num, filts.den, fs)) >>> wz, hz = signal.freqz(filtz.num, filtz.den) >>> ws, hs = signal.freqs(filts.num, filts.den, worN=fs*wz)

>>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hz).clip(1e-15)), label=r'$|H(j \omega)|$') >>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hs).clip(1e-15)), label=r'$|H_z(e^j \omega)|$') >>> plt.legend() >>> plt.xlabel('Frequency Hz') >>> plt.ylabel('Magnitude dB') >>> plt.grid()

Return a digital IIR filter from an analog one using a bilinear transform.

Transform a set of poles and zeros from the analog s-plane to the digital z-plane using Tustin's method, which substitutes ``(z-1) / (z+1)`` for ``s``, maintaining the shape of the frequency response.

Parameters ---------- z : array_like Zeros of the analog filter transfer function. p : array_like Poles of the analog filter transfer function. k : float System gain of the analog filter transfer function. fs : float Sample rate, as ordinary frequency (e.g., hertz). No prewarping is done in this function.

Returns ------- z : ndarray Zeros of the transformed digital filter transfer function. p : ndarray Poles of the transformed digital filter transfer function. k : float System gain of the transformed digital filter.

See Also -------- lp2lp_zpk, lp2hp_zpk, lp2bp_zpk, lp2bs_zpk bilinear

Notes ----- .. versionadded:: 1.1.0

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> fs = 100 >>> bf = 2 * np.pi * np.array(7, 13) >>> filts = signal.lti( *signal.butter(4, bf, btype='bandpass', analog=True, output='zpk')) >>> filtz = signal.lti( *signal.bilinear_zpk(filts.zeros, filts.poles, filts.gain, fs)) >>> wz, hz = signal.freqz_zpk(filtz.zeros, filtz.poles, filtz.gain) >>> ws, hs = signal.freqs_zpk(filts.zeros, filts.poles, filts.gain, worN=fs*wz) >>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hz).clip(1e-15)), label=r'$|H(j \omega)|$') >>> plt.semilogx(wz*fs/(2*np.pi), 20*np.log10(np.abs(hs).clip(1e-15)), label=r'$|H_z(e^j \omega)|$') >>> plt.legend() >>> plt.xlabel('Frequency Hz') >>> plt.ylabel('Magnitude dB') >>> plt.grid()

Return a Blackman window.

The Blackman window is a taper formed by using the first three terms of a summation of cosines. It was designed to have close to the minimal leakage possible. It is close to optimal, only slightly worse than a Kaiser window.

.. warning:: scipy.signal.blackman is deprecated, use scipy.signal.windows.blackman instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Notes ----- The Blackman window is defined as

.. math:: w(n) = 0.42 - 0.5 \cos(2\pi n/M) + 0.08 \cos(4\pi n/M)

The 'exact Blackman' window was designed to null out the third and fourth sidelobes, but has discontinuities at the boundaries, resulting in a 6 dB/oct fall-off. This window is an approximation of the 'exact' window, which does not null the sidelobes as well, but is smooth at the edges, improving the fall-off rate to 18 dB/oct. 3_

Most references to the Blackman window come from the signal processing literature, where it is used as one of many windowing functions for smoothing values. It is also known as an apodization (which means 'removing the foot', i.e. smoothing discontinuities at the beginning and end of the sampled signal) or tapering function. It is known as a 'near optimal' tapering function, almost as good (by some measures) as the Kaiser window.

References ---------- .. 1 Blackman, R.B. and Tukey, J.W., (1958) The measurement of power spectra, Dover Publications, New York. .. 2 Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Upper Saddle River, NJ: Prentice-Hall, 1999, pp. 468-471. .. 3 Harris, Fredric J. (Jan 1978). 'On the use of Windows for Harmonic Analysis with the Discrete Fourier Transform'. Proceedings of the IEEE 66 (1): 51-83. :doi:`10.1109/PROC.1978.10837`.

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.blackman(51) >>> plt.plot(window) >>> plt.title('Blackman window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = np.abs(fftshift(A / abs(A).max())) >>> response = 20 * np.log10(np.maximum(response, 1e-10)) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Blackman window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Return a minimum 4-term Blackman-Harris window.

.. warning:: scipy.signal.blackmanharris is deprecated, use scipy.signal.windows.blackmanharris instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.blackmanharris(51) >>> plt.plot(window) >>> plt.title('Blackman-Harris window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Blackman-Harris window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Calculate Bode magnitude and phase data of a continuous-time system.

Parameters ---------- system : an instance of the LTI class or a tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1 (instance of `lti`) * 2 (num, den) * 3 (zeros, poles, gain) * 4 (A, B, C, D)

w : array_like, optional Array of frequencies (in rad/s). Magnitude and phase data is calculated for every value in this array. If not given a reasonable set will be calculated. n : int, optional Number of frequency points to compute if `w` is not given. The `n` frequencies are logarithmically spaced in an interval chosen to include the influence of the poles and zeros of the system.

Returns ------- w : 1D ndarray Frequency array rad/s mag : 1D ndarray Magnitude array dB phase : 1D ndarray Phase array deg

Notes ----- If (num, den) is passed in for ``system``, coefficients for both the numerator and denominator should be specified in descending exponent order (e.g. ``s^2 + 3s + 5`` would be represented as ``1, 3, 5``).

.. versionadded:: 0.11.0

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> sys = signal.TransferFunction(1, 1, 1) >>> w, mag, phase = signal.bode(sys)

>>> plt.figure() >>> plt.semilogx(w, mag) # Bode magnitude plot >>> plt.figure() >>> plt.semilogx(w, phase) # Bode phase plot >>> plt.show()

Return a Bohman window.

.. warning:: scipy.signal.bohman is deprecated, use scipy.signal.windows.bohman instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.bohman(51) >>> plt.plot(window) >>> plt.title('Bohman window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Bohman window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Return a boxcar or rectangular window.

Also known as a rectangular window or Dirichlet window, this is equivalent to no window at all.

.. warning:: scipy.signal.boxcar is deprecated, use scipy.signal.windows.boxcar instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional Whether the window is symmetric. (Has no effect for boxcar.)

Returns ------- w : ndarray The window, with the maximum value normalized to 1.

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.boxcar(51) >>> plt.plot(window) >>> plt.title('Boxcar window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the boxcar window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

B-spline basis function of order n.

Notes ----- Uses numpy.piecewise and automatic function-generator.

Return (z,p,k) for analog prototype of Nth-order Butterworth filter.

The filter will have an angular (e.g., rad/s) cutoff frequency of 1.

See Also -------- butter : Filter design function using this prototype

Butterworth digital and analog filter design.

Design an Nth-order digital or analog Butterworth filter and return the filter coefficients.

Parameters ---------- N : int The order of the filter. Wn : array_like The critical frequency or frequencies. For lowpass and highpass filters, Wn is a scalar; for bandpass and bandstop filters, Wn is a length-2 sequence.

For a Butterworth filter, this is the point at which the gain drops to 1/sqrt(2) that of the passband (the '-3 dB point').

For digital filters, `Wn` are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`Wn` is thus in half-cycles / sample.)

For analog filters, `Wn` is an angular frequency (e.g. rad/s). btype : 'lowpass', 'highpass', 'bandpass', 'bandstop', optional The type of filter. Default is 'lowpass'. analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. output : 'ba', 'zpk', 'sos', optional Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or second-order sections ('sos'). Default is 'ba' for backwards compatibility, but 'sos' should be used for general-purpose filtering. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- b, a : ndarray, ndarray Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. Only returned if ``output='ba'``. z, p, k : ndarray, ndarray, float Zeros, poles, and system gain of the IIR filter transfer function. Only returned if ``output='zpk'``. sos : ndarray Second-order sections representation of the IIR filter. Only returned if ``output=='sos'``.

See Also -------- buttord, buttap

Notes ----- The Butterworth filter has maximally flat frequency response in the passband.

The ``'sos'`` output parameter was added in 0.16.0.

If the transfer function form ``b, a`` is requested, numerical problems can occur since the conversion between roots and the polynomial coefficients is a numerically sensitive operation, even for N >= 4. It is recommended to work with the SOS representation.

Examples -------- Design an analog filter and plot its frequency response, showing the critical points:

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> b, a = signal.butter(4, 100, 'low', analog=True) >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Butterworth filter frequency response') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.show()

Generate a signal made up of 10 Hz and 20 Hz, sampled at 1 kHz

>>> t = np.linspace(0, 1, 1000, False) # 1 second >>> sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) >>> fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) >>> ax1.plot(t, sig) >>> ax1.set_title('10 Hz and 20 Hz sinusoids') >>> ax1.axis(0, 1, -2, 2)

Design a digital high-pass filter at 15 Hz to remove the 10 Hz tone, and apply it to the signal. (It's recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (``ba``) format):

>>> sos = signal.butter(10, 15, 'hp', fs=1000, output='sos') >>> filtered = signal.sosfilt(sos, sig) >>> ax2.plot(t, filtered) >>> ax2.set_title('After 15 Hz high-pass filter') >>> ax2.axis(0, 1, -2, 2) >>> ax2.set_xlabel('Time seconds') >>> plt.tight_layout() >>> plt.show()

Butterworth filter order selection.

Return the order of the lowest order digital or analog Butterworth filter that loses no more than `gpass` dB in the passband and has at least `gstop` dB attenuation in the stopband.

Parameters ---------- wp, ws : float Passband and stopband edge frequencies.

For digital filters, these are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`wp` and `ws` are thus in half-cycles / sample.) For example:

• Lowpass: wp = 0.2, ws = 0.3
• Highpass: wp = 0.3, ws = 0.2
• Bandpass: wp = 0.2, 0.5, ws = 0.1, 0.6
• Bandstop: wp = 0.1, 0.6, ws = 0.2, 0.5

For analog filters, `wp` and `ws` are angular frequencies (e.g., rad/s). gpass : float The maximum loss in the passband (dB). gstop : float The minimum attenuation in the stopband (dB). analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- ord : int The lowest order for a Butterworth filter which meets specs. wn : ndarray or float The Butterworth natural frequency (i.e. the '3dB frequency'). Should be used with `butter` to give filter results. If `fs` is specified, this is in the same units, and `fs` must also be passed to `butter`.

See Also -------- butter : Filter design using order and critical points cheb1ord : Find order and critical points from passband and stopband spec cheb2ord, ellipord iirfilter : General filter design using order and critical frequencies iirdesign : General filter design using passband and stopband spec

Examples -------- Design an analog bandpass filter with passband within 3 dB from 20 to 50 rad/s, while rejecting at least -40 dB below 14 and above 60 rad/s. Plot its frequency response, showing the passband and stopband constraints in gray.

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> N, Wn = signal.buttord(20, 50, 14, 60, 3, 40, True) >>> b, a = signal.butter(N, Wn, 'band', True) >>> w, h = signal.freqs(b, a, np.logspace(1, 2, 500)) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Butterworth bandpass filter fit to constraints') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.grid(which='both', axis='both') >>> plt.fill(1, 14, 14, 1, -40, -40, 99, 99, '0.9', lw=0) # stop >>> plt.fill(20, 20, 50, 50, -99, -3, -3, -99, '0.9', lw=0) # pass >>> plt.fill(60, 60, 1e9, 1e9, 99, -40, -40, 99, '0.9', lw=0) # stop >>> plt.axis(10, 100, -60, 3) >>> plt.show()

Return (x, phi, psi) at dyadic points ``K/2**J`` from filter coefficients.

Parameters ---------- hk : array_like Coefficients of low-pass filter. J : int, optional Values will be computed at grid points ``K/2**J``. Default is 7.

Returns ------- x : ndarray The dyadic points ``K/2**J`` for ``K=0...N * (2**J)-1`` where ``len(hk) = len(gk) = N+1``. phi : ndarray The scaling function ``phi(x)`` at `x`: ``phi(x) = sum(hk * phi(2x-k))``, where k is from 0 to N. psi : ndarray, optional The wavelet function ``psi(x)`` at `x`: ``phi(x) = sum(gk * phi(2x-k))``, where k is from 0 to N. `psi` is only returned if `gk` is not None.

Notes ----- The algorithm uses the vector cascade algorithm described by Strang and Nguyen in 'Wavelets and Filter Banks'. It builds a dictionary of values and slices for quick reuse. Then inserts vectors into final vector at the end.

Return (z,p,k) for Nth-order Chebyshev type I analog lowpass filter.

The returned filter prototype has `rp` decibels of ripple in the passband.

The filter's angular (e.g. rad/s) cutoff frequency is normalized to 1, defined as the point at which the gain first drops below ``-rp``.

See Also -------- cheby1 : Filter design function using this prototype

Chebyshev type I filter order selection.

Return the order of the lowest order digital or analog Chebyshev Type I filter that loses no more than `gpass` dB in the passband and has at least `gstop` dB attenuation in the stopband.

Parameters ---------- wp, ws : float Passband and stopband edge frequencies.

For digital filters, these are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`wp` and `ws` are thus in half-cycles / sample.) For example:

• Lowpass: wp = 0.2, ws = 0.3
• Highpass: wp = 0.3, ws = 0.2
• Bandpass: wp = 0.2, 0.5, ws = 0.1, 0.6
• Bandstop: wp = 0.1, 0.6, ws = 0.2, 0.5

For analog filters, `wp` and `ws` are angular frequencies (e.g., rad/s). gpass : float The maximum loss in the passband (dB). gstop : float The minimum attenuation in the stopband (dB). analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- ord : int The lowest order for a Chebyshev type I filter that meets specs. wn : ndarray or float The Chebyshev natural frequency (the '3dB frequency') for use with `cheby1` to give filter results. If `fs` is specified, this is in the same units, and `fs` must also be passed to `cheby1`.

See Also -------- cheby1 : Filter design using order and critical points buttord : Find order and critical points from passband and stopband spec cheb2ord, ellipord iirfilter : General filter design using order and critical frequencies iirdesign : General filter design using passband and stopband spec

Examples -------- Design a digital lowpass filter such that the passband is within 3 dB up to 0.2*(fs/2), while rejecting at least -40 dB above 0.3*(fs/2). Plot its frequency response, showing the passband and stopband constraints in gray.

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> N, Wn = signal.cheb1ord(0.2, 0.3, 3, 40) >>> b, a = signal.cheby1(N, 3, Wn, 'low') >>> w, h = signal.freqz(b, a) >>> plt.semilogx(w / np.pi, 20 * np.log10(abs(h))) >>> plt.title('Chebyshev I lowpass filter fit to constraints') >>> plt.xlabel('Normalized frequency') >>> plt.ylabel('Amplitude dB') >>> plt.grid(which='both', axis='both') >>> plt.fill(.01, 0.2, 0.2, .01, -3, -3, -99, -99, '0.9', lw=0) # stop >>> plt.fill(0.3, 0.3, 2, 2, 9, -40, -40, 9, '0.9', lw=0) # pass >>> plt.axis(0.08, 1, -60, 3) >>> plt.show()

Return (z,p,k) for Nth-order Chebyshev type I analog lowpass filter.

The returned filter prototype has `rs` decibels of ripple in the stopband.

The filter's angular (e.g. rad/s) cutoff frequency is normalized to 1, defined as the point at which the gain first reaches ``-rs``.

See Also -------- cheby2 : Filter design function using this prototype

Chebyshev type II filter order selection.

Return the order of the lowest order digital or analog Chebyshev Type II filter that loses no more than `gpass` dB in the passband and has at least `gstop` dB attenuation in the stopband.

Parameters ---------- wp, ws : float Passband and stopband edge frequencies.

For digital filters, these are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`wp` and `ws` are thus in half-cycles / sample.) For example:

• Lowpass: wp = 0.2, ws = 0.3
• Highpass: wp = 0.3, ws = 0.2
• Bandpass: wp = 0.2, 0.5, ws = 0.1, 0.6
• Bandstop: wp = 0.1, 0.6, ws = 0.2, 0.5

For analog filters, `wp` and `ws` are angular frequencies (e.g., rad/s). gpass : float The maximum loss in the passband (dB). gstop : float The minimum attenuation in the stopband (dB). analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- ord : int The lowest order for a Chebyshev type II filter that meets specs. wn : ndarray or float The Chebyshev natural frequency (the '3dB frequency') for use with `cheby2` to give filter results. If `fs` is specified, this is in the same units, and `fs` must also be passed to `cheby2`.

See Also -------- cheby2 : Filter design using order and critical points buttord : Find order and critical points from passband and stopband spec cheb1ord, ellipord iirfilter : General filter design using order and critical frequencies iirdesign : General filter design using passband and stopband spec

Examples -------- Design a digital bandstop filter which rejects -60 dB from 0.2*(fs/2) to 0.5*(fs/2), while staying within 3 dB below 0.1*(fs/2) or above 0.6*(fs/2). Plot its frequency response, showing the passband and stopband constraints in gray.

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> N, Wn = signal.cheb2ord(0.1, 0.6, 0.2, 0.5, 3, 60) >>> b, a = signal.cheby2(N, 60, Wn, 'stop') >>> w, h = signal.freqz(b, a) >>> plt.semilogx(w / np.pi, 20 * np.log10(abs(h))) >>> plt.title('Chebyshev II bandstop filter fit to constraints') >>> plt.xlabel('Normalized frequency') >>> plt.ylabel('Amplitude dB') >>> plt.grid(which='both', axis='both') >>> plt.fill(.01, .1, .1, .01, -3, -3, -99, -99, '0.9', lw=0) # stop >>> plt.fill(.2, .2, .5, .5, 9, -60, -60, 9, '0.9', lw=0) # pass >>> plt.fill(.6, .6, 2, 2, -99, -3, -3, -99, '0.9', lw=0) # stop >>> plt.axis(0.06, 1, -80, 3) >>> plt.show()

Return a Dolph-Chebyshev window.

.. warning:: scipy.signal.chebwin is deprecated, use scipy.signal.windows.chebwin instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. at : float Attenuation (in dB). sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value always normalized to 1

Notes ----- This window optimizes for the narrowest main lobe width for a given order `M` and sidelobe equiripple attenuation `at`, using Chebyshev polynomials. It was originally developed by Dolph to optimize the directionality of radio antenna arrays.

Unlike most windows, the Dolph-Chebyshev is defined in terms of its frequency response:

.. math:: W(k) = \frac \cos\{M \cos^{-1\beta \cos(\frac{\pi k}{M})}

}

\cosh[M \cosh^{-1(\beta)

where

.. math:: \beta = \cosh \left \frac{1}{M} \cosh^{-1}(10^\frac{A}{20}) \right

and 0 <= abs(k) <= M-1. A is the attenuation in decibels (`at`).

The time domain window is then generated using the IFFT, so power-of-two `M` are the fastest to generate, and prime number `M` are the slowest.

The equiripple condition in the frequency domain creates impulses in the time domain, which appear at the ends of the window.

References ---------- .. 1 C. Dolph, 'A current distribution for broadside arrays which optimizes the relationship between beam width and side-lobe level', Proceedings of the IEEE, Vol. 34, Issue 6 .. 2 Peter Lynch, 'The Dolph-Chebyshev Window: A Simple Optimal Filter', American Meteorological Society (April 1997) http://mathsci.ucd.ie/~plynch/Publications/Dolph.pdf .. 3 F. J. Harris, 'On the use of windows for harmonic analysis with the discrete Fourier transforms', Proceedings of the IEEE, Vol. 66, No. 1, January 1978

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.chebwin(51, at=100) >>> plt.plot(window) >>> plt.title('Dolph-Chebyshev window (100 dB)') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the Dolph-Chebyshev window (100 dB)') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Chebyshev type I digital and analog filter design.

Design an Nth-order digital or analog Chebyshev type I filter and return the filter coefficients.

Parameters ---------- N : int The order of the filter. rp : float The maximum ripple allowed below unity gain in the passband. Specified in decibels, as a positive number. Wn : array_like A scalar or length-2 sequence giving the critical frequencies. For Type I filters, this is the point in the transition band at which the gain first drops below -`rp`.

For digital filters, `Wn` are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`Wn` is thus in half-cycles / sample.)

For analog filters, `Wn` is an angular frequency (e.g., rad/s). btype : 'lowpass', 'highpass', 'bandpass', 'bandstop', optional The type of filter. Default is 'lowpass'. analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. output : 'ba', 'zpk', 'sos', optional Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or second-order sections ('sos'). Default is 'ba' for backwards compatibility, but 'sos' should be used for general-purpose filtering. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- b, a : ndarray, ndarray Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. Only returned if ``output='ba'``. z, p, k : ndarray, ndarray, float Zeros, poles, and system gain of the IIR filter transfer function. Only returned if ``output='zpk'``. sos : ndarray Second-order sections representation of the IIR filter. Only returned if ``output=='sos'``.

See Also -------- cheb1ord, cheb1ap

Notes ----- The Chebyshev type I filter maximizes the rate of cutoff between the frequency response's passband and stopband, at the expense of ripple in the passband and increased ringing in the step response.

Type I filters roll off faster than Type II (`cheby2`), but Type II filters do not have any ripple in the passband.

The equiripple passband has N maxima or minima (for example, a 5th-order filter has 3 maxima and 2 minima). Consequently, the DC gain is unity for odd-order filters, or -rp dB for even-order filters.

The ``'sos'`` output parameter was added in 0.16.0.

Examples -------- Design an analog filter and plot its frequency response, showing the critical points:

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> b, a = signal.cheby1(4, 5, 100, 'low', analog=True) >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Chebyshev Type I frequency response (rp=5)') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.axhline(-5, color='green') # rp >>> plt.show()

Generate a signal made up of 10 Hz and 20 Hz, sampled at 1 kHz

>>> t = np.linspace(0, 1, 1000, False) # 1 second >>> sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) >>> fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) >>> ax1.plot(t, sig) >>> ax1.set_title('10 Hz and 20 Hz sinusoids') >>> ax1.axis(0, 1, -2, 2)

Design a digital high-pass filter at 15 Hz to remove the 10 Hz tone, and apply it to the signal. (It's recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (``ba``) format):

>>> sos = signal.cheby1(10, 1, 15, 'hp', fs=1000, output='sos') >>> filtered = signal.sosfilt(sos, sig) >>> ax2.plot(t, filtered) >>> ax2.set_title('After 15 Hz high-pass filter') >>> ax2.axis(0, 1, -2, 2) >>> ax2.set_xlabel('Time seconds') >>> plt.tight_layout() >>> plt.show()

Chebyshev type II digital and analog filter design.

Design an Nth-order digital or analog Chebyshev type II filter and return the filter coefficients.

Parameters ---------- N : int The order of the filter. rs : float The minimum attenuation required in the stop band. Specified in decibels, as a positive number. Wn : array_like A scalar or length-2 sequence giving the critical frequencies. For Type II filters, this is the point in the transition band at which the gain first reaches -`rs`.

For digital filters, `Wn` are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`Wn` is thus in half-cycles / sample.)

For analog filters, `Wn` is an angular frequency (e.g., rad/s). btype : 'lowpass', 'highpass', 'bandpass', 'bandstop', optional The type of filter. Default is 'lowpass'. analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. output : 'ba', 'zpk', 'sos', optional Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or second-order sections ('sos'). Default is 'ba' for backwards compatibility, but 'sos' should be used for general-purpose filtering. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- b, a : ndarray, ndarray Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. Only returned if ``output='ba'``. z, p, k : ndarray, ndarray, float Zeros, poles, and system gain of the IIR filter transfer function. Only returned if ``output='zpk'``. sos : ndarray Second-order sections representation of the IIR filter. Only returned if ``output=='sos'``.

See Also -------- cheb2ord, cheb2ap

Notes ----- The Chebyshev type II filter maximizes the rate of cutoff between the frequency response's passband and stopband, at the expense of ripple in the stopband and increased ringing in the step response.

Type II filters do not roll off as fast as Type I (`cheby1`).

The ``'sos'`` output parameter was added in 0.16.0.

Examples -------- Design an analog filter and plot its frequency response, showing the critical points:

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> b, a = signal.cheby2(4, 40, 100, 'low', analog=True) >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Chebyshev Type II frequency response (rs=40)') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.axhline(-40, color='green') # rs >>> plt.show()

Generate a signal made up of 10 Hz and 20 Hz, sampled at 1 kHz

>>> t = np.linspace(0, 1, 1000, False) # 1 second >>> sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) >>> fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) >>> ax1.plot(t, sig) >>> ax1.set_title('10 Hz and 20 Hz sinusoids') >>> ax1.axis(0, 1, -2, 2)

Design a digital high-pass filter at 17 Hz to remove the 10 Hz tone, and apply it to the signal. (It's recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (``ba``) format):

>>> sos = signal.cheby2(12, 20, 17, 'hp', fs=1000, output='sos') >>> filtered = signal.sosfilt(sos, sig) >>> ax2.plot(t, filtered) >>> ax2.set_title('After 17 Hz high-pass filter') >>> ax2.axis(0, 1, -2, 2) >>> ax2.set_xlabel('Time seconds') >>> plt.show()

Check whether the Constant OverLap Add (COLA) constraint is met

Parameters ---------- window : str or tuple or array_like Desired window to use. If `window` is a string or tuple, it is passed to `get_window` to generate the window values, which are DFT-even by default. See `get_window` for a list of windows and required parameters. If `window` is array_like it will be used directly as the window and its length must be nperseg. nperseg : int Length of each segment. noverlap : int Number of points to overlap between segments. tol : float, optional The allowed variance of a bin's weighted sum from the median bin sum.

Returns ------- verdict : bool `True` if chosen combination satisfies COLA within `tol`, `False` otherwise

See Also -------- check_NOLA: Check whether the Nonzero Overlap Add (NOLA) constraint is met stft: Short Time Fourier Transform istft: Inverse Short Time Fourier Transform

Notes ----- In order to enable inversion of an STFT via the inverse STFT in `istft`, it is sufficient that the signal windowing obeys the constraint of 'Constant OverLap Add' (COLA). This ensures that every point in the input data is equally weighted, thereby avoiding aliasing and allowing full reconstruction.

Some examples of windows that satisfy COLA:

• Rectangular window at overlap of 0, 1/2, 2/3, 3/4, ...
• Bartlett window at overlap of 1/2, 3/4, 5/6, ...
• Hann window at 1/2, 2/3, 3/4, ...
• Any Blackman family window at 2/3 overlap
• Any window with ``noverlap = nperseg-1``

A very comprehensive list of other windows may be found in 2_, wherein the COLA condition is satisfied when the 'Amplitude Flatness' is unity.

.. versionadded:: 0.19.0

References ---------- .. 1 Julius O. Smith III, 'Spectral Audio Signal Processing', W3K Publishing, 2011,ISBN 978-0-9745607-3-1. .. 2 G. Heinzel, A. Ruediger and R. Schilling, 'Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new at-top windows', 2002, http://hdl.handle.net/11858/00-001M-0000-0013-557A-5

Examples -------- >>> from scipy import signal

Confirm COLA condition for rectangular window of 75% (3/4) overlap:

>>> signal.check_COLA(signal.boxcar(100), 100, 75) True

COLA is not true for 25% (1/4) overlap, though:

>>> signal.check_COLA(signal.boxcar(100), 100, 25) False

'Symmetrical' Hann window (for filter design) is not COLA:

>>> signal.check_COLA(signal.hann(120, sym=True), 120, 60) False

'Periodic' or 'DFT-even' Hann window (for FFT analysis) is COLA for overlap of 1/2, 2/3, 3/4, etc.:

>>> signal.check_COLA(signal.hann(120, sym=False), 120, 60) True

>>> signal.check_COLA(signal.hann(120, sym=False), 120, 80) True

>>> signal.check_COLA(signal.hann(120, sym=False), 120, 90) True

Check whether the Nonzero Overlap Add (NOLA) constraint is met

Parameters ---------- window : str or tuple or array_like Desired window to use. If `window` is a string or tuple, it is passed to `get_window` to generate the window values, which are DFT-even by default. See `get_window` for a list of windows and required parameters. If `window` is array_like it will be used directly as the window and its length must be nperseg. nperseg : int Length of each segment. noverlap : int Number of points to overlap between segments. tol : float, optional The allowed variance of a bin's weighted sum from the median bin sum.

Returns ------- verdict : bool `True` if chosen combination satisfies the NOLA constraint within `tol`, `False` otherwise

See Also -------- check_COLA: Check whether the Constant OverLap Add (COLA) constraint is met stft: Short Time Fourier Transform istft: Inverse Short Time Fourier Transform

Notes ----- In order to enable inversion of an STFT via the inverse STFT in `istft`, the signal windowing must obey the constraint of 'nonzero overlap add' (NOLA):

.. math:: \sum_

w^

n-tH \ne 0

for all :math:`n`, where :math:`w` is the window function, :math:`t` is the frame index, and :math:`H` is the hop size (:math:`H` = `nperseg` - `noverlap`).

This ensures that the normalization factors in the denominator of the overlap-add inversion equation are not zero. Only very pathological windows will fail the NOLA constraint.

.. versionadded:: 1.2.0

References ---------- .. 1 Julius O. Smith III, 'Spectral Audio Signal Processing', W3K Publishing, 2011,ISBN 978-0-9745607-3-1. .. 2 G. Heinzel, A. Ruediger and R. Schilling, 'Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new at-top windows', 2002, http://hdl.handle.net/11858/00-001M-0000-0013-557A-5

Examples -------- >>> from scipy import signal

Confirm NOLA condition for rectangular window of 75% (3/4) overlap:

>>> signal.check_NOLA(signal.boxcar(100), 100, 75) True

NOLA is also true for 25% (1/4) overlap:

>>> signal.check_NOLA(signal.boxcar(100), 100, 25) True

'Symmetrical' Hann window (for filter design) is also NOLA:

>>> signal.check_NOLA(signal.hann(120, sym=True), 120, 60) True

As long as there is overlap, it takes quite a pathological window to fail NOLA:

>>> w = np.ones(64, dtype='float') >>> w::2 = 0 >>> signal.check_NOLA(w, 64, 32) False

If there is not enough overlap, a window with zeros at the ends will not work:

>>> signal.check_NOLA(signal.hann(64), 64, 0) False >>> signal.check_NOLA(signal.hann(64), 64, 1) False >>> signal.check_NOLA(signal.hann(64), 64, 2) True

Frequency-swept cosine generator.

In the following, 'Hz' should be interpreted as 'cycles per unit'; there is no requirement here that the unit is one second. The important distinction is that the units of rotation are cycles, not radians. Likewise, `t` could be a measurement of space instead of time.

Parameters ---------- t : array_like Times at which to evaluate the waveform. f0 : float Frequency (e.g. Hz) at time t=0. t1 : float Time at which `f1` is specified. f1 : float Frequency (e.g. Hz) of the waveform at time `t1`. 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 ------- y : ndarray 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.

See Also -------- sweep_poly

Notes ----- 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.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)``

f0 and f1 must be nonzero.

Examples -------- The following will be used in the examples:

>>> from scipy.signal import chirp, spectrogram >>> import matplotlib.pyplot as plt

For the first example, we'll plot the waveform for a linear chirp from 6 Hz to 1 Hz over 10 seconds:

>>> t = np.linspace(0, 10, 1500) >>> w = chirp(t, f0=6, f1=1, t1=10, method='linear') >>> plt.plot(t, w) >>> plt.title('Linear Chirp, f(0)=6, f(10)=1') >>> plt.xlabel('t (sec)') >>> plt.show()

For the remaining examples, we'll use higher frequency ranges, and demonstrate the result using `scipy.signal.spectrogram`. We'll use a 4 second interval sampled at 7200 Hz.

>>> fs = 7200 >>> T = 4 >>> t = np.arange(0, int(T*fs)) / fs

We'll use this function to plot the spectrogram in each example.

>>> def plot_spectrogram(title, w, fs): ... ff, tt, Sxx = spectrogram(w, fs=fs, nperseg=256, nfft=576) ... plt.pcolormesh(tt, ff:145, Sxx:145, cmap='gray_r', shading='gouraud') ... plt.title(title) ... plt.xlabel('t (sec)') ... plt.ylabel('Frequency (Hz)') ... plt.grid() ...

Quadratic chirp from 1500 Hz to 250 Hz (vertex of the parabolic curve of the frequency is at t=0):

>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic') >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f(T)=250', w, fs) >>> plt.show()

Quadratic chirp from 1500 Hz to 250 Hz (vertex of the parabolic curve of the frequency is at t=T):

>>> w = chirp(t, f0=1500, f1=250, t1=T, method='quadratic', ... vertex_zero=False) >>> plot_spectrogram(f'Quadratic Chirp, f(0)=1500, f(T)=250\n' + ... '(vertex_zero=False)', w, fs) >>> plt.show()

Logarithmic chirp from 1500 Hz to 250 Hz:

>>> w = chirp(t, f0=1500, f1=250, t1=T, method='logarithmic') >>> plot_spectrogram(f'Logarithmic Chirp, f(0)=1500, f(T)=250', w, fs) >>> plt.show()

Hyperbolic chirp from 1500 Hz to 250 Hz:

>>> w = chirp(t, f0=1500, f1=250, t1=T, method='hyperbolic') >>> plot_spectrogram(f'Hyperbolic Chirp, f(0)=1500, f(T)=250', w, fs) >>> plt.show()

Find the fastest convolution/correlation method.

This primarily exists to be called during the ``method='auto'`` option in `convolve` and `correlate`. It can also be used to determine the value of ``method`` for many different convolutions of the same dtype/shape. In addition, it supports timing the convolution to adapt the value of ``method`` to a particular set of inputs and/or hardware.

Parameters ---------- in1 : array_like The first argument passed into the convolution function. in2 : array_like The second argument passed into the convolution function. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear convolution of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. measure : bool, optional If True, run and time the convolution of `in1` and `in2` with both methods and return the fastest. If False (default), predict the fastest method using precomputed values.

Returns ------- method : str A string indicating which convolution method is fastest, either 'direct' or 'fft' times : dict, optional A dictionary containing the times (in seconds) needed for each method. This value is only returned if ``measure=True``.

See Also -------- convolve correlate

Notes ----- Generally, this method is 99% accurate for 2D signals and 85% accurate for 1D signals for randomly chosen input sizes. For precision, use ``measure=True`` to find the fastest method by timing the convolution. This can be used to avoid the minimal overhead of finding the fastest ``method`` later, or to adapt the value of ``method`` to a particular set of inputs.

Experiments were run on an Amazon EC2 r5a.2xlarge machine to test this function. These experiments measured the ratio between the time required when using ``method='auto'`` and the time required for the fastest method (i.e., ``ratio = time_auto / min(time_fft, time_direct)``). In these experiments, we found:

* There is a 95% chance of this ratio being less than 1.5 for 1D signals and a 99% chance of being less than 2.5 for 2D signals. * The ratio was always less than 2.5/5 for 1D/2D signals respectively. * This function is most inaccurate for 1D convolutions that take between 1 and 10 milliseconds with ``method='direct'``. A good proxy for this (at least in our experiments) is ``1e6 <= in1.size * in2.size <= 1e7``.

The 2D results almost certainly generalize to 3D/4D/etc because the implementation is the same (the 1D implementation is different).

All the numbers above are specific to the EC2 machine. However, we did find that this function generalizes fairly decently across hardware. The speed tests were of similar quality (and even slightly better) than the same tests performed on the machine to tune this function's numbers (a mid-2014 15-inch MacBook Pro with 16GB RAM and a 2.5GHz Intel i7 processor).

There are cases when `fftconvolve` supports the inputs but this function returns `direct` (e.g., to protect against floating point integer precision).

.. versionadded:: 0.19

Examples -------- Estimate the fastest method for a given input:

>>> from scipy import signal >>> img = np.random.rand(32, 32) >>> filter = np.random.rand(8, 8) >>> method = signal.choose_conv_method(img, filter, mode='same') >>> method 'fft'

This can then be applied to other arrays of the same dtype and shape:

>>> img2 = np.random.rand(32, 32) >>> filter2 = np.random.rand(8, 8) >>> corr2 = signal.correlate(img2, filter2, mode='same', method=method) >>> conv2 = signal.convolve(img2, filter2, mode='same', method=method)

The output of this function (``method``) works with `correlate` and `convolve`.

Sort roots based on magnitude.

Parameters ---------- p : array_like The roots to sort, as a 1-D array.

Returns ------- p_sorted : ndarray Sorted roots. indx : ndarray Array of indices needed to sort the input `p`.

Examples -------- >>> from scipy import signal >>> vals = 1, 4, 1+1.j, 3 >>> p_sorted, indx = signal.cmplx_sort(vals) >>> p_sorted array(1.+0.j, 1.+1.j, 3.+0.j, 4.+0.j) >>> indx array(0, 2, 3, 1)

Estimate the magnitude squared coherence estimate, Cxy, of discrete-time signals X and Y using Welch's method.

``Cxy = abs(Pxy)**2/(Pxx*Pyy)``, where `Pxx` and `Pyy` are power spectral density estimates of X and Y, and `Pxy` is the cross spectral density estimate of X and Y.

Parameters ---------- x : array_like Time series of measurement values y : array_like Time series of measurement values fs : float, optional Sampling frequency of the `x` and `y` time series. Defaults to 1.0. window : str or tuple or array_like, optional Desired window to use. If `window` is a string or tuple, it is passed to `get_window` to generate the window values, which are DFT-even by default. See `get_window` for a list of windows and required parameters. If `window` is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window. nperseg : int, optional Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window. noverlap: int, optional Number of points to overlap between segments. If `None`, ``noverlap = nperseg // 2``. Defaults to `None`. nfft : int, optional Length of the FFT used, if a zero padded FFT is desired. If `None`, the FFT length is `nperseg`. Defaults to `None`. detrend : str or function or `False`, optional Specifies how to detrend each segment. If `detrend` is a string, it is passed as the `type` argument to the `detrend` function. If it is a function, it takes a segment and returns a detrended segment. If `detrend` is `False`, no detrending is done. Defaults to 'constant'. axis : int, optional Axis along which the coherence is computed for both inputs; the default is over the last axis (i.e. ``axis=-1``).

Returns ------- f : ndarray Array of sample frequencies. Cxy : ndarray Magnitude squared coherence of x and y.

See Also -------- periodogram: Simple, optionally modified periodogram lombscargle: Lomb-Scargle periodogram for unevenly sampled data welch: Power spectral density by Welch's method. csd: Cross spectral density by Welch's method.

Notes -------- An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.

.. versionadded:: 0.16.0

References ---------- .. 1 P. Welch, 'The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms', IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967. .. 2 Stoica, Petre, and Randolph Moses, 'Spectral Analysis of Signals' Prentice Hall, 2005

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

Generate two test signals with some common features.

>>> fs = 10e3 >>> N = 1e5 >>> amp = 20 >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = np.arange(N) / fs >>> b, a = signal.butter(2, 0.25, 'low') >>> x = np.random.normal(scale=np.sqrt(noise_power), size=time.shape) >>> y = signal.lfilter(b, a, x) >>> x += amp*np.sin(2*np.pi*freq*time) >>> y += np.random.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)

Compute and plot the coherence.

>>> f, Cxy = signal.coherence(x, y, fs, nperseg=1024) >>> plt.semilogy(f, Cxy) >>> plt.xlabel('frequency Hz') >>> plt.ylabel('Coherence') >>> plt.show()

Transform a continuous to a discrete state-space system.

Parameters ---------- system : a tuple describing the system or an instance of `lti` The following gives the number of elements in the tuple and the interpretation:

* 1: (instance of `lti`) * 2: (num, den) * 3: (zeros, poles, gain) * 4: (A, B, C, D)

dt : float The discretization time step. method : str, optional Which method to use:

* gbt: generalized bilinear transformation * bilinear: Tustin's approximation ('gbt' with alpha=0.5) * euler: Euler (or forward differencing) method ('gbt' with alpha=0) * backward_diff: Backwards differencing ('gbt' with alpha=1.0) * zoh: zero-order hold (default) * foh: first-order hold ( *versionadded: 1.3.0* ) * impulse: equivalent impulse response ( *versionadded: 1.3.0* )

alpha : float within 0, 1, optional The generalized bilinear transformation weighting parameter, which should only be specified with method='gbt', and is ignored otherwise

Returns ------- sysd : tuple containing the discrete system Based on the input type, the output will be of the form

* (num, den, dt) for transfer function input * (zeros, poles, gain, dt) for zeros-poles-gain input * (A, B, C, D, dt) for state-space system input

Notes ----- By default, the routine uses a Zero-Order Hold (zoh) method to perform the transformation. Alternatively, a generalized bilinear transformation may be used, which includes the common Tustin's bilinear approximation, an Euler's method technique, or a backwards differencing technique.

The Zero-Order Hold (zoh) method is based on 1_, the generalized bilinear approximation is based on 2_ and 3_, the First-Order Hold (foh) method is based on 4_.

References ---------- .. 1 https://en.wikipedia.org/wiki/Discretization#Discretization_of_linear_state_space_models

.. 2 http://techteach.no/publications/discretetime_signals_systems/discrete.pdf

.. 3 G. Zhang, X. Chen, and T. Chen, Digital redesign via the generalized bilinear transformation, Int. J. Control, vol. 82, no. 4, pp. 741-754, 2009. (https://www.mypolyuweb.hk/~magzhang/Research/ZCC09_IJC.pdf)

.. 4 G. F. Franklin, J. D. Powell, and M. L. Workman, Digital control of dynamic systems, 3rd ed. Menlo Park, Calif: Addison-Wesley, pp. 204-206, 1998.

Convolve two N-dimensional arrays.

Convolve `in1` and `in2`, with the output size determined by the `mode` argument.

Parameters ---------- in1 : array_like First input. in2 : array_like Second input. Should have the same number of dimensions as `in1`. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear convolution of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. In 'valid' mode, either `in1` or `in2` must be at least as large as the other in every dimension. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. method : str 'auto', 'direct', 'fft', optional A string indicating which method to use to calculate the convolution.

``direct`` The convolution is determined directly from sums, the definition of convolution. ``fft`` The Fourier Transform is used to perform the convolution by calling `fftconvolve`. ``auto`` Automatically chooses direct or Fourier method based on an estimate of which is faster (default). See Notes for more detail.

.. versionadded:: 0.19.0

Returns ------- convolve : array An N-dimensional array containing a subset of the discrete linear convolution of `in1` with `in2`.

See Also -------- numpy.polymul : performs polynomial multiplication (same operation, but also accepts poly1d objects) choose_conv_method : chooses the fastest appropriate convolution method fftconvolve : Always uses the FFT method. oaconvolve : Uses the overlap-add method to do convolution, which is generally faster when the input arrays are large and significantly different in size.

Notes ----- By default, `convolve` and `correlate` use ``method='auto'``, which calls `choose_conv_method` to choose the fastest method using pre-computed values (`choose_conv_method` can also measure real-world timing with a keyword argument). Because `fftconvolve` relies on floating point numbers, there are certain constraints that may force `method=direct` (more detail in `choose_conv_method` docstring).

Examples -------- Smooth a square pulse using a Hann window:

>>> from scipy import signal >>> sig = np.repeat(0., 1., 0., 100) >>> win = signal.hann(50) >>> filtered = signal.convolve(sig, win, mode='same') / sum(win)

>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_win, ax_filt) = plt.subplots(3, 1, sharex=True) >>> ax_orig.plot(sig) >>> ax_orig.set_title('Original pulse') >>> ax_orig.margins(0, 0.1) >>> ax_win.plot(win) >>> ax_win.set_title('Filter impulse response') >>> ax_win.margins(0, 0.1) >>> ax_filt.plot(filtered) >>> ax_filt.set_title('Filtered signal') >>> ax_filt.margins(0, 0.1) >>> fig.tight_layout() >>> fig.show()

Convolve two 2-dimensional arrays.

Convolve `in1` and `in2` with output size determined by `mode`, and boundary conditions determined by `boundary` and `fillvalue`.

Parameters ---------- in1 : array_like First input. in2 : array_like Second input. Should have the same number of dimensions as `in1`. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear convolution of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. In 'valid' mode, either `in1` or `in2` must be at least as large as the other in every dimension. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. boundary : str 'fill', 'wrap', 'symm', optional A flag indicating how to handle boundaries:

``fill`` pad input arrays with fillvalue. (default) ``wrap`` circular boundary conditions. ``symm`` symmetrical boundary conditions.

fillvalue : scalar, optional Value to fill pad input arrays with. Default is 0.

Returns ------- out : ndarray A 2-dimensional array containing a subset of the discrete linear convolution of `in1` with `in2`.

Examples -------- Compute the gradient of an image by 2D convolution with a complex Scharr operator. (Horizontal operator is real, vertical is imaginary.) Use symmetric boundary condition to avoid creating edges at the image boundaries.

>>> from scipy import signal >>> from scipy import misc >>> ascent = misc.ascent() >>> scharr = np.array([ -3-3j, 0-10j, +3 -3j], ... [-10+0j, 0+ 0j, +10 +0j], ... [ -3+3j, 0+10j, +3 +3j]) # Gx + j*Gy >>> grad = signal.convolve2d(ascent, scharr, boundary='symm', mode='same')

>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_mag, ax_ang) = plt.subplots(3, 1, figsize=(6, 15)) >>> ax_orig.imshow(ascent, cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_mag.imshow(np.absolute(grad), cmap='gray') >>> ax_mag.set_title('Gradient magnitude') >>> ax_mag.set_axis_off() >>> ax_ang.imshow(np.angle(grad), cmap='hsv') # hsv is cyclic, like angles >>> ax_ang.set_title('Gradient orientation') >>> ax_ang.set_axis_off() >>> fig.show()

Cross-correlate two N-dimensional arrays.

Cross-correlate `in1` and `in2`, with the output size determined by the `mode` argument.

Parameters ---------- in1 : array_like First input. in2 : array_like Second input. Should have the same number of dimensions as `in1`. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear cross-correlation of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. In 'valid' mode, either `in1` or `in2` must be at least as large as the other in every dimension. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. method : str 'auto', 'direct', 'fft', optional A string indicating which method to use to calculate the correlation.

``direct`` The correlation is determined directly from sums, the definition of correlation. ``fft`` The Fast Fourier Transform is used to perform the correlation more quickly (only available for numerical arrays.) ``auto`` Automatically chooses direct or Fourier method based on an estimate of which is faster (default). See `convolve` Notes for more detail.

.. versionadded:: 0.19.0

Returns ------- correlate : array An N-dimensional array containing a subset of the discrete linear cross-correlation of `in1` with `in2`.

See Also -------- choose_conv_method : contains more documentation on `method`.

Notes ----- The correlation z of two d-dimensional arrays x and y is defined as::

z...,k,... = sum..., i_l, ... x..., i_l,... * conj(y..., i_l - k,...)

This way, if x and y are 1-D arrays and ``z = correlate(x, y, 'full')`` then

.. math::

zk = (x * y)(k - N + 1) = \sum_l=0^ ||x||-1x_l y_l-k+N-1^*

for :math:`k = 0, 1, ..., ||x|| + ||y|| - 2`

where :math:`||x||` is the length of ``x``, :math:`N = \max(||x||,||y||)`, and :math:`y_m` is 0 when m is outside the range of y.

``method='fft'`` only works for numerical arrays as it relies on `fftconvolve`. In certain cases (i.e., arrays of objects or when rounding integers can lose precision), ``method='direct'`` is always used.

When using 'same' mode with even-length inputs, the outputs of `correlate` and `correlate2d` differ: There is a 1-index offset between them.

Examples -------- Implement a matched filter using cross-correlation, to recover a signal that has passed through a noisy channel.

>>> from scipy import signal >>> sig = np.repeat(0., 1., 1., 0., 1., 0., 0., 1., 128) >>> sig_noise = sig + np.random.randn(len(sig)) >>> corr = signal.correlate(sig_noise, np.ones(128), mode='same') / 128

>>> import matplotlib.pyplot as plt >>> clock = np.arange(64, len(sig), 128) >>> fig, (ax_orig, ax_noise, ax_corr) = plt.subplots(3, 1, sharex=True) >>> ax_orig.plot(sig) >>> ax_orig.plot(clock, sigclock, 'ro') >>> ax_orig.set_title('Original signal') >>> ax_noise.plot(sig_noise) >>> ax_noise.set_title('Signal with noise') >>> ax_corr.plot(corr) >>> ax_corr.plot(clock, corrclock, 'ro') >>> ax_corr.axhline(0.5, ls=':') >>> ax_corr.set_title('Cross-correlated with rectangular pulse') >>> ax_orig.margins(0, 0.1) >>> fig.tight_layout() >>> fig.show()

Cross-correlate two 2-dimensional arrays.

Cross correlate `in1` and `in2` with output size determined by `mode`, and boundary conditions determined by `boundary` and `fillvalue`.

Parameters ---------- in1 : array_like First input. in2 : array_like Second input. Should have the same number of dimensions as `in1`. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear cross-correlation of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. In 'valid' mode, either `in1` or `in2` must be at least as large as the other in every dimension. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. boundary : str 'fill', 'wrap', 'symm', optional A flag indicating how to handle boundaries:

``fill`` pad input arrays with fillvalue. (default) ``wrap`` circular boundary conditions. ``symm`` symmetrical boundary conditions.

fillvalue : scalar, optional Value to fill pad input arrays with. Default is 0.

Returns ------- correlate2d : ndarray A 2-dimensional array containing a subset of the discrete linear cross-correlation of `in1` with `in2`.

Notes ----- When using 'same' mode with even-length inputs, the outputs of `correlate` and `correlate2d` differ: There is a 1-index offset between them.

Examples -------- Use 2D cross-correlation to find the location of a template in a noisy image:

>>> from scipy import signal >>> from scipy import misc >>> face = misc.face(gray=True) - misc.face(gray=True).mean() >>> template = np.copy(face300:365, 670:750) # right eye >>> template -= template.mean() >>> face = face + np.random.randn( *face.shape) * 50 # add noise >>> corr = signal.correlate2d(face, template, boundary='symm', mode='same') >>> y, x = np.unravel_index(np.argmax(corr), corr.shape) # find the match

>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_template, ax_corr) = plt.subplots(3, 1, ... figsize=(6, 15)) >>> ax_orig.imshow(face, cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_template.imshow(template, cmap='gray') >>> ax_template.set_title('Template') >>> ax_template.set_axis_off() >>> ax_corr.imshow(corr, cmap='gray') >>> ax_corr.set_title('Cross-correlation') >>> ax_corr.set_axis_off() >>> ax_orig.plot(x, y, 'ro') >>> fig.show()

Return a window with a simple cosine shape.

.. warning:: scipy.signal.cosine is deprecated, use scipy.signal.windows.cosine instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Notes -----

.. versionadded:: 0.13.0

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.cosine(51) >>> plt.plot(window) >>> plt.title('Cosine window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the cosine window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample') >>> plt.show()

Estimate the cross power spectral density, Pxy, using Welch's method.

Parameters ---------- x : array_like Time series of measurement values y : array_like Time series of measurement values fs : float, optional Sampling frequency of the `x` and `y` time series. Defaults to 1.0. window : str or tuple or array_like, optional Desired window to use. If `window` is a string or tuple, it is passed to `get_window` to generate the window values, which are DFT-even by default. See `get_window` for a list of windows and required parameters. If `window` is array_like it will be used directly as the window and its length must be nperseg. Defaults to a Hann window. nperseg : int, optional Length of each segment. Defaults to None, but if window is str or tuple, is set to 256, and if window is array_like, is set to the length of the window. noverlap: int, optional Number of points to overlap between segments. If `None`, ``noverlap = nperseg // 2``. Defaults to `None`. nfft : int, optional Length of the FFT used, if a zero padded FFT is desired. If `None`, the FFT length is `nperseg`. Defaults to `None`. detrend : str or function or `False`, optional Specifies how to detrend each segment. If `detrend` is a string, it is passed as the `type` argument to the `detrend` function. If it is a function, it takes a segment and returns a detrended segment. If `detrend` is `False`, no detrending is done. Defaults to 'constant'. return_onesided : bool, optional If `True`, return a one-sided spectrum for real data. If `False` return a two-sided spectrum. Defaults to `True`, but for complex data, a two-sided spectrum is always returned. scaling : 'density', 'spectrum' , optional Selects between computing the cross spectral density ('density') where `Pxy` has units of V**2/Hz and computing the cross spectrum ('spectrum') where `Pxy` has units of V**2, if `x` and `y` are measured in V and `fs` is measured in Hz. Defaults to 'density' axis : int, optional Axis along which the CSD is computed for both inputs; the default is over the last axis (i.e. ``axis=-1``). average : 'mean', 'median' , optional Method to use when averaging periodograms. Defaults to 'mean'.

.. versionadded:: 1.2.0

Returns ------- f : ndarray Array of sample frequencies. Pxy : ndarray Cross spectral density or cross power spectrum of x,y.

See Also -------- periodogram: Simple, optionally modified periodogram lombscargle: Lomb-Scargle periodogram for unevenly sampled data welch: Power spectral density by Welch's method. Equivalent to csd(x,x) coherence: Magnitude squared coherence by Welch's method.

Notes -------- By convention, Pxy is computed with the conjugate FFT of X multiplied by the FFT of Y.

If the input series differ in length, the shorter series will be zero-padded to match.

An appropriate amount of overlap will depend on the choice of window and on your requirements. For the default Hann window an overlap of 50% is a reasonable trade off between accurately estimating the signal power, while not over counting any of the data. Narrower windows may require a larger overlap.

.. versionadded:: 0.16.0

References ---------- .. 1 P. Welch, 'The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms', IEEE Trans. Audio Electroacoust. vol. 15, pp. 70-73, 1967. .. 2 Rabiner, Lawrence R., and B. Gold. 'Theory and Application of Digital Signal Processing' Prentice-Hall, pp. 414-419, 1975

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

Generate two test signals with some common features.

>>> fs = 10e3 >>> N = 1e5 >>> amp = 20 >>> freq = 1234.0 >>> noise_power = 0.001 * fs / 2 >>> time = np.arange(N) / fs >>> b, a = signal.butter(2, 0.25, 'low') >>> x = np.random.normal(scale=np.sqrt(noise_power), size=time.shape) >>> y = signal.lfilter(b, a, x) >>> x += amp*np.sin(2*np.pi*freq*time) >>> y += np.random.normal(scale=0.1*np.sqrt(noise_power), size=time.shape)

Compute and plot the magnitude of the cross spectral density.

>>> f, Pxy = signal.csd(x, y, fs, nperseg=1024) >>> plt.semilogy(f, np.abs(Pxy)) >>> plt.xlabel('frequency Hz') >>> plt.ylabel('CSD V**2/Hz') >>> plt.show()

Compute cubic spline coefficients for rank-1 array.

Find the cubic spline coefficients for a 1-D signal assuming mirror-symmetric boundary conditions. To obtain the signal back from the spline representation mirror-symmetric-convolve these coefficients with a length 3 FIR window 1.0, 4.0, 1.0/ 6.0 .

Parameters ---------- signal : ndarray A rank-1 array representing samples of a signal. lamb : float, optional Smoothing coefficient, default is 0.0.

Returns ------- c : ndarray Cubic spline coefficients.

Evaluate a spline at the new set of points.

`dx` is the old sample-spacing while `x0` was the old origin. In other-words the old-sample points (knot-points) for which the `cj` represent spline coefficients were at equally-spaced points of:

oldx = x0 + j*dx j=0...N-1, with N=len(cj)

Edges are handled using mirror-symmetric boundary conditions.

A cubic B-spline.

This is a special case of `bspline`, and equivalent to ``bspline(x, 3)``.

Continuous wavelet transform.

Performs a continuous wavelet transform on `data`, using the `wavelet` function. A CWT performs a convolution with `data` using the `wavelet` function, which is characterized by a width parameter and length parameter. The `wavelet` function is allowed to be complex.

Parameters ---------- data : (N,) ndarray data on which to perform the transform. wavelet : function Wavelet function, which should take 2 arguments. The first argument is the number of points that the returned vector will have (len(wavelet(length,width)) == length). The second is a width parameter, defining the size of the wavelet (e.g. standard deviation of a gaussian). See `ricker`, which satisfies these requirements. widths : (M,) sequence Widths to use for transform. dtype : data-type, optional The desired data type of output. Defaults to ``float64`` if the output of `wavelet` is real and ``complex128`` if it is complex.

.. versionadded:: 1.4.0

kwargs Keyword arguments passed to wavelet function.

.. versionadded:: 1.4.0

Returns ------- cwt: (M, N) ndarray Will have shape of (len(widths), len(data)).

Notes -----

.. versionadded:: 1.4.0

For non-symmetric, complex-valued wavelets, the input signal is convolved with the time-reversed complex-conjugate of the wavelet data 1.

::

length = min(10 * widthii, len(data)) cwtii,: = signal.convolve(data, np.conj(wavelet(length, widthii, **kwargs))::-1, mode='same')

References ---------- .. 1 S. Mallat, 'A Wavelet Tour of Signal Processing (3rd Edition)', Academic Press, 2009.

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt >>> t = np.linspace(-1, 1, 200, endpoint=False) >>> sig = np.cos(2 * np.pi * 7 * t) + signal.gausspulse(t - 0.4, fc=2) >>> widths = np.arange(1, 31) >>> cwtmatr = signal.cwt(sig, signal.ricker, widths) >>> plt.imshow(cwtmatr, extent=-1, 1, 1, 31, cmap='PRGn', aspect='auto', ... vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max()) >>> plt.show()

The coefficients for the FIR low-pass filter producing Daubechies wavelets.

p>=1 gives the order of the zero at f=1/2. There are 2p filter coefficients.

Parameters ---------- p : int Order of the zero at f=1/2, can have values from 1 to 34.

Returns ------- daub : ndarray Return

Calculate Bode magnitude and phase data of a discrete-time system.

Parameters ---------- system : an instance of the LTI class or a tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1 (instance of `dlti`) * 2 (num, den, dt) * 3 (zeros, poles, gain, dt) * 4 (A, B, C, D, dt)

w : array_like, optional Array of frequencies (in radians/sample). Magnitude and phase data is calculated for every value in this array. If not given a reasonable set will be calculated. n : int, optional Number of frequency points to compute if `w` is not given. The `n` frequencies are logarithmically spaced in an interval chosen to include the influence of the poles and zeros of the system.

Returns ------- w : 1D ndarray Frequency array rad/time_unit mag : 1D ndarray Magnitude array dB phase : 1D ndarray Phase array deg

Notes ----- If (num, den) is passed in for ``system``, coefficients for both the numerator and denominator should be specified in descending exponent order (e.g. ``z^2 + 3z + 5`` would be represented as ``1, 3, 5``).

.. versionadded:: 0.18.0

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

Transfer function: H(z) = 1 / (z^2 + 2z + 3)

>>> sys = signal.TransferFunction(1, 1, 2, 3, dt=0.05)

Equivalent: sys.bode()

>>> w, mag, phase = signal.dbode(sys)

>>> plt.figure() >>> plt.semilogx(w, mag) # Bode magnitude plot >>> plt.figure() >>> plt.semilogx(w, phase) # Bode phase plot >>> plt.show()

Downsample the signal after applying an anti-aliasing filter.

By default, an order 8 Chebyshev type I filter is used. A 30 point FIR filter with Hamming window is used if `ftype` is 'fir'.

Parameters ---------- x : array_like The signal to be downsampled, as an N-dimensional array. q : int The downsampling factor. When using IIR downsampling, it is recommended to call `decimate` multiple times for downsampling factors higher than 13. n : int, optional The order of the filter (1 less than the length for 'fir'). Defaults to 8 for 'iir' and 20 times the downsampling factor for 'fir'. ftype : str 'iir', 'fir' or ``dlti`` instance, optional If 'iir' or 'fir', specifies the type of lowpass filter. If an instance of an `dlti` object, uses that object to filter before downsampling. axis : int, optional The axis along which to decimate. zero_phase : bool, optional Prevent phase shift by filtering with `filtfilt` instead of `lfilter` when using an IIR filter, and shifting the outputs back by the filter's group delay when using an FIR filter. The default value of ``True`` is recommended, since a phase shift is generally not desired.

.. versionadded:: 0.18.0

Returns ------- y : ndarray The down-sampled signal.

See Also -------- resample : Resample up or down using the FFT method. resample_poly : Resample using polyphase filtering and an FIR filter.

Notes ----- The ``zero_phase`` keyword was added in 0.18.0. The possibility to use instances of ``dlti`` as ``ftype`` was added in 0.18.0.

Deconvolves ``divisor`` out of ``signal`` using inverse filtering.

Returns the quotient and remainder such that ``signal = convolve(divisor, quotient) + remainder``

Parameters ---------- signal : array_like Signal data, typically a recorded signal divisor : array_like Divisor data, typically an impulse response or filter that was applied to the original signal

Returns ------- quotient : ndarray Quotient, typically the recovered original signal remainder : ndarray Remainder

Examples -------- Deconvolve a signal that's been filtered:

>>> from scipy import signal >>> original = 0, 1, 0, 0, 1, 1, 0, 0 >>> impulse_response = 2, 1 >>> recorded = signal.convolve(impulse_response, original) >>> recorded array(0, 2, 1, 0, 2, 3, 1, 0, 0) >>> recovered, remainder = signal.deconvolve(recorded, impulse_response) >>> recovered array( 0., 1., 0., 0., 1., 1., 0., 0.)

See Also -------- numpy.polydiv : performs polynomial division (same operation, but also accepts poly1d objects)

Remove linear trend along axis from data.

Parameters ---------- data : array_like The input data. axis : int, optional The axis along which to detrend the data. By default this is the last axis (-1). type : 'linear', 'constant', optional The type of detrending. If ``type == 'linear'`` (default), the result of a linear least-squares fit to `data` is subtracted from `data`. If ``type == 'constant'``, only the mean of `data` is subtracted. bp : array_like of ints, optional A sequence of break points. If given, an individual linear fit is performed for each part of `data` between two break points. Break points are specified as indices into `data`. This parameter only has an effect when ``type == 'linear'``. overwrite_data : bool, optional If True, perform in place detrending and avoid a copy. Default is False

Returns ------- ret : ndarray The detrended input data.

Examples -------- >>> from scipy import signal >>> randgen = np.random.RandomState(9) >>> npoints = 1000 >>> noise = randgen.randn(npoints) >>> x = 3 + 2*np.linspace(0, 1, npoints) + noise >>> (signal.detrend(x) - noise).max() < 0.01 True

Calculate the frequency response of a discrete-time system.

Parameters ---------- system : an instance of the `dlti` class or a tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1 (instance of `dlti`) * 2 (numerator, denominator, dt) * 3 (zeros, poles, gain, dt) * 4 (A, B, C, D, dt)

w : array_like, optional Array of frequencies (in radians/sample). Magnitude and phase data is calculated for every value in this array. If not given a reasonable set will be calculated. n : int, optional Number of frequency points to compute if `w` is not given. The `n` frequencies are logarithmically spaced in an interval chosen to include the influence of the poles and zeros of the system. whole : bool, optional Normally, if 'w' is not given, frequencies are computed from 0 to the Nyquist frequency, pi radians/sample (upper-half of unit-circle). If `whole` is True, compute frequencies from 0 to 2*pi radians/sample.

Returns ------- w : 1D ndarray Frequency array radians/sample H : 1D ndarray Array of complex magnitude values

Notes ----- If (num, den) is passed in for ``system``, coefficients for both the numerator and denominator should be specified in descending exponent order (e.g. ``z^2 + 3z + 5`` would be represented as ``1, 3, 5``).

.. versionadded:: 0.18.0

Examples -------- Generating the Nyquist plot of a transfer function

>>> from scipy import signal >>> import matplotlib.pyplot as plt

Transfer function: H(z) = 1 / (z^2 + 2z + 3)

>>> sys = signal.TransferFunction(1, 1, 2, 3, dt=0.05)

>>> w, H = signal.dfreqresp(sys)

>>> plt.figure() >>> plt.plot(H.real, H.imag, 'b') >>> plt.plot(H.real, -H.imag, 'r') >>> plt.show()

Impulse response of discrete-time system.

Parameters ---------- system : tuple of array_like or instance of `dlti` A tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1: (instance of `dlti`) * 3: (num, den, dt) * 4: (zeros, poles, gain, dt) * 5: (A, B, C, D, dt)

x0 : array_like, optional Initial state-vector. Defaults to zero. t : array_like, optional Time points. Computed if not given. n : int, optional The number of time points to compute (if `t` is not given).

Returns ------- tout : ndarray Time values for the output, as a 1-D array. yout : tuple of ndarray Impulse response of system. Each element of the tuple represents the output of the system based on an impulse in each input.

See Also -------- impulse, dstep, dlsim, cont2discrete

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> butter = signal.dlti( *signal.butter(3, 0.5)) >>> t, y = signal.dimpulse(butter, n=25) >>> plt.step(t, np.squeeze(y)) >>> plt.grid() >>> plt.xlabel('n samples') >>> plt.ylabel('Amplitude')

Simulate output of a discrete-time linear system.

Parameters ---------- system : tuple of array_like or instance of `dlti` A tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1: (instance of `dlti`) * 3: (num, den, dt) * 4: (zeros, poles, gain, dt) * 5: (A, B, C, D, dt)

u : array_like An input array describing the input at each time `t` (interpolation is assumed between given times). If there are multiple inputs, then each column of the rank-2 array represents an input. t : array_like, optional The time steps at which the input is defined. If `t` is given, it must be the same length as `u`, and the final value in `t` determines the number of steps returned in the output. x0 : array_like, optional The initial conditions on the state vector (zero by default).

Returns ------- tout : ndarray Time values for the output, as a 1-D array. yout : ndarray System response, as a 1-D array. xout : ndarray, optional Time-evolution of the state-vector. Only generated if the input is a `StateSpace` system.

See Also -------- lsim, dstep, dimpulse, cont2discrete

Examples -------- A simple integrator transfer function with a discrete time step of 1.0 could be implemented as:

>>> from scipy import signal >>> tf = (1.0,, 1.0, -1.0, 1.0) >>> t_in = 0.0, 1.0, 2.0, 3.0 >>> u = np.asarray(0.0, 0.0, 1.0, 1.0) >>> t_out, y = signal.dlsim(tf, u, t=t_in) >>> y.T array([ 0., 0., 0., 1.])

Step response of discrete-time system.

Parameters ---------- system : tuple of array_like A tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1: (instance of `dlti`) * 3: (num, den, dt) * 4: (zeros, poles, gain, dt) * 5: (A, B, C, D, dt)

x0 : array_like, optional Initial state-vector. Defaults to zero. t : array_like, optional Time points. Computed if not given. n : int, optional The number of time points to compute (if `t` is not given).

Returns ------- tout : ndarray Output time points, as a 1-D array. yout : tuple of ndarray Step response of system. Each element of the tuple represents the output of the system based on a step response to each input.

See Also -------- step, dimpulse, dlsim, cont2discrete

Examples -------- >>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> butter = signal.dlti( *signal.butter(3, 0.5)) >>> t, y = signal.dstep(butter, n=25) >>> plt.step(t, np.squeeze(y)) >>> plt.grid() >>> plt.xlabel('n samples') >>> plt.ylabel('Amplitude')

Elliptic (Cauer) digital and analog filter design.

Design an Nth-order digital or analog elliptic filter and return the filter coefficients.

Parameters ---------- N : int The order of the filter. rp : float The maximum ripple allowed below unity gain in the passband. Specified in decibels, as a positive number. rs : float The minimum attenuation required in the stop band. Specified in decibels, as a positive number. Wn : array_like A scalar or length-2 sequence giving the critical frequencies. For elliptic filters, this is the point in the transition band at which the gain first drops below -`rp`.

For digital filters, `Wn` are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`Wn` is thus in half-cycles / sample.)

For analog filters, `Wn` is an angular frequency (e.g., rad/s). btype : 'lowpass', 'highpass', 'bandpass', 'bandstop', optional The type of filter. Default is 'lowpass'. analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. output : 'ba', 'zpk', 'sos', optional Type of output: numerator/denominator ('ba'), pole-zero ('zpk'), or second-order sections ('sos'). Default is 'ba' for backwards compatibility, but 'sos' should be used for general-purpose filtering. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- b, a : ndarray, ndarray Numerator (`b`) and denominator (`a`) polynomials of the IIR filter. Only returned if ``output='ba'``. z, p, k : ndarray, ndarray, float Zeros, poles, and system gain of the IIR filter transfer function. Only returned if ``output='zpk'``. sos : ndarray Second-order sections representation of the IIR filter. Only returned if ``output=='sos'``.

See Also -------- ellipord, ellipap

Notes ----- Also known as Cauer or Zolotarev filters, the elliptical filter maximizes the rate of transition between the frequency response's passband and stopband, at the expense of ripple in both, and increased ringing in the step response.

As `rp` approaches 0, the elliptical filter becomes a Chebyshev type II filter (`cheby2`). As `rs` approaches 0, it becomes a Chebyshev type I filter (`cheby1`). As both approach 0, it becomes a Butterworth filter (`butter`).

The equiripple passband has N maxima or minima (for example, a 5th-order filter has 3 maxima and 2 minima). Consequently, the DC gain is unity for odd-order filters, or -rp dB for even-order filters.

The ``'sos'`` output parameter was added in 0.16.0.

Examples -------- Design an analog filter and plot its frequency response, showing the critical points:

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> b, a = signal.ellip(4, 5, 40, 100, 'low', analog=True) >>> w, h = signal.freqs(b, a) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Elliptic filter frequency response (rp=5, rs=40)') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.margins(0, 0.1) >>> plt.grid(which='both', axis='both') >>> plt.axvline(100, color='green') # cutoff frequency >>> plt.axhline(-40, color='green') # rs >>> plt.axhline(-5, color='green') # rp >>> plt.show()

Generate a signal made up of 10 Hz and 20 Hz, sampled at 1 kHz

>>> t = np.linspace(0, 1, 1000, False) # 1 second >>> sig = np.sin(2*np.pi*10*t) + np.sin(2*np.pi*20*t) >>> fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) >>> ax1.plot(t, sig) >>> ax1.set_title('10 Hz and 20 Hz sinusoids') >>> ax1.axis(0, 1, -2, 2)

Design a digital high-pass filter at 17 Hz to remove the 10 Hz tone, and apply it to the signal. (It's recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (``ba``) format):

>>> sos = signal.ellip(8, 1, 100, 17, 'hp', fs=1000, output='sos') >>> filtered = signal.sosfilt(sos, sig) >>> ax2.plot(t, filtered) >>> ax2.set_title('After 17 Hz high-pass filter') >>> ax2.axis(0, 1, -2, 2) >>> ax2.set_xlabel('Time seconds') >>> plt.tight_layout() >>> plt.show()

Return (z,p,k) of Nth-order elliptic analog lowpass filter.

The filter is a normalized prototype that has `rp` decibels of ripple in the passband and a stopband `rs` decibels down.

The filter's angular (e.g., rad/s) cutoff frequency is normalized to 1, defined as the point at which the gain first drops below ``-rp``.

See Also -------- ellip : Filter design function using this prototype

References ---------- .. 1 Lutova, Tosic, and Evans, 'Filter Design for Signal Processing', Chapters 5 and 12.

Elliptic (Cauer) filter order selection.

Return the order of the lowest order digital or analog elliptic filter that loses no more than `gpass` dB in the passband and has at least `gstop` dB attenuation in the stopband.

Parameters ---------- wp, ws : float Passband and stopband edge frequencies.

For digital filters, these are in the same units as `fs`. By default, `fs` is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (`wp` and `ws` are thus in half-cycles / sample.) For example:

• Lowpass: wp = 0.2, ws = 0.3
• Highpass: wp = 0.3, ws = 0.2
• Bandpass: wp = 0.2, 0.5, ws = 0.1, 0.6
• Bandstop: wp = 0.1, 0.6, ws = 0.2, 0.5

For analog filters, `wp` and `ws` are angular frequencies (e.g., rad/s). gpass : float The maximum loss in the passband (dB). gstop : float The minimum attenuation in the stopband (dB). analog : bool, optional When True, return an analog filter, otherwise a digital filter is returned. fs : float, optional The sampling frequency of the digital system.

.. versionadded:: 1.2.0

Returns ------- ord : int The lowest order for an Elliptic (Cauer) filter that meets specs. wn : ndarray or float The Chebyshev natural frequency (the '3dB frequency') for use with `ellip` to give filter results. If `fs` is specified, this is in the same units, and `fs` must also be passed to `ellip`.

See Also -------- ellip : Filter design using order and critical points buttord : Find order and critical points from passband and stopband spec cheb1ord, cheb2ord iirfilter : General filter design using order and critical frequencies iirdesign : General filter design using passband and stopband spec

Examples -------- Design an analog highpass filter such that the passband is within 3 dB above 30 rad/s, while rejecting -60 dB at 10 rad/s. Plot its frequency response, showing the passband and stopband constraints in gray.

>>> from scipy import signal >>> import matplotlib.pyplot as plt

>>> N, Wn = signal.ellipord(30, 10, 3, 60, True) >>> b, a = signal.ellip(N, 3, 60, Wn, 'high', True) >>> w, h = signal.freqs(b, a, np.logspace(0, 3, 500)) >>> plt.semilogx(w, 20 * np.log10(abs(h))) >>> plt.title('Elliptical highpass filter fit to constraints') >>> plt.xlabel('Frequency radians / second') >>> plt.ylabel('Amplitude dB') >>> plt.grid(which='both', axis='both') >>> plt.fill(.1, 10, 10, .1, 1e4, 1e4, -60, -60, '0.9', lw=0) # stop >>> plt.fill(30, 30, 1e9, 1e9, -99, -3, -3, -99, '0.9', lw=0) # pass >>> plt.axis(1, 300, -80, 3) >>> plt.show()

Return an exponential (or Poisson) window.

.. warning:: scipy.signal.exponential is deprecated, use scipy.signal.windows.exponential instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. center : float, optional Parameter defining the center location of the window function. The default value if not given is ``center = (M-1) / 2``. This parameter must take its default value for symmetric windows. tau : float, optional Parameter defining the decay. For ``center = 0`` use ``tau = -(M-1) / ln(x)`` if ``x`` is the fraction of the window remaining at the end. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Notes ----- The Exponential window is defined as

.. math:: w(n) = e^

|n-center| / \tau

}

References ---------- S. Gade and H. Herlufsen, 'Windows to FFT analysis (Part I)', Technical Review 3, Bruel & Kjaer, 1987.

Examples -------- Plot the symmetric window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> M = 51 >>> tau = 3.0 >>> window = signal.exponential(M, tau=tau) >>> plt.plot(window) >>> plt.title('Exponential Window (tau=3.0)') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -35, 0) >>> plt.title('Frequency response of the Exponential window (tau=3.0)') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

This function can also generate non-symmetric windows:

>>> tau2 = -(M-1) / np.log(0.01) >>> window2 = signal.exponential(M, 0, tau2, False) >>> plt.figure() >>> plt.plot(window2) >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

Convolve two N-dimensional arrays using FFT.

Convolve `in1` and `in2` using the fast Fourier transform method, with the output size determined by the `mode` argument.

This is generally much faster than `convolve` for large arrays (n > ~500), but can be slower when only a few output values are needed, and can only output float arrays (int or object array inputs will be cast to float).

As of v0.19, `convolve` automatically chooses this method or the direct method based on an estimation of which is faster.

Parameters ---------- in1 : array_like First input. in2 : array_like Second input. Should have the same number of dimensions as `in1`. mode : str 'full', 'valid', 'same', optional A string indicating the size of the output:

``full`` The output is the full discrete linear convolution of the inputs. (Default) ``valid`` The output consists only of those elements that do not rely on the zero-padding. In 'valid' mode, either `in1` or `in2` must be at least as large as the other in every dimension. ``same`` The output is the same size as `in1`, centered with respect to the 'full' output. axes : int or array_like of ints or None, optional Axes over which to compute the convolution. The default is over all axes.

Returns ------- out : array An N-dimensional array containing a subset of the discrete linear convolution of `in1` with `in2`.

See Also -------- convolve : Uses the direct convolution or FFT convolution algorithm depending on which is faster. oaconvolve : Uses the overlap-add method to do convolution, which is generally faster when the input arrays are large and significantly different in size.

Examples -------- Autocorrelation of white noise is an impulse.

>>> from scipy import signal >>> sig = np.random.randn(1000) >>> autocorr = signal.fftconvolve(sig, sig::-1, mode='full')

>>> import matplotlib.pyplot as plt >>> fig, (ax_orig, ax_mag) = plt.subplots(2, 1) >>> ax_orig.plot(sig) >>> ax_orig.set_title('White noise') >>> ax_mag.plot(np.arange(-len(sig)+1,len(sig)), autocorr) >>> ax_mag.set_title('Autocorrelation') >>> fig.tight_layout() >>> fig.show()

Gaussian blur implemented using FFT convolution. Notice the dark borders around the image, due to the zero-padding beyond its boundaries. The `convolve2d` function allows for other types of image boundaries, but is far slower.

>>> from scipy import misc >>> face = misc.face(gray=True) >>> kernel = np.outer(signal.gaussian(70, 8), signal.gaussian(70, 8)) >>> blurred = signal.fftconvolve(face, kernel, mode='same')

>>> fig, (ax_orig, ax_kernel, ax_blurred) = plt.subplots(3, 1, ... figsize=(6, 15)) >>> ax_orig.imshow(face, cmap='gray') >>> ax_orig.set_title('Original') >>> ax_orig.set_axis_off() >>> ax_kernel.imshow(kernel, cmap='gray') >>> ax_kernel.set_title('Gaussian kernel') >>> ax_kernel.set_axis_off() >>> ax_blurred.imshow(blurred, cmap='gray') >>> ax_blurred.set_title('Blurred') >>> ax_blurred.set_axis_off() >>> fig.show()

Apply a digital filter forward and backward to a signal.

This function applies a linear digital filter twice, once forward and once backwards. The combined filter has zero phase and a filter order twice that of the original.

The function provides options for handling the edges of the signal.

The function `sosfiltfilt` (and filter design using ``output='sos'``) should be preferred over `filtfilt` for most filtering tasks, as second-order sections have fewer numerical problems.

Parameters ---------- b : (N,) array_like The numerator coefficient vector of the filter. a : (N,) array_like The denominator coefficient vector of the filter. If ``a0`` is not 1, then both `a` and `b` are normalized by ``a0``. x : array_like The array of data to be filtered. axis : int, optional The axis of `x` to which the filter is applied. Default is -1. padtype : str or None, optional Must be 'odd', 'even', 'constant', or None. This determines the type of extension to use for the padded signal to which the filter is applied. If `padtype` is None, no padding is used. The default is 'odd'. padlen : int or None, optional The number of elements by which to extend `x` at both ends of `axis` before applying the filter. This value must be less than ``x.shapeaxis - 1``. ``padlen=0`` implies no padding. The default value is ``3 * max(len(a), len(b))``. method : str, optional Determines the method for handling the edges of the signal, either 'pad' or 'gust'. When `method` is 'pad', the signal is padded; the type of padding is determined by `padtype` and `padlen`, and `irlen` is ignored. When `method` is 'gust', Gustafsson's method is used, and `padtype` and `padlen` are ignored. irlen : int or None, optional When `method` is 'gust', `irlen` specifies the length of the impulse response of the filter. If `irlen` is None, no part of the impulse response is ignored. For a long signal, specifying `irlen` can significantly improve the performance of the filter.

Returns ------- y : ndarray The filtered output with the same shape as `x`.

See Also -------- sosfiltfilt, lfilter_zi, lfilter, lfiltic, savgol_filter, sosfilt

Notes ----- When `method` is 'pad', the function pads the data along the given axis in one of three ways: odd, even or constant. The odd and even extensions have the corresponding symmetry about the end point of the data. The constant extension extends the data with the values at the end points. On both the forward and backward passes, the initial condition of the filter is found by using `lfilter_zi` and scaling it by the end point of the extended data.

When `method` is 'gust', Gustafsson's method 1_ is used. Initial conditions are chosen for the forward and backward passes so that the forward-backward filter gives the same result as the backward-forward filter.

The option to use Gustaffson's method was added in scipy version 0.16.0.

References ---------- .. 1 F. Gustaffson, 'Determining the initial states in forward-backward filtering', Transactions on Signal Processing, Vol. 46, pp. 988-992, 1996.

Examples -------- The examples will use several functions from `scipy.signal`.

>>> from scipy import signal >>> import matplotlib.pyplot as plt

First we create a one second signal that is the sum of two pure sine waves, with frequencies 5 Hz and 250 Hz, sampled at 2000 Hz.

>>> t = np.linspace(0, 1.0, 2001) >>> xlow = np.sin(2 * np.pi * 5 * t) >>> xhigh = np.sin(2 * np.pi * 250 * t) >>> x = xlow + xhigh

Now create a lowpass Butterworth filter with a cutoff of 0.125 times the Nyquist frequency, or 125 Hz, and apply it to ``x`` with `filtfilt`. The result should be approximately ``xlow``, with no phase shift.

>>> b, a = signal.butter(8, 0.125) >>> y = signal.filtfilt(b, a, x, padlen=150) >>> np.abs(y - xlow).max() 9.1086182074789912e-06

We get a fairly clean result for this artificial example because the odd extension is exact, and with the moderately long padding, the filter's transients have dissipated by the time the actual data is reached. In general, transient effects at the edges are unavoidable.

The following example demonstrates the option ``method='gust'``.

First, create a filter.

>>> b, a = signal.ellip(4, 0.01, 120, 0.125) # Filter to be applied. >>> np.random.seed(123456)

`sig` is a random input signal to be filtered.

>>> n = 60 >>> sig = np.random.randn(n)**3 + 3*np.random.randn(n).cumsum()

Apply `filtfilt` to `sig`, once using the Gustafsson method, and once using padding, and plot the results for comparison.

>>> fgust = signal.filtfilt(b, a, sig, method='gust') >>> fpad = signal.filtfilt(b, a, sig, padlen=50) >>> plt.plot(sig, 'k-', label='input') >>> plt.plot(fgust, 'b-', linewidth=4, label='gust') >>> plt.plot(fpad, 'c-', linewidth=1.5, label='pad') >>> plt.legend(loc='best') >>> plt.show()

The `irlen` argument can be used to improve the performance of Gustafsson's method.

Estimate the impulse response length of the filter.

>>> z, p, k = signal.tf2zpk(b, a) >>> eps = 1e-9 >>> r = np.max(np.abs(p)) >>> approx_impulse_len = int(np.ceil(np.log(eps) / np.log(r))) >>> approx_impulse_len 137

Apply the filter to a longer signal, with and without the `irlen` argument. The difference between `y1` and `y2` is small. For long signals, using `irlen` gives a significant performance improvement.

>>> x = np.random.randn(5000) >>> y1 = signal.filtfilt(b, a, x, method='gust') >>> y2 = signal.filtfilt(b, a, x, method='gust', irlen=approx_impulse_len) >>> print(np.max(np.abs(y1 - y2))) 1.80056858312e-10

Find peaks inside a signal based on peak properties.

This function takes a 1-D array and finds all local maxima by simple comparison of neighboring values. Optionally, a subset of these peaks can be selected by specifying conditions for a peak's properties.

Parameters ---------- x : sequence A signal with peaks. height : number or ndarray or sequence, optional Required height of peaks. Either a number, ``None``, an array matching `x` or a 2-element sequence of the former. The first element is always interpreted as the minimal and the second, if supplied, as the maximal required height. threshold : number or ndarray or sequence, optional Required threshold of peaks, the vertical distance to its neighboring samples. Either a number, ``None``, an array matching `x` or a 2-element sequence of the former. The first element is always interpreted as the minimal and the second, if supplied, as the maximal required threshold. distance : number, optional Required minimal horizontal distance (>= 1) in samples between neighbouring peaks. Smaller peaks are removed first until the condition is fulfilled for all remaining peaks. prominence : number or ndarray or sequence, optional Required prominence of peaks. Either a number, ``None``, an array matching `x` or a 2-element sequence of the former. The first element is always interpreted as the minimal and the second, if supplied, as the maximal required prominence. width : number or ndarray or sequence, optional Required width of peaks in samples. Either a number, ``None``, an array matching `x` or a 2-element sequence of the former. The first element is always interpreted as the minimal and the second, if supplied, as the maximal required width. wlen : int, optional Used for calculation of the peaks prominences, thus it is only used if one of the arguments `prominence` or `width` is given. See argument `wlen` in `peak_prominences` for a full description of its effects. rel_height : float, optional Used for calculation of the peaks width, thus it is only used if `width` is given. See argument `rel_height` in `peak_widths` for a full description of its effects. plateau_size : number or ndarray or sequence, optional Required size of the flat top of peaks in samples. Either a number, ``None``, an array matching `x` or a 2-element sequence of the former. The first element is always interpreted as the minimal and the second, if supplied as the maximal required plateau size.

.. versionadded:: 1.2.0

Returns ------- peaks : ndarray Indices of peaks in `x` that satisfy all given conditions. properties : dict A dictionary containing properties of the returned peaks which were calculated as intermediate results during evaluation of the specified conditions:

* 'peak_heights' If `height` is given, the height of each peak in `x`. * 'left_thresholds', 'right_thresholds' If `threshold` is given, these keys contain a peaks vertical distance to its neighbouring samples. * 'prominences', 'right_bases', 'left_bases' If `prominence` is given, these keys are accessible. See `peak_prominences` for a description of their content. * 'width_heights', 'left_ips', 'right_ips' If `width` is given, these keys are accessible. See `peak_widths` for a description of their content. * 'plateau_sizes', left_edges', 'right_edges' If `plateau_size` is given, these keys are accessible and contain the indices of a peak's edges (edges are still part of the plateau) and the calculated plateau sizes.

.. versionadded:: 1.2.0

To calculate and return properties without excluding peaks, provide the open interval ``(None, None)`` as a value to the appropriate argument (excluding `distance`).

Warns ----- PeakPropertyWarning Raised if a peak's properties have unexpected values (see `peak_prominences` and `peak_widths`).

Warnings -------- This function may return unexpected results for data containing NaNs. To avoid this, NaNs should either be removed or replaced.

See Also -------- find_peaks_cwt Find peaks using the wavelet transformation. peak_prominences Directly calculate the prominence of peaks. peak_widths Directly calculate the width of peaks.

Notes ----- In the context of this function, a peak or local maximum is defined as any sample whose two direct neighbours have a smaller amplitude. For flat peaks (more than one sample of equal amplitude wide) the index of the middle sample is returned (rounded down in case the number of samples is even). For noisy signals the peak locations can be off because the noise might change the position of local maxima. In those cases consider smoothing the signal before searching for peaks or use other peak finding and fitting methods (like `find_peaks_cwt`).

Some additional comments on specifying conditions:

* Almost all conditions (excluding `distance`) can be given as half-open or closed intervals, e.g., ``1`` or ``(1, None)`` defines the half-open interval :math:`1, \infty` while ``(None, 1)`` defines the interval :math:`-\infty, 1`. The open interval ``(None, None)`` can be specified as well, which returns the matching properties without exclusion of peaks. * The border is always included in the interval used to select valid peaks. * For several conditions the interval borders can be specified with arrays matching `x` in shape which enables dynamic constrains based on the sample position. * The conditions are evaluated in the following order: `plateau_size`, `height`, `threshold`, `distance`, `prominence`, `width`. In most cases this order is the fastest one because faster operations are applied first to reduce the number of peaks that need to be evaluated later. * While indices in `peaks` are guaranteed to be at least `distance` samples apart, edges of flat peaks may be closer than the allowed `distance`. * Use `wlen` to reduce the time it takes to evaluate the conditions for `prominence` or `width` if `x` is large or has many local maxima (see `peak_prominences`).

.. versionadded:: 1.1.0

Examples -------- To demonstrate this function's usage we use a signal `x` supplied with SciPy (see `scipy.misc.electrocardiogram`). Let's find all peaks (local maxima) in `x` whose amplitude lies above 0.

>>> import matplotlib.pyplot as plt >>> from scipy.misc import electrocardiogram >>> from scipy.signal import find_peaks >>> x = electrocardiogram()2000:4000 >>> peaks, _ = find_peaks(x, height=0) >>> plt.plot(x) >>> plt.plot(peaks, xpeaks, 'x') >>> plt.plot(np.zeros_like(x), '--', color='gray') >>> plt.show()

We can select peaks below 0 with ``height=(None, 0)`` or use arrays matching `x` in size to reflect a changing condition for different parts of the signal.

>>> border = np.sin(np.linspace(0, 3 * np.pi, x.size)) >>> peaks, _ = find_peaks(x, height=(-border, border)) >>> plt.plot(x) >>> plt.plot(-border, '--', color='gray') >>> plt.plot(border, ':', color='gray') >>> plt.plot(peaks, xpeaks, 'x') >>> plt.show()

Another useful condition for periodic signals can be given with the `distance` argument. In this case, we can easily select the positions of QRS complexes within the electrocardiogram (ECG) by demanding a distance of at least 150 samples.

>>> peaks, _ = find_peaks(x, distance=150) >>> np.diff(peaks) array(186, 180, 177, 171, 177, 169, 167, 164, 158, 162, 172) >>> plt.plot(x) >>> plt.plot(peaks, xpeaks, 'x') >>> plt.show()

Especially for noisy signals peaks can be easily grouped by their prominence (see `peak_prominences`). E.g., we can select all peaks except for the mentioned QRS complexes by limiting the allowed prominence to 0.6.

>>> peaks, properties = find_peaks(x, prominence=(None, 0.6)) >>> properties'prominences'.max() 0.5049999999999999 >>> plt.plot(x) >>> plt.plot(peaks, xpeaks, 'x') >>> plt.show()

And, finally, let's examine a different section of the ECG which contains beat forms of different shape. To select only the atypical heart beats, we combine two conditions: a minimal prominence of 1 and width of at least 20 samples.

>>> x = electrocardiogram()17000:18000 >>> peaks, properties = find_peaks(x, prominence=1, width=20) >>> properties'prominences', properties'widths' (array(1.495, 2.3 ), array(36.93773946, 39.32723577)) >>> plt.plot(x) >>> plt.plot(peaks, xpeaks, 'x') >>> plt.vlines(x=peaks, ymin=xpeaks - properties'prominences', ... ymax = xpeaks, color = 'C1') >>> plt.hlines(y=properties'width_heights', xmin=properties'left_ips', ... xmax=properties'right_ips', color = 'C1') >>> plt.show()

Find peaks in a 1-D array with wavelet transformation.

The general approach is to smooth `vector` by convolving it with `wavelet(width)` for each width in `widths`. Relative maxima which appear at enough length scales, and with sufficiently high SNR, are accepted.

Parameters ---------- vector : ndarray 1-D array in which to find the peaks. widths : sequence 1-D array of widths to use for calculating the CWT matrix. In general, this range should cover the expected width of peaks of interest. wavelet : callable, optional Should take two parameters and return a 1-D array to convolve with `vector`. The first parameter determines the number of points of the returned wavelet array, the second parameter is the scale (`width`) of the wavelet. Should be normalized and symmetric. Default is the ricker wavelet. max_distances : ndarray, optional At each row, a ridge line is only connected if the relative max at rown is within ``max_distancesn`` from the relative max at ``rown+1``. Default value is ``widths/4``. gap_thresh : float, optional If a relative maximum is not found within `max_distances`, there will be a gap. A ridge line is discontinued if there are more than `gap_thresh` points without connecting a new relative maximum. Default is the first value of the widths array i.e. widths0. min_length : int, optional Minimum length a ridge line needs to be acceptable. Default is ``cwt.shape0 / 4``, ie 1/4-th the number of widths. min_snr : float, optional Minimum SNR ratio. Default 1. The signal is the value of the cwt matrix at the shortest length scale (``cwt0, loc``), the noise is the `noise_perc`th percentile of datapoints contained within a window of `window_size` around ``cwt0, loc``. noise_perc : float, optional When calculating the noise floor, percentile of data points examined below which to consider noise. Calculated using `stats.scoreatpercentile`. Default is 10. window_size : int, optional Size of window to use to calculate noise floor. Default is ``cwt.shape1 / 20``.

Returns ------- peaks_indices : ndarray Indices of the locations in the `vector` where peaks were found. The list is sorted.

See Also -------- cwt Continuous wavelet transform. find_peaks Find peaks inside a signal based on peak properties.

Notes ----- This approach was designed for finding sharp peaks among noisy data, however with proper parameter selection it should function well for different peak shapes.

The algorithm is as follows: 1. Perform a continuous wavelet transform on `vector`, for the supplied `widths`. This is a convolution of `vector` with `wavelet(width)` for each width in `widths`. See `cwt`. 2. Identify 'ridge lines' in the cwt matrix. These are relative maxima at each row, connected across adjacent rows. See identify_ridge_lines 3. Filter the ridge_lines using filter_ridge_lines.

.. versionadded:: 0.11.0

References ---------- .. 1 Bioinformatics (2006) 22 (17): 2059-2065. :doi:`10.1093/bioinformatics/btl355`

Examples -------- >>> from scipy import signal >>> xs = np.arange(0, np.pi, 0.05) >>> data = np.sin(xs) >>> peakind = signal.find_peaks_cwt(data, np.arange(1,10)) >>> peakind, xspeakind, datapeakind (32, array( 1.6), array( 0.9995736))

Find array of frequencies for computing the response of an analog filter.

Parameters ---------- num, den : array_like, 1-D The polynomial coefficients of the numerator and denominator of the transfer function of the filter or LTI system, where the coefficients are ordered from highest to lowest degree. Or, the roots of the transfer function numerator and denominator (i.e., zeroes and poles). N : int The length of the array to be computed. kind : str 'ba', 'zp', optional Specifies whether the numerator and denominator are specified by their polynomial coefficients ('ba'), or their roots ('zp').

Returns ------- w : (N,) ndarray A 1-D array of frequencies, logarithmically spaced.

Examples -------- Find a set of nine frequencies that span the 'interesting part' of the frequency response for the filter with the transfer function

H(s) = s / (s^2 + 8s + 25)

>>> from scipy import signal >>> signal.findfreqs(1, 0, 1, 8, 25, N=9) array( 1.00000000e-02, 3.16227766e-02, 1.00000000e-01, 3.16227766e-01, 1.00000000e+00, 3.16227766e+00, 1.00000000e+01, 3.16227766e+01, 1.00000000e+02)

FIR filter design using least-squares error minimization.

Calculate the filter coefficients for the linear-phase finite impulse response (FIR) filter which has the best approximation to the desired frequency response described by `bands` and `desired` in the least squares sense (i.e., the integral of the weighted mean-squared error within the specified bands is minimized).

Parameters ---------- numtaps : int The number of taps in the FIR filter. `numtaps` must be odd. bands : array_like A monotonic nondecreasing sequence containing the band edges in Hz. All elements must be non-negative and less than or equal to the Nyquist frequency given by `nyq`. desired : array_like A sequence the same size as `bands` containing the desired gain at the start and end point of each band. weight : array_like, optional A relative weighting to give to each band region when solving the least squares problem. `weight` has to be half the size of `bands`. nyq : float, optional *Deprecated. Use `fs` instead.* Nyquist frequency. Each frequency in `bands` must be between 0 and `nyq` (inclusive). Default is 1. fs : float, optional The sampling frequency of the signal. Each frequency in `bands` must be between 0 and ``fs/2`` (inclusive). Default is 2.

Returns ------- coeffs : ndarray Coefficients of the optimal (in a least squares sense) FIR filter.

See also -------- firwin firwin2 minimum_phase remez

Notes ----- This implementation follows the algorithm given in 1_. As noted there, least squares design has multiple advantages:

1. Optimal in a least-squares sense. 2. Simple, non-iterative method. 3. The general solution can obtained by solving a linear system of equations. 4. Allows the use of a frequency dependent weighting function.

This function constructs a Type I linear phase FIR filter, which contains an odd number of `coeffs` satisfying for :math:`n < numtaps`:

.. math:: coeffs(n) = coeffs(numtaps - 1 - n)

The odd number of coefficients and filter symmetry avoid boundary conditions that could otherwise occur at the Nyquist and 0 frequencies (e.g., for Type II, III, or IV variants).

.. versionadded:: 0.18

References ---------- .. 1 Ivan Selesnick, Linear-Phase Fir Filter Design By Least Squares. OpenStax CNX. Aug 9, 2005. http://cnx.org/contents/eb1ecb35-03a9-4610-ba87-41cd771c95f2@7

Examples -------- We want to construct a band-pass filter. Note that the behavior in the frequency ranges between our stop bands and pass bands is unspecified, and thus may overshoot depending on the parameters of our filter:

>>> from scipy import signal >>> import matplotlib.pyplot as plt >>> fig, axs = plt.subplots(2) >>> fs = 10.0 # Hz >>> desired = (0, 0, 1, 1, 0, 0) >>> for bi, bands in enumerate(((0, 1, 2, 3, 4, 5), (0, 1, 2, 4, 4.5, 5))): ... fir_firls = signal.firls(73, bands, desired, fs=fs) ... fir_remez = signal.remez(73, bands, desired::2, fs=fs) ... fir_firwin2 = signal.firwin2(73, bands, desired, fs=fs) ... hs = list() ... ax = axsbi ... for fir in (fir_firls, fir_remez, fir_firwin2): ... freq, response = signal.freqz(fir) ... hs.append(ax.semilogy(0.5*fs*freq/np.pi, np.abs(response))0) ... for band, gains in zip(zip(bands::2, bands1::2), ... zip(desired::2, desired1::2)): ... ax.semilogy(band, np.maximum(gains, 1e-7), 'k--', linewidth=2) ... if bi == 0: ... ax.legend(hs, ('firls', 'remez', 'firwin2'), ... loc='lower center', frameon=False) ... else: ... ax.set_xlabel('Frequency (Hz)') ... ax.grid(True) ... ax.set(title='Band-pass %d-%d Hz' % bands2:4, ylabel='Magnitude') ... >>> fig.tight_layout() >>> plt.show()

FIR filter design using the window method.

This function computes the coefficients of a finite impulse response filter. The filter will have linear phase; it will be Type I if `numtaps` is odd and Type II if `numtaps` is even.

Type II filters always have zero response at the Nyquist frequency, so a ValueError exception is raised if firwin is called with `numtaps` even and having a passband whose right end is at the Nyquist frequency.

Parameters ---------- numtaps : int Length of the filter (number of coefficients, i.e. the filter order + 1). `numtaps` must be odd if a passband includes the Nyquist frequency. cutoff : float or 1-D array_like Cutoff frequency of filter (expressed in the same units as `fs`) OR an array of cutoff frequencies (that is, band edges). In the latter case, the frequencies in `cutoff` should be positive and monotonically increasing between 0 and `fs/2`. The values 0 and `fs/2` must not be included in `cutoff`. width : float or None, optional If `width` is not None, then assume it is the approximate width of the transition region (expressed in the same units as `fs`) for use in Kaiser FIR filter design. In this case, the `window` argument is ignored. window : string or tuple of string and parameter values, optional Desired window to use. See `scipy.signal.get_window` for a list of windows and required parameters. pass_zero : True, False, 'bandpass', 'lowpass', 'highpass', 'bandstop', optional If True, the gain at the frequency 0 (i.e., the 'DC gain') is 1. If False, the DC gain is 0. Can also be a string argument for the desired filter type (equivalent to ``btype`` in IIR design functions).

.. versionadded:: 1.3.0 Support for string arguments. scale : bool, optional Set to True to scale the coefficients so that the frequency response is exactly unity at a certain frequency. That frequency is either:

• 0 (DC) if the first passband starts at 0 (i.e. pass_zero is True)
• `fs/2` (the Nyquist frequency) if the first passband ends at `fs/2` (i.e the filter is a single band highpass filter); center of first passband otherwise

nyq : float, optional *Deprecated. Use `fs` instead.* This is the Nyquist frequency. Each frequency in `cutoff` must be between 0 and `nyq`. Default is 1. fs : float, optional The sampling frequency of the signal. Each frequency in `cutoff` must be between 0 and ``fs/2``. Default is 2.

Returns ------- h : (numtaps,) ndarray Coefficients of length `numtaps` FIR filter.

Raises ------ ValueError If any value in `cutoff` is less than or equal to 0 or greater than or equal to ``fs/2``, if the values in `cutoff` are not strictly monotonically increasing, or if `numtaps` is even but a passband includes the Nyquist frequency.

See Also -------- firwin2 firls minimum_phase remez

Examples -------- Low-pass from 0 to f:

>>> from scipy import signal >>> numtaps = 3 >>> f = 0.1 >>> signal.firwin(numtaps, f) array( 0.06799017, 0.86401967, 0.06799017)

Use a specific window function:

>>> signal.firwin(numtaps, f, window='nuttall') array( 3.56607041e-04, 9.99286786e-01, 3.56607041e-04)

High-pass ('stop' from 0 to f):

>>> signal.firwin(numtaps, f, pass_zero=False) array(-0.00859313, 0.98281375, -0.00859313)

Band-pass:

>>> f1, f2 = 0.1, 0.2 >>> signal.firwin(numtaps, f1, f2, pass_zero=False) array( 0.06301614, 0.88770441, 0.06301614)

Band-stop:

>>> signal.firwin(numtaps, f1, f2) array(-0.00801395, 1.0160279 , -0.00801395)

Multi-band (passbands are 0, f1, f2, f3 and f4, 1):

>>> f3, f4 = 0.3, 0.4 >>> signal.firwin(numtaps, f1, f2, f3, f4) array(-0.01376344, 1.02752689, -0.01376344)

Multi-band (passbands are f1, f2 and f3,f4):

>>> signal.firwin(numtaps, f1, f2, f3, f4, pass_zero=False) array( 0.04890915, 0.91284326, 0.04890915)

FIR filter design using the window method.

From the given frequencies `freq` and corresponding gains `gain`, this function constructs an FIR filter with linear phase and (approximately) the given frequency response.

Parameters ---------- numtaps : int The number of taps in the FIR filter. `numtaps` must be less than `nfreqs`. freq : array_like, 1-D The frequency sampling points. Typically 0.0 to 1.0 with 1.0 being Nyquist. The Nyquist frequency is half `fs`. The values in `freq` must be nondecreasing. A value can be repeated once to implement a discontinuity. The first value in `freq` must be 0, and the last value must be ``fs/2``. Values 0 and ``fs/2`` must not be repeated. gain : array_like The filter gains at the frequency sampling points. Certain constraints to gain values, depending on the filter type, are applied, see Notes for details. nfreqs : int, optional The size of the interpolation mesh used to construct the filter. For most efficient behavior, this should be a power of 2 plus 1 (e.g, 129, 257, etc). The default is one more than the smallest power of 2 that is not less than `numtaps`. `nfreqs` must be greater than `numtaps`. window : string or (string, float) or float, or None, optional Window function to use. Default is 'hamming'. See `scipy.signal.get_window` for the complete list of possible values. If None, no window function is applied. nyq : float, optional *Deprecated. Use `fs` instead.* This is the Nyquist frequency. Each frequency in `freq` must be between 0 and `nyq`. Default is 1. antisymmetric : bool, optional Whether resulting impulse response is symmetric/antisymmetric. See Notes for more details. fs : float, optional The sampling frequency of the signal. Each frequency in `cutoff` must be between 0 and ``fs/2``. Default is 2.

Returns ------- taps : ndarray The filter coefficients of the FIR filter, as a 1-D array of length `numtaps`.

See also -------- firls firwin minimum_phase remez

Notes ----- From the given set of frequencies and gains, the desired response is constructed in the frequency domain. The inverse FFT is applied to the desired response to create the associated convolution kernel, and the first `numtaps` coefficients of this kernel, scaled by `window`, are returned.

The FIR filter will have linear phase. The type of filter is determined by the value of 'numtaps` and `antisymmetric` flag. There are four possible combinations:

• odd `numtaps`, `antisymmetric` is False, type I filter is produced
• even `numtaps`, `antisymmetric` is False, type II filter is produced
• odd `numtaps`, `antisymmetric` is True, type III filter is produced
• even `numtaps`, `antisymmetric` is True, type IV filter is produced

Magnitude response of all but type I filters are subjects to following constraints:

• type II -- zero at the Nyquist frequency
• type III -- zero at zero and Nyquist frequencies
• type IV -- zero at zero frequency

.. versionadded:: 0.9.0

References ---------- .. 1 Oppenheim, A. V. and Schafer, R. W., 'Discrete-Time Signal Processing', Prentice-Hall, Englewood Cliffs, New Jersey (1989). (See, for example, Section 7.4.)

.. 2 Smith, Steven W., 'The Scientist and Engineer's Guide to Digital Signal Processing', Ch. 17. http://www.dspguide.com/ch17/1.htm

Examples -------- A lowpass FIR filter with a response that is 1 on 0.0, 0.5, and that decreases linearly on 0.5, 1.0 from 1 to 0:

>>> from scipy import signal >>> taps = signal.firwin2(150, 0.0, 0.5, 1.0, 1.0, 1.0, 0.0) >>> print(taps72:78) -0.02286961 -0.06362756 0.57310236 0.57310236 -0.06362756 -0.02286961

Return a flat top window.

.. warning:: scipy.signal.flattop is deprecated, use scipy.signal.windows.flattop instead.

Parameters ---------- M : int Number of points in the output window. If zero or less, an empty array is returned. sym : bool, optional When True (default), generates a symmetric window, for use in filter design. When False, generates a periodic window, for use in spectral analysis.

Returns ------- w : ndarray The window, with the maximum value normalized to 1 (though the value 1 does not appear if `M` is even and `sym` is True).

Notes ----- Flat top windows are used for taking accurate measurements of signal amplitude in the frequency domain, with minimal scalloping error from the center of a frequency bin to its edges, compared to others. This is a 5th-order cosine window, with the 5 terms optimized to make the main lobe maximally flat. 1_

References ---------- .. 1 D'Antona, Gabriele, and A. Ferrero, 'Digital Signal Processing for Measurement Systems', Springer Media, 2006, p. 70 :doi:`10.1007/0-387-28666-7`.

Examples -------- Plot the window and its frequency response:

>>> from scipy import signal >>> from scipy.fft import fft, fftshift >>> import matplotlib.pyplot as plt

>>> window = signal.flattop(51) >>> plt.plot(window) >>> plt.title('Flat top window') >>> plt.ylabel('Amplitude') >>> plt.xlabel('Sample')

>>> plt.figure() >>> A = fft(window, 2048) / (len(window)/2.0) >>> freq = np.linspace(-0.5, 0.5, len(A)) >>> response = 20 * np.log10(np.abs(fftshift(A / abs(A).max()))) >>> plt.plot(freq, response) >>> plt.axis(-0.5, 0.5, -120, 0) >>> plt.title('Frequency response of the flat top window') >>> plt.ylabel('Normalized magnitude dB') >>> plt.xlabel('Normalized frequency cycles per sample')

Calculate the frequency response of a continuous-time system.

Parameters ---------- system : an instance of the `lti` class or a tuple describing the system. The following gives the number of elements in the tuple and the interpretation:

* 1 (instance of `lti`) * 2 (num, den) * 3 (zeros, poles, gain) * 4 (A, B, C, D)

w : array_like, optional Array of frequencies (in rad/s). Magnitude and phase data is calculated for every value in this array. If not given, a reasonable set will be calculated. n : int, optional Number of frequency points to compute if `w` is not given. The `n` frequencies are logarithmically spaced in an interval chosen to include the influence of the poles and zeros of the system.

Returns ------- w : 1D ndarray Frequency array rad/s H : 1D ndarray Array of complex magnitude values

Notes ----- If (num, den) is passed in for ``system``, coefficients for both the numerator and denominator should be specified in descending exponent order (e.g. ``s^2 + 3s + 5`` would be represented as ``1, 3, 5``).

Examples -------- Generating the Nyquist plot of a transfer function

>>> from scipy import signal >>> import matplotlib.pyplot as plt

Transfer function: H(s) = 5 / (s-1)^3

>>> s1 = signal.ZerosPolesGain(, 1, 1, 1, 5)

>>> w, H = signal.freqresp(s1)

>>> plt.figure() >>> plt.plot(H.real, H.imag, 'b') >>> plt.plot(H.real, -H.imag, 'r') >>> plt.show()

Compute frequency response of analog filter.

Given the M-order numerator `b` and N-order denominator `a` of an analog filter, compute its frequency response::

b0*(jw)**M + b1*(jw)**(M-1) + ... + bM H(w) = ---------------------------------------------- a0*(jw)**N + a1*(jw)**(N-1) + ... + aN

Parameters ---------- b : array_like Numerator of a linear filter. a : array_like Denominator of a linear filter. worN : None, int, array_like, optional If None, then compute at 200 frequencies around the interesting parts of the response curve (determined by pole-zero locations). If a single integer, then compute at that many frequencies. Otherwise, compute the response at the angular frequencies (e.g., rad/s) given in `worN`. plot : callable, optional A callable that takes two arguments. If given, the return parameters `w` and `h` are passed to plot. Useful for plotting the frequency response inside `freqs`.