image-match/python/py/Lib/site-packages/scipy/signal/windows/_windows.py

"""The suite of window functions."""

import math
import numbers
import operator
import warnings
from scipy._lib import doccer

from scipy import linalg, special, fft as sp_fft
from scipy._lib.array_api_compat import numpy as np_compat
from scipy._lib._array_api import array_namespace, xp_device
from scipy._lib._array_api import xp_capabilities
from scipy._lib import array_api_extra as xpx

__all__ = ['boxcar', 'triang', 'parzen', 'bohman', 'blackman', 'nuttall',
           'blackmanharris', 'flattop', 'bartlett', 'barthann',
           'hamming', 'kaiser', 'kaiser_bessel_derived', 'gaussian',
           'general_cosine', 'general_gaussian', 'general_hamming',
           'chebwin', 'cosine', 'hann', 'exponential', 'tukey', 'taylor',
           'dpss', 'get_window', 'lanczos']


def _len_guards(M):
    """Handle small or incorrect window lengths"""
    if int(M) != M or M < 0:
        raise ValueError('Window length M must be a non-negative integer')
    return M <= 1


def _extend(M, sym):
    """Extend window by 1 sample if needed for DFT-even symmetry"""
    if not sym:
        return M + 1, True
    else:
        return M, False


def _truncate(w, needed):
    """Truncate window by 1 sample if needed for DFT-even symmetry"""
    if needed:
        return w[:-1]
    else:
        return w


def _namespace(xp):
    """A shim for the `device` arg of `np.asarray(x, device=device)` and acos/arccos.

    Will be able to replace with `np_compat if xp is None else xp` when we drop
    support for numpy 1.x and cupy 13.x
    """
    return np_compat if xp is None else array_namespace(xp.empty(0))


def _general_cosine_impl(M, a, xp, device, sym=True):
    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    fac = xp.linspace(-xp.pi, xp.pi, M, dtype=xp.float64, device=device)
    w = xp.zeros(M, dtype=xp.float64, device=device)
    for k in range(a.shape[0]):
        w += a[k] * xp.cos(k * fac)

    return _truncate(w, needs_trunc)


@xp_capabilities()
def general_cosine(M, a, sym=True):
    r"""
    Generic weighted sum of cosine terms window

    Parameters
    ----------
    M : int
        Number of points in the output window
    a : array_like
        Sequence of weighting coefficients. This uses the convention of being
        centered on the origin, so these will typically all be positive
        numbers, not alternating sign.
    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 array of window values.

    References
    ----------
    .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE
           Transactions on Acoustics, Speech, and Signal Processing, vol. 29,
           no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
    .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the
           Discrete Fourier transform (DFT), including a comprehensive list of
           window functions and some new flat-top windows", February 15, 2002
           https://holometer.fnal.gov/GH_FFT.pdf

    Examples
    --------
    Heinzel describes a flat-top window named "HFT90D" with formula: [2]_

    .. math::  w_j = 1 - 1.942604 \cos(z) + 1.340318 \cos(2z)
               - 0.440811 \cos(3z) + 0.043097 \cos(4z)

    where

    .. math::  z = \frac{2 \pi j}{N}, j = 0...N - 1

    Since this uses the convention of starting at the origin, to reproduce the
    window, we need to convert every other coefficient to a positive number:

    >>> HFT90D = [1, 1.942604, 1.340318, 0.440811, 0.043097]

    The paper states that the highest sidelobe is at -90.2 dB.  Reproduce
    Figure 42 by plotting the window and its frequency response, and confirm
    the sidelobe level in red:

    >>> import numpy as np
    >>> from scipy.signal.windows import general_cosine
    >>> from scipy.fft import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> window = general_cosine(1000, HFT90D, sym=False)
    >>> plt.plot(window)
    >>> plt.title("HFT90D window")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")

    >>> plt.figure()
    >>> A = fft(window, 10000) / (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([-50/1000, 50/1000, -140, 0])
    >>> plt.title("Frequency response of the HFT90D window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")
    >>> plt.axhline(-90.2, color='red')
    >>> plt.show()
    """
    xp = array_namespace(a)
    a = xp.asarray(a)
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def boxcar(M, sym=True, *, xp=None, device=None):
    """Return a boxcar or rectangular window.

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

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    sym : bool, optional
        Whether the window is symmetric. (Has no effect for boxcar.)
    %(xp_device_snippet)s

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

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

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

    >>> window = signal.windows.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]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    w = xp.ones(M, dtype=xp.float64, device=device)

    return _truncate(w, needs_trunc)


@xp_capabilities()
def triang(M, sym=True, *, xp=None, device=None):
    """Return a triangular window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

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

    See Also
    --------
    bartlett : A triangular window that touches zero

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

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

    >>> window = signal.windows.triang(51)
    >>> plt.plot(window)
    >>> plt.title("Triangular 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 triangular window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(1, (M + 1) // 2 + 1, dtype=xp.float64, device=device)
    if M % 2 == 0:
        w = (2 * n - 1.0) / M
        w = xp.concat([w, xp.flip(w)])
    else:
        w = 2 * n / (M + 1.0)
        w = xp.concat([w, xp.flip(w[:-1])])

    return _truncate(w, needs_trunc)


@xp_capabilities()
def parzen(M, sym=True, *, xp=None, device=None):
    """Return a Parzen window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

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

    References
    ----------
    .. [1] E. Parzen, "Mathematical Considerations in the Estimation of
           Spectra", Technometrics,  Vol. 3, No. 2 (May, 1961), pp. 167-190

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

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

    >>> window = signal.windows.parzen(51)
    >>> plt.plot(window)
    >>> plt.title("Parzen 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 Parzen window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(-(M - 1) / 2.0, (M - 1) / 2.0 + 0.5, 1.0,
                  dtype=xp.float64, device=device)
    w = xp.where(abs(n) <= (M - 1) / 4.0,
                 (1 - 6 * (abs(n) / (M / 2.0)) ** 2.0 +
                  6 * (abs(n) / (M / 2.0)) ** 3.0),
                 2 * (1 - abs(n) / (M / 2.0)) ** 3.0)
    return _truncate(w, needs_trunc)


@xp_capabilities()
def bohman(M, sym=True, *, xp=None, device=None):
    """Return a Bohman window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

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

    >>> plt.figure()
    >>> A = fft(window, 2047) / (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]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    fac = abs(xp.linspace(-1, 1, M, dtype=xp.float64, device=device)[1:-1])
    w = (1 - fac) * xp.cos(xp.pi * fac) + 1.0 / xp.pi * xp.sin(xp.pi * fac)
    one = xp.zeros(1, dtype=xp.float64, device=device)
    w = xp.concat([one, w, one])

    return _truncate(w, needs_trunc)


@xp_capabilities()
def blackman(M, sym=True, *, xp=None, device=None):
    r"""
    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.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

    >>> window = signal.windows.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]")

    """
    # Docstring adapted from NumPy's blackman function
    xp = _namespace(xp)
    a = xp.asarray([0.42, 0.50, 0.08], dtype=xp.float64, device=device)
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def nuttall(M, sym=True, *, xp=None, device=None):
    """Return a minimum 4-term Blackman-Harris window according to Nuttall.

    This variation is called "Nuttall4c" by Heinzel. [2]_

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

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

    References
    ----------
    .. [1] A. Nuttall, "Some windows with very good sidelobe behavior," IEEE
           Transactions on Acoustics, Speech, and Signal Processing, vol. 29,
           no. 1, pp. 84-91, Feb 1981. :doi:`10.1109/TASSP.1981.1163506`.
    .. [2] Heinzel G. et al., "Spectrum and spectral density estimation by the
           Discrete Fourier transform (DFT), including a comprehensive list of
           window functions and some new flat-top windows", February 15, 2002
           https://holometer.fnal.gov/GH_FFT.pdf

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

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

    >>> window = signal.windows.nuttall(51)
    >>> plt.plot(window)
    >>> plt.title("Nuttall 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 Nuttall window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)
    a = xp.asarray(
        [0.3635819, 0.4891775, 0.1365995, 0.0106411], dtype=xp.float64, device=device
    )
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def blackmanharris(M, sym=True, *, xp=None, device=None):
    """Return a minimum 4-term Blackman-Harris window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

    >>> window = signal.windows.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]")

    """
    xp = _namespace(xp)
    a = xp.asarray(
        [0.35875, 0.48829, 0.14128, 0.01168], dtype=xp.float64, device=device
    )
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def flattop(M, sym=True, *, xp=None, device=None):
    """Return a flat top window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

    >>> window = signal.windows.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]")

    """
    xp = _namespace(xp)
    a = xp.asarray(
        [0.21557895, 0.41663158, 0.277263158, 0.083578947, 0.006947368],
        dtype=xp.float64, device=device
    )
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def bartlett(M, sym=True, *, xp=None, device=None):
    r"""
    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.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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{2}{M-1} \left(
              \frac{M-1}{2} - \left|n - \frac{M-1}{2}\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:

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

    >>> window = signal.windows.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]")

    """
    # Docstring adapted from NumPy's bartlett function
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device)

    # cf https://github.com/data-apis/array-api-strict/issues/77
    w = xp.where(n <= (M - 1) / 2.0,
                 2.0 * n / (M - 1), 2.0 - 2.0 * n / (M - 1))

    return _truncate(w, needs_trunc)


@xp_capabilities()
def hann(M, sym=True, *, xp=None, device=None):
    r"""
    Return a Hann window.

    The Hann window is a taper formed by using a raised cosine or sine-squared
    with ends that touch zero.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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 Hann window is defined as

    .. math::  w(n) = 0.5 - 0.5 \cos\left(\frac{2\pi{n}}{M-1}\right)
               \qquad 0 \leq n \leq M-1

    The window was named for Julius von Hann, an Austrian meteorologist. It is
    also known as the Cosine Bell. It is sometimes erroneously referred to as
    the "Hanning" window, from the use of "hann" as a verb in the original
    paper and confusion with the very similar Hamming window.

    Most references to the Hann 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.

    References
    ----------
    .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
           spectra, Dover Publications, New York.
    .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics",
           The University of Alberta Press, 1975, pp. 106-108.
    .. [3] Wikipedia, "Window function",
           https://en.wikipedia.org/wiki/Window_function
    .. [4] W.H. Press,  B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
           "Numerical Recipes", Cambridge University Press, 1986, page 425.

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

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

    >>> window = signal.windows.hann(51)
    >>> plt.plot(window)
    >>> plt.title("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 = 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 Hann window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    # Docstring adapted from NumPy's hanning function
    return general_hamming(M, 0.5, sym, xp=xp, device=device)


@xp_capabilities()
def tukey(M, alpha=0.5, sym=True, *, xp=None, device=None):
    r"""Return a Tukey window, also known as a tapered cosine window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    alpha : float, optional
        Shape parameter of the Tukey window, representing the fraction of the
        window inside the cosine tapered region.
        If zero, the Tukey window is equivalent to a rectangular window.
        If one, the Tukey window is equivalent to a Hann window.
    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.
    %(xp_device_snippet)s

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

    References
    ----------
    .. [1] 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`
    .. [2] Wikipedia, "Window function",
           https://en.wikipedia.org/wiki/Window_function#Tukey_window

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

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

    >>> window = signal.windows.tukey(51)
    >>> plt.plot(window)
    >>> plt.title("Tukey window")
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")
    >>> plt.ylim([0, 1.1])

    >>> 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 Tukey window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)

    if alpha <= 0:
        return xp.ones(M, dtype=xp.float64, device=device)
    elif alpha >= 1.0:
        return hann(M, sym=sym, xp=xp, device=device)

    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device)
    width = int(math.floor(alpha*(M-1)/2.0))
    n1 = n[0:width+1]
    n2 = n[width+1:M-width-1]
    n3 = n[M-width-1:]

    w1 = 0.5 * (1 + xp.cos(xp.pi * (-1 + 2.0*n1/alpha/(M-1))))
    w2 = xp.ones(n2.shape, device=device)
    w3 = 0.5 * (1 + xp.cos(xp.pi * (-2.0/alpha + 1 + 2.0*n3/alpha/(M-1))))

    w = xp.concat((w1, w2, w3))

    return _truncate(w, needs_trunc)


@xp_capabilities()
def barthann(M, sym=True, *, xp=None, device=None):
    """Return a modified Bartlett-Hann window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

    >>> window = signal.windows.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]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device)
    fac = abs(n / (M - 1.0) - 0.5)
    w = 0.62 - 0.48 * fac + 0.38 * xp.cos(2 * xp.pi * fac)

    return _truncate(w, needs_trunc)


@xp_capabilities()
def general_hamming(M, alpha, sym=True, *, xp=None, device=None):
    r"""Return a generalized Hamming window.

    The generalized Hamming window is constructed by multiplying a rectangular
    window by one period of a cosine function [1]_.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    alpha : float
        The window coefficient, :math:`\alpha`
    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.
    %(xp_device_snippet)s

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

    See Also
    --------
    hamming, hann

    Notes
    -----
    The generalized Hamming window is defined as

    .. math:: w(n) = \alpha - \left(1 - \alpha\right)
              \cos\left(\frac{2\pi{n}}{M-1}\right) \qquad 0 \leq n \leq M-1

    Both the common Hamming window and Hann window are special cases of the
    generalized Hamming window with :math:`\alpha` = 0.54 and :math:`\alpha` =
    0.5, respectively [2]_.

    References
    ----------
    .. [1] DSPRelated, "Generalized Hamming Window Family",
           https://www.dsprelated.com/freebooks/sasp/Generalized_Hamming_Window_Family.html
    .. [2] Wikipedia, "Window function",
           https://en.wikipedia.org/wiki/Window_function
    .. [3] Riccardo Piantanida ESA, "Sentinel-1 Level 1 Detailed Algorithm
           Definition",
           https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Level-1-Detailed-Algorithm-Definition
    .. [4] Matthieu Bourbigot ESA, "Sentinel-1 Product Definition",
           https://sentinel.esa.int/documents/247904/1877131/Sentinel-1-Product-Definition

    Examples
    --------
    The Sentinel-1A/B Instrument Processing Facility uses generalized Hamming
    windows in the processing of spaceborne Synthetic Aperture Radar (SAR)
    data [3]_. The facility uses various values for the :math:`\alpha`
    parameter based on operating mode of the SAR instrument. Some common
    :math:`\alpha` values include 0.75, 0.7 and 0.52 [4]_. As an example, we
    plot these different windows.

    >>> import numpy as np
    >>> from scipy.signal.windows import general_hamming
    >>> from scipy.fft import fft, fftshift
    >>> import matplotlib.pyplot as plt

    >>> fig1, spatial_plot = plt.subplots()
    >>> spatial_plot.set_title("Generalized Hamming Windows")
    >>> spatial_plot.set_ylabel("Amplitude")
    >>> spatial_plot.set_xlabel("Sample")

    >>> fig2, freq_plot = plt.subplots()
    >>> freq_plot.set_title("Frequency Responses")
    >>> freq_plot.set_ylabel("Normalized magnitude [dB]")
    >>> freq_plot.set_xlabel("Normalized frequency [cycles per sample]")

    >>> for alpha in [0.75, 0.7, 0.52]:
    ...     window = general_hamming(41, alpha)
    ...     spatial_plot.plot(window, label="{:.2f}".format(alpha))
    ...     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())))
    ...     freq_plot.plot(freq, response, label="{:.2f}".format(alpha))
    >>> freq_plot.legend(loc="upper right")
    >>> spatial_plot.legend(loc="upper right")

    """
    xp = _namespace(xp)
    a = xp.asarray([alpha, 1. - alpha], dtype=xp.float64, device=device)
    device = xp_device(a)
    return _general_cosine_impl(M, a, xp, device, sym=sym)


@xp_capabilities()
def hamming(M, sym=True, *, xp=None, device=None):
    r"""Return a Hamming window.

    The Hamming window is a taper formed by using a raised cosine with
    non-zero endpoints, optimized to minimize the nearest side lobe.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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 Hamming window is defined as

    .. math::  w(n) = 0.54 - 0.46 \cos\left(\frac{2\pi{n}}{M-1}\right)
               \qquad 0 \leq n \leq M-1

    The Hamming was named for R. W. Hamming, an associate of J. W. Tukey and
    is described in Blackman and Tukey. It was recommended for smoothing the
    truncated autocovariance function in the time domain.
    Most references to the Hamming 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.

    References
    ----------
    .. [1] Blackman, R.B. and Tukey, J.W., (1958) The measurement of power
           spectra, Dover Publications, New York.
    .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
           University of Alberta Press, 1975, pp. 109-110.
    .. [3] Wikipedia, "Window function",
           https://en.wikipedia.org/wiki/Window_function
    .. [4] W.H. Press,  B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling,
           "Numerical Recipes", Cambridge University Press, 1986, page 425.

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

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

    >>> window = signal.windows.hamming(51)
    >>> plt.plot(window)
    >>> plt.title("Hamming 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 Hamming window")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    # Docstring adapted from NumPy's hamming function
    return general_hamming(M, 0.54, sym, xp=xp, device=device)


@xp_capabilities()
def kaiser(M, beta, sym=True, *, xp=None, device=None):
    r"""Return a Kaiser window.

    The Kaiser window is a taper formed by using a Bessel function.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    beta : float
        Shape parameter, determines trade-off between main-lobe width and
        side lobe level. As beta gets large, the window narrows.
    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.
    %(xp_device_snippet)s

    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 Kaiser window is defined as

    .. math::  w(n) = I_0\left( \beta \sqrt{1-\frac{4n^2}{(M-1)^2}}
               \right)/I_0(\beta)

    with

    .. math:: \quad -\frac{M-1}{2} \leq n \leq \frac{M-1}{2},

    where :math:`I_0` is the modified zeroth-order Bessel function.

    The Kaiser was named for Jim Kaiser, who discovered a simple approximation
    to the DPSS window based on Bessel functions.
    The Kaiser window is a very good approximation to the discrete prolate
    spheroidal sequence, or Slepian window, which is the transform which
    maximizes the energy in the main lobe of the window relative to total
    energy.

    The Kaiser can approximate other windows by varying the beta parameter.
    (Some literature uses alpha = beta/pi.) [4]_

    ====  =======================
    beta  Window shape
    ====  =======================
    0     Rectangular
    5     Similar to a Hamming
    6     Similar to a Hann
    8.6   Similar to a Blackman
    ====  =======================

    A beta value of 14 is probably a good starting point. Note that as beta
    gets large, the window narrows, and so the number of samples needs to be
    large enough to sample the increasingly narrow spike, otherwise NaNs will
    be returned.

    Most references to the Kaiser 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.

    References
    ----------
    .. [1] J. F. Kaiser, "Digital Filters" - Ch 7 in "Systems analysis by
           digital computer", Editors: F.F. Kuo and J.F. Kaiser, p 218-285.
           John Wiley and Sons, New York, (1966).
    .. [2] E.R. Kanasewich, "Time Sequence Analysis in Geophysics", The
           University of Alberta Press, 1975, pp. 177-178.
    .. [3] Wikipedia, "Window function",
           https://en.wikipedia.org/wiki/Window_function
    .. [4] F. J. Harris, "On the use of windows for harmonic analysis with the
           discrete Fourier transform," Proceedings of the IEEE, vol. 66,
           no. 1, pp. 51-83, Jan. 1978. :doi:`10.1109/PROC.1978.10837`.

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

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

    >>> window = signal.windows.kaiser(51, beta=14)
    >>> plt.plot(window)
    >>> plt.title(r"Kaiser window ($\beta$=14)")
    >>> 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(r"Frequency response of the Kaiser window ($\beta$=14)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    # Docstring adapted from NumPy's kaiser function
    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device)
    alpha = (M - 1) / 2.0
    w = (special.i0(beta * xp.sqrt(1 - ((n - alpha) / alpha) ** 2.0)) /
         special.i0(xp.asarray(beta, dtype=xp.float64)))

    return _truncate(w, needs_trunc)


@xp_capabilities()
def kaiser_bessel_derived(M, beta, *, sym=True, xp=None, device=None):
    """Return a Kaiser-Bessel derived window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
        Note that this window is only defined for an even
        number of points.
    beta : float
        Kaiser window shape parameter.
    sym : bool, optional
        This parameter only exists to comply with the interface offered by
        the other window functions and to be callable by `get_window`.
        When True (default), generates a symmetric window, for use in filter
        design.
    %(xp_device_snippet)s

    Returns
    -------
    w : ndarray
        The window, normalized to fulfil the Princen-Bradley condition.

    See Also
    --------
    kaiser

    Notes
    -----
    It is designed to be suitable for use with the modified discrete cosine
    transform (MDCT) and is mainly used in audio signal processing and
    audio coding.

    .. versionadded:: 1.9.0

    References
    ----------
    .. [1] Bosi, Marina, and Richard E. Goldberg. Introduction to Digital
           Audio Coding and Standards. Dordrecht: Kluwer, 2003.
    .. [2] Wikipedia, "Kaiser window",
           https://en.wikipedia.org/wiki/Kaiser_window

    Examples
    --------
    Plot the Kaiser-Bessel derived window based on the wikipedia
    reference [2]_:

    >>> import numpy as np
    >>> from scipy import signal
    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots()
    >>> N = 50
    >>> for alpha in [0.64, 2.55, 7.64, 31.83]:
    ...     ax.plot(signal.windows.kaiser_bessel_derived(2*N, np.pi*alpha),
    ...             label=f"{alpha=}")
    >>> ax.grid(True)
    >>> ax.set_title("Kaiser-Bessel derived window")
    >>> ax.set_ylabel("Amplitude")
    >>> ax.set_xlabel("Sample")
    >>> ax.set_xticks([0, N, 2*N-1])
    >>> ax.set_xticklabels(["0", "N", "2N+1"])  # doctest: +SKIP
    >>> ax.set_yticks([0.0, 0.2, 0.4, 0.6, 0.707, 0.8, 1.0])
    >>> fig.legend(loc="center")
    >>> fig.tight_layout()
    >>> fig.show()
    """
    xp = _namespace(xp)

    if not sym:
        raise ValueError(
            "Kaiser-Bessel Derived windows are only defined for symmetric "
            "shapes"
        )
    elif M < 1:
        return xp.asarray([])
    elif M % 2:
        raise ValueError(
            "Kaiser-Bessel Derived windows are only defined for even number "
            "of points"
        )

    kaiser_window = kaiser(M // 2 + 1, beta, xp=xp, device=device)
    csum = xp.cumulative_sum(kaiser_window)
    half_window = xp.sqrt(csum[:-1] / csum[-1])
    w = xp.concat((half_window, xp.flip(half_window)), axis=0)
    return xp.asarray(w, device=device)


@xp_capabilities()
def gaussian(M, std, sym=True, *, xp=None, device=None):
    r"""Return a Gaussian window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    std : float
        The standard deviation, sigma.
    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.
    %(xp_device_snippet)s

    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 Gaussian window is defined as

    .. math::  w(n) = e^{ -\frac{1}{2}\left(\frac{n}{\sigma}\right)^2 }

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

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

    >>> window = signal.windows.gaussian(51, std=7)
    >>> plt.plot(window)
    >>> plt.title(r"Gaussian window ($\sigma$=7)")
    >>> 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(r"Frequency response of the Gaussian window ($\sigma$=7)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device) - (M - 1.0) / 2.0
    sig2 = 2 * std * std
    w = xp.exp(-n ** 2 / sig2)

    return _truncate(w, needs_trunc)


@xp_capabilities()
def general_gaussian(M, p, sig, sym=True, *, xp=None, device=None):
    r"""Return a window with a generalized Gaussian shape.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    p : float
        Shape parameter.  p = 1 is identical to `gaussian`, p = 0.5 is
        the same shape as the Laplace distribution.
    sig : float
        The standard deviation, sigma.
    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.
    %(xp_device_snippet)s

    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 generalized Gaussian window is defined as

    .. math::  w(n) = e^{ -\frac{1}{2}\left|\frac{n}{\sigma}\right|^{2p} }

    the half-power point is at

    .. math::  (2 \log(2))^{1/(2 p)} \sigma

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

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

    >>> window = signal.windows.general_gaussian(51, p=1.5, sig=7)
    >>> plt.plot(window)
    >>> plt.title(r"Generalized Gaussian window (p=1.5, $\sigma$=7)")
    >>> 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(r"Freq. resp. of the gen. Gaussian "
    ...           r"window (p=1.5, $\sigma$=7)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    n = xp.arange(0, M, dtype=xp.float64, device=device) - (M - 1.0) / 2.0
    w = xp.exp(-0.5 * abs(n / sig) ** (2 * p))

    return _truncate(w, needs_trunc)


# `chebwin` contributed by Kumar Appaiah.
@xp_capabilities(skip_backends=(("jax.numpy", "item assignment"),
                                ("dask.array", "data-dependent output shapes")))
def chebwin(M, at, sym=True, *, xp=None, device=None):
    r"""Return a Dolph-Chebyshev window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

    >>> window = signal.windows.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]")

    """
    xp = _namespace(xp)

    if abs(at) < 45:
        warnings.warn("This window is not suitable for spectral analysis "
                      "for attenuation values lower than about 45dB because "
                      "the equivalent noise bandwidth of a Chebyshev window "
                      "does not grow monotonically with increasing sidelobe "
                      "attenuation when the attenuation is smaller than "
                      "about 45 dB.",
                      stacklevel=2)
    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    # compute the parameter beta
    order = M - 1.0
    _val = xp.asarray(10 ** (abs(at) / 20.), dtype=xp.float64, device=device)
    beta = xp.cosh(1.0 / order * xp.acosh(_val))
    k = xp.arange(M, dtype=xp.float64, device=device)
    x = beta * xp.cos(xp.pi * k / M)
    # Find the window's DFT coefficients
    # Use analytic definition of Chebyshev polynomial instead of expansion
    # from scipy.special. Using the expansion in scipy.special leads to errors.
    p = xp.zeros_like(x)
    p[x > 1] = xp.cosh(order * xp.acosh(x[x > 1]))
    p[x < -1] = (2 * (M % 2) - 1) * xp.cosh(order * xp.acosh(-x[x < -1]))
    p[abs(x) <= 1] = xp.cos(order * xp.acos(x[abs(x) <= 1]))

    # Appropriate IDFT and filling up
    # depending on even/odd M
    if M % 2:
        w = xp.real(sp_fft.fft(p))
        n = (M + 1) // 2
        w = w[:n]
        w = xp.concat((xp.flip(w[1:n]), w))
    else:
        p = p * xp.exp(1j * xp.pi / M * xp.arange(M, dtype=xp.float64, device=device))
        w = xp.real(sp_fft.fft(p))
        n = M // 2 + 1
        w = xp.concat((xp.flip(w[1:n]), w[1:n]))
    w = w / xp.max(w)

    return _truncate(w, needs_trunc)


@xp_capabilities()
def cosine(M, sym=True, *, xp=None, device=None):
    """Return a window with a simple cosine shape.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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:

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

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

    >>> plt.figure()
    >>> A = fft(window, 2047) / (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()

    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    w = xp.sin(xp.pi / M * (xp.arange(M, dtype=xp.float64, device=device) + .5))

    return _truncate(w, needs_trunc)


@xp_capabilities()
def exponential(M, center=None, tau=1., sym=True, *, xp=None, device=None):
    r"""Return an exponential (or Poisson) window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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
    ----------
    .. [1] 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:

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

    >>> M = 51
    >>> tau = 3.0
    >>> window = signal.windows.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.windows.exponential(M, 0, tau2, False)
    >>> plt.figure()
    >>> plt.plot(window2)
    >>> plt.ylabel("Amplitude")
    >>> plt.xlabel("Sample")
    """
    xp = _namespace(xp)

    if sym and center is not None:
        raise ValueError("If sym==True, center must be None.")
    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    if center is None:
        center = (M-1) / 2

    n = xp.arange(0, M, dtype=xp.float64, device=device)
    w = xp.exp(-abs(n-center) / tau)

    return _truncate(w, needs_trunc)


@xp_capabilities(skip_backends=[("jax.numpy", "item assignment")])
def taylor(M, nbar=4, sll=30, norm=True, sym=True, *, xp=None, device=None):
    """
    Return a Taylor window.

    The Taylor window taper function approximates the Dolph-Chebyshev window's
    constant sidelobe level for a parameterized number of near-in sidelobes,
    but then allows a taper beyond [2]_.

    The SAR (synthetic aperture radar) community commonly uses Taylor
    weighting for image formation processing because it provides strong,
    selectable sidelobe suppression with minimum broadening of the
    mainlobe [1]_.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    nbar : int, optional
        Number of nearly constant level sidelobes adjacent to the mainlobe.
    sll : float, optional
        Desired suppression of sidelobe level in decibels (dB) relative to the
        DC gain of the mainlobe. This should be a positive number.
    norm : bool, optional
        When True (default), divides the window by the largest (middle) value
        for odd-length windows or the value that would occur between the two
        repeated middle values for even-length windows such that all values
        are less than or equal to 1. When False the DC gain will remain at 1
        (0 dB) and the sidelobes will be `sll` dB down.
    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.
    %(xp_device_snippet)s

    Returns
    -------
    out : array
        The window. When `norm` is True (default), the maximum value is
        normalized to 1 (though the value 1 does not appear if `M` is
        even and `sym` is True).

    See Also
    --------
    chebwin, kaiser, bartlett, blackman, hamming, hann

    References
    ----------
    .. [1] W. Carrara, R. Goodman, and R. Majewski, "Spotlight Synthetic
           Aperture Radar: Signal Processing Algorithms" Pages 512-513,
           July 1995.
    .. [2] Armin Doerry, "Catalog of Window Taper Functions for
           Sidelobe Control", 2017.
           https://www.researchgate.net/profile/Armin_Doerry/publication/316281181_Catalog_of_Window_Taper_Functions_for_Sidelobe_Control/links/58f92cb2a6fdccb121c9d54d/Catalog-of-Window-Taper-Functions-for-Sidelobe-Control.pdf

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

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

    >>> window = signal.windows.taylor(51, nbar=20, sll=100, norm=False)
    >>> plt.plot(window)
    >>> plt.title("Taylor 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 Taylor window (100 dB)")
    >>> plt.ylabel("Normalized magnitude [dB]")
    >>> plt.xlabel("Normalized frequency [cycles per sample]")

    """  # noqa: E501
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    # Original text uses a negative sidelobe level parameter and then negates
    # it in the calculation of B. To keep consistent with other methods we
    # assume the sidelobe level parameter to be positive.
    B = xp.asarray(10**(sll / 20), device=device)
    A = xp.acosh(B) / xp.pi
    s2 = nbar**2 / (A**2 + (nbar - 0.5)**2)
    ma = xp.arange(1, nbar, dtype=xp.float64, device=device)

    Fm = xp.empty(nbar - 1, dtype=xp.float64, device=device)
    signs = xp.empty_like(ma)
    signs[::2] = 1
    signs[1::2] = -1
    m2 = ma*ma
    for mi, m in enumerate(ma):
        numer = signs[mi] * xp.prod(1 - m2[mi]/s2/(A**2 + (ma - 0.5)**2))
        denom = 2 * xp.prod(1 - m2[mi]/m2[:mi]) * xp.prod(1 - m2[mi]/m2[mi+1:])
        Fm[mi] = numer / denom

    def W(n):
        return 1 + 2*xp.matmul(Fm, xp.cos(
            2*xp.pi*ma[:, xp.newaxis]*(n-M/2.+0.5)/M))

    w = W(xp.arange(M, dtype=xp.float64, device=device))

    # normalize (Note that this is not described in the original text [1])
    if norm:
        scale = 1.0 / W((M - 1) / 2)
        w *= scale

    return _truncate(w, needs_trunc)


@xp_capabilities(np_only=True, reason='banded linear algebra is numpy-only')
def dpss(M, NW, Kmax=None, sym=True, norm=None, return_ratios=False,
         *, xp=None, device=None):
    """
    Compute the Discrete Prolate Spheroidal Sequences (DPSS).

    DPSS (or Slepian sequences) are often used in multitaper power spectral
    density estimation (see [1]_). The first window in the sequence can be
    used to maximize the energy concentration in the main lobe, and is also
    called the Slepian window.

    Parameters
    ----------
    M : int
        Window length.
    NW : float
        Standardized half bandwidth corresponding to ``2*NW = BW/f0 = BW*M*dt``
        where ``dt`` is taken as 1.
    Kmax : int | None, optional
        Number of DPSS windows to return (orders ``0`` through ``Kmax-1``).
        If None (default), return only a single window of shape ``(M,)``
        instead of an array of windows of shape ``(Kmax, M)``.
    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.
    norm : {2, 'approximate', 'subsample'} | None, optional
        If 'approximate' or 'subsample', then the windows are normalized by the
        maximum, and a correction scale-factor for even-length windows
        is applied either using ``M**2/(M**2+NW)`` ("approximate") or
        a FFT-based subsample shift ("subsample"), see Notes for details.
        If None, then "approximate" is used when ``Kmax=None`` and 2 otherwise
        (which uses the l2 norm).
    return_ratios : bool, optional
        If True, also return the concentration ratios in addition to the
        windows.
    %(xp_device_snippet)s

    Returns
    -------
    v : ndarray, shape (Kmax, M) or (M,)
        The DPSS windows. Will be 1D if `Kmax` is None.
    r : ndarray, shape (Kmax,) or float, optional
        The concentration ratios for the windows. Only returned if
        `return_ratios` evaluates to True. Will be 0D if `Kmax` is None.

    Notes
    -----
    This computation uses the tridiagonal eigenvector formulation given
    in [2]_.

    The default normalization for ``Kmax=None``, i.e. window-generation mode,
    simply using the l-infinity norm would create a window with two unity
    values, which creates slight normalization differences between even and odd
    orders. The approximate correction of ``M**2/float(M**2+NW)`` for even
    sample numbers is used to counteract this effect (see Examples below).

    For very long signals (e.g., 1e6 elements), it can be useful to compute
    windows orders of magnitude shorter and use interpolation (e.g.,
    `scipy.interpolate.interp1d`) to obtain tapers of length `M`,
    but this in general will not preserve orthogonality between the tapers.

    .. versionadded:: 1.1

    References
    ----------
    .. [1] Percival DB, Walden WT. Spectral Analysis for Physical Applications:
       Multitaper and Conventional Univariate Techniques.
       Cambridge University Press; 1993.
    .. [2] Slepian, D. Prolate spheroidal wave functions, Fourier analysis, and
       uncertainty V: The discrete case. Bell System Technical Journal,
       Volume 57 (1978), 1371430.
    .. [3] Kaiser, JF, Schafer RW. On the Use of the I0-Sinh Window for
       Spectrum Analysis. IEEE Transactions on Acoustics, Speech and
       Signal Processing. ASSP-28 (1): 105-107; 1980.

    Examples
    --------
    We can compare the window to `kaiser`, which was invented as an alternative
    that was easier to calculate [3]_ (example adapted from
    `here <https://ccrma.stanford.edu/~jos/sasp/Kaiser_DPSS_Windows_Compared.html>`_):

    >>> import numpy as np
    >>> import matplotlib.pyplot as plt
    >>> from scipy.signal import windows, freqz
    >>> M = 51
    >>> fig, axes = plt.subplots(3, 2, figsize=(5, 7))
    >>> for ai, alpha in enumerate((1, 3, 5)):
    ...     win_dpss = windows.dpss(M, alpha)
    ...     beta = alpha*np.pi
    ...     win_kaiser = windows.kaiser(M, beta)
    ...     for win, c in ((win_dpss, 'k'), (win_kaiser, 'r')):
    ...         win /= win.sum()
    ...         axes[ai, 0].plot(win, color=c, lw=1.)
    ...         axes[ai, 0].set(xlim=[0, M-1], title=rf'$\\alpha$ = {alpha}',
    ...                         ylabel='Amplitude')
    ...         w, h = freqz(win)
    ...         axes[ai, 1].plot(w, 20 * np.log10(np.abs(h)), color=c, lw=1.)
    ...         axes[ai, 1].set(xlim=[0, np.pi],
    ...                         title=rf'$\\beta$ = {beta:0.2f}',
    ...                         ylabel='Magnitude (dB)')
    >>> for ax in axes.ravel():
    ...     ax.grid(True)
    >>> axes[2, 1].legend(['DPSS', 'Kaiser'])
    >>> fig.tight_layout()
    >>> plt.show()

    And here are examples of the first four windows, along with their
    concentration ratios:

    >>> M = 512
    >>> NW = 2.5
    >>> win, eigvals = windows.dpss(M, NW, 4, return_ratios=True)
    >>> fig, ax = plt.subplots(1)
    >>> ax.plot(win.T, linewidth=1.)
    >>> ax.set(xlim=[0, M-1], ylim=[-0.1, 0.1], xlabel='Samples',
    ...        title=f'DPSS, {M:d}, {NW:0.1f}')
    >>> ax.legend([f'win[{ii}] ({ratio:0.4f})'
    ...            for ii, ratio in enumerate(eigvals)])
    >>> fig.tight_layout()
    >>> plt.show()

    Using a standard :math:`l_{\\infty}` norm would produce two unity values
    for even `M`, but only one unity value for odd `M`. This produces uneven
    window power that can be counteracted by the approximate correction
    ``M**2/float(M**2+NW)``, which can be selected by using
    ``norm='approximate'`` (which is the same as ``norm=None`` when
    ``Kmax=None``, as is the case here). Alternatively, the slower
    ``norm='subsample'`` can be used, which uses subsample shifting in the
    frequency domain (FFT) to compute the correction:

    >>> Ms = np.arange(1, 41)
    >>> factors = (50, 20, 10, 5, 2.0001)
    >>> energy = np.empty((3, len(Ms), len(factors)))
    >>> for mi, M in enumerate(Ms):
    ...     for fi, factor in enumerate(factors):
    ...         NW = M / float(factor)
    ...         # Corrected using empirical approximation (default)
    ...         win = windows.dpss(M, NW)
    ...         energy[0, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
    ...         # Corrected using subsample shifting
    ...         win = windows.dpss(M, NW, norm='subsample')
    ...         energy[1, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
    ...         # Uncorrected (using l-infinity norm)
    ...         win /= win.max()
    ...         energy[2, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
    >>> fig, ax = plt.subplots(1)
    >>> hs = ax.plot(Ms, energy[2], '-o', markersize=4,
    ...              markeredgecolor='none')
    >>> leg = [hs[-1]]
    >>> for hi, hh in enumerate(hs):
    ...     h1 = ax.plot(Ms, energy[0, :, hi], '-o', markersize=4,
    ...                  color=hh.get_color(), markeredgecolor='none',
    ...                  alpha=0.66)
    ...     h2 = ax.plot(Ms, energy[1, :, hi], '-o', markersize=4,
    ...                  color=hh.get_color(), markeredgecolor='none',
    ...                  alpha=0.33)
    ...     if hi == len(hs) - 1:
    ...         leg.insert(0, h1[0])
    ...         leg.insert(0, h2[0])
    >>> ax.set(xlabel='M (samples)', ylabel=r'Power / $\\sqrt{M}$')
    >>> ax.legend(leg, ['Uncorrected', r'Corrected: $\\frac{M^2}{M^2+NW}$',
    ...                 'Corrected (subsample)'])
    >>> fig.tight_layout()

    """
    xp = _namespace(xp)

    if norm is None:
        norm = 'approximate' if Kmax is None else 2
    known_norms = (2, 'approximate', 'subsample')
    if norm not in known_norms:
        raise ValueError(f'norm must be one of {known_norms}, got {norm}')
    if Kmax is None:
        singleton = True
        Kmax = 1
    else:
        singleton = False
    if _len_guards(M):
        if not return_ratios:
            return xp.ones(M, dtype=xp.float64)
        elif singleton:
            return xp.ones(M, dtype=xp.float64), 1.
        else:
            return xp.ones(M, dtype=xp.float64), xp.ones(1, dtype=xp.float64)
    Kmax = operator.index(Kmax)
    if not 0 < Kmax <= M:
        raise ValueError('Kmax must be greater than 0 and less than M')
    if NW >= M/2.:
        raise ValueError('NW must be less than M/2.')
    if NW <= 0:
        raise ValueError('NW must be positive')
    M, needs_trunc = _extend(M, sym)
    W = float(NW) / M
    nidx = xp.arange(M, dtype=xp.float64, device=device)

    # Here we want to set up an optimization problem to find a sequence
    # whose energy is maximally concentrated within band [-W,W].
    # Thus, the measure lambda(T,W) is the ratio between the energy within
    # that band, and the total energy. This leads to the eigen-system
    # (A - (l1)I)v = 0, where the eigenvector corresponding to the largest
    # eigenvalue is the sequence with maximally concentrated energy. The
    # collection of eigenvectors of this system are called Slepian
    # sequences, or discrete prolate spheroidal sequences (DPSS). Only the
    # first K, K = 2NW/dt orders of DPSS will exhibit good spectral
    # concentration
    # [see https://en.wikipedia.org/wiki/Spectral_concentration_problem]

    # Here we set up an alternative symmetric tri-diagonal eigenvalue
    # problem such that
    # (B - (l2)I)v = 0, and v are our DPSS (but eigenvalues l2 != l1)
    # the main diagonal = ([M-1-2*t]/2)**2 cos(2PIW), t=[0,1,2,...,M-1]
    # and the first off-diagonal = t(M-t)/2, t=[1,2,...,M-1]
    # [see Percival and Walden, 1993]
    d = ((M - 1 - 2 * nidx) / 2.) ** 2 * xp.cos(xp.asarray(2 * xp.pi * W))
    e = nidx[1:] * (M - nidx[1:]) / 2.

    # only calculate the highest Kmax eigenvalues
    w, windows = linalg.eigh_tridiagonal(
        d, e, select='i', select_range=(M - Kmax, M - 1))
    w = w[::-1]
    windows = windows[:, ::-1].T

    # By convention (Percival and Walden, 1993 pg 379)
    # * symmetric tapers (k=0,2,4,...) should have a positive average.
    fix_even = (windows[::2, ...].sum(axis=1) < 0)
    for i, f in enumerate(fix_even):
        if f:
            windows[2 * i] *= -1
    # * antisymmetric tapers should begin with a positive lobe
    #   (this depends on the definition of "lobe", here we'll take the first
    #   point above the numerical noise, which should be good enough for
    #   sufficiently smooth functions, and more robust than relying on an
    #   algorithm that uses max(abs(w)), which is susceptible to numerical
    #   noise problems)
    thresh = max(1e-7, 1. / M)
    for i, w in enumerate(windows[1::2, ...]):
        if w[w * w > thresh][0] < 0:
            windows[2 * i + 1] *= -1

    # Now find the eigenvalues of the original spectral concentration problem
    # Use the autocorr sequence technique from Percival and Walden, 1993 pg 390
    if return_ratios:
        dpss_rxx = _fftautocorr(xp.asarray(windows))
        r = 4 * W * xpx.sinc(xp.asarray(2 * W * nidx), xp=xp)
        r[0] = 2 * W
        ratios = xp.matmul(dpss_rxx, r)
        if singleton:
            ratios = ratios[0]
        ratios = xp.asarray(ratios, device=device)
    # Deal with sym and Kmax=None
    if norm != 2:
        windows /= windows.max()
        if M % 2 == 0:
            if norm == 'approximate':
                correction = M**2 / float(M**2 + NW)
            else:
                s = sp_fft.rfft(windows[0])
                shift = -(1 - 1./M) * xp.arange(1, M//2 + 1, dtype=xp.float64)
                s[1:] *= 2 * xp.exp(-1j * xp.pi * shift)
                correction = M / s.real.sum()
            windows *= correction
    # else we're already l2 normed, so do nothing
    if needs_trunc:
        windows = windows[:, :-1]
    if singleton:
        windows = windows[0]
    windows = xp.asarray(windows, device=device)
    return (windows, ratios) if return_ratios else windows


@xp_capabilities()
def lanczos(M, *, sym=True, xp=None, device=None):
    r"""Return a Lanczos window also known as a sinc window.

    Parameters
    ----------
    M : int
        Number of points in the output window. If zero, an empty array
        is returned. An exception is thrown when it is negative.
    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.
    %(xp_device_snippet)s

    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 Lanczos window is defined as

    .. math::  w(n) = sinc \left( \frac{2n}{M - 1} - 1 \right)

    where

    .. math::  sinc(x) = \frac{\sin(\pi x)}{\pi x}

    The Lanczos window has reduced Gibbs oscillations and is widely used for
    filtering climate timeseries with good properties in the physical and
    spectral domains.

    .. versionadded:: 1.10

    References
    ----------
    .. [1] Lanczos, C., and Teichmann, T. (1957). Applied analysis.
           Physics Today, 10, 44.
    .. [2] Duchon C. E. (1979) Lanczos Filtering in One and Two Dimensions.
           Journal of Applied Meteorology, Vol 18, pp 1016-1022.
    .. [3] Thomson, R. E. and Emery, W. J. (2014) Data Analysis Methods in
           Physical Oceanography (Third Edition), Elsevier, pp 593-637.
    .. [4] Wikipedia, "Window function",
           http://en.wikipedia.org/wiki/Window_function

    Examples
    --------
    Plot the window

    >>> import numpy as np
    >>> from scipy.signal.windows import lanczos
    >>> from scipy.fft import fft, fftshift
    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(1)
    >>> window = lanczos(51)
    >>> ax.plot(window)
    >>> ax.set_title("Lanczos window")
    >>> ax.set_ylabel("Amplitude")
    >>> ax.set_xlabel("Sample")
    >>> fig.tight_layout()
    >>> plt.show()

    and its frequency response:

    >>> fig, ax = plt.subplots(1)
    >>> 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())))
    >>> ax.plot(freq, response)
    >>> ax.set_xlim(-0.5, 0.5)
    >>> ax.set_ylim(-120, 0)
    >>> ax.set_title("Frequency response of the lanczos window")
    >>> ax.set_ylabel("Normalized magnitude [dB]")
    >>> ax.set_xlabel("Normalized frequency [cycles per sample]")
    >>> fig.tight_layout()
    >>> plt.show()
    """
    xp = _namespace(xp)

    if _len_guards(M):
        return xp.ones(M, dtype=xp.float64, device=device)
    M, needs_trunc = _extend(M, sym)

    # To make sure that the window is symmetric, we concatenate the right hand
    # half of the window and the flipped one which is the left hand half of
    # the window.
    def _calc_right_side_lanczos(n, m):
        return xpx.sinc(2. * xp.arange(n, m, dtype=xp.float64) / (m - 1) - 1.0, xp=xp)

    if M % 2 == 0:
        wh = _calc_right_side_lanczos(M/2, M)
        w = xp.concat([xp.flip(wh), wh])
    else:
        wh = _calc_right_side_lanczos((M+1)/2, M)
        w = xp.concat([xp.flip(wh), xp.ones(1), wh])

    return _truncate(w, needs_trunc)


def _fftautocorr(x):
    """Compute the autocorrelation of a real array and crop the result."""
    N = x.shape[-1]
    use_N = sp_fft.next_fast_len(2*N-1)
    x_fft = sp_fft.rfft(x, use_N, axis=-1)
    cxy = sp_fft.irfft(x_fft * x_fft.conj(), n=use_N)[:, :N]
    # Or equivalently (but in most cases slower):
    # cxy = xp.asarray([xp.convolve(xx, yy[::-1], mode='full')
    #                   for xx, yy in zip(x, x)])[:, N-1:2*N-1]
    return cxy

_WIN_FUNC_DATA = { # Format: {(name0, name1, ...): (function, needs parameters)
    ('barthann', 'brthan', 'bth'): (barthann, False),
    ('bartlett', 'bart', 'brt'): (bartlett, False),
    ('blackman', 'black', 'blk'): (blackman, False),
    ('blackmanharris', 'blackharr', 'bkh'): (blackmanharris, False),
    ('bohman', 'bman', 'bmn'): (bohman, False),
    ('boxcar', 'box', 'ones', 'rect', 'rectangular'): (boxcar, False),
    ('chebwin', 'cheb'): (chebwin, True),
    ('cosine', 'halfcosine'): (cosine, False),
    ('dpss',): (dpss, True),
    ('exponential', 'poisson'): (exponential, 'OPTIONAL'),
    ('flattop', 'flat', 'flt'): (flattop, False),
    ('gaussian', 'gauss', 'gss'): (gaussian, True),
    ('general cosine', 'general_cosine'): (general_cosine, True),
    ('general gaussian', 'general_gaussian', 'general gauss',
     'general_gauss', 'ggs'): (general_gaussian, True),
    ('general hamming', 'general_hamming'): (general_hamming, True),
    ('hamming', 'hamm', 'ham'): (hamming, False),
    ('hann', 'han'): (hann, False),
    ('kaiser', 'ksr'): (kaiser, True),
    ('kaiser bessel derived', 'kaiser_bessel_derived',
     'kbd'): (kaiser_bessel_derived, True),
    ('lanczos', 'sinc'): (lanczos, False),
    ('nuttall', 'nutl', 'nut'): (nuttall, False),
    ('parzen', 'parz', 'par'): (parzen, False),
    ('taylor', 'taylorwin'): (taylor, 'OPTIONAL'),
    ('triangle', 'triang', 'tri'): (triang, False),
    ('tukey', 'tuk'): (tukey, 'OPTIONAL'), }
_WIN_FUNCS = dict()
for nn_, v_ in _WIN_FUNC_DATA.items():
    _WIN_FUNCS.update({n_: v_ for n_ in nn_})


@xp_capabilities()
def get_window(window, Nx, fftbins=True, *, xp=None, device=None):
    r"""Convenience function for creating various windows.

    This function is a wrapper for the window functions provided in the
    `scipy.signal.windows` namespace.

    Parameters
    ----------
    window : str | tuple | float
        Either a string with the window name or a tuple consisting of window name and
        window parameters. If it is a float, a `~scipy.singal.kaiser` window is
        created with `window` being the shape parameter. Consult the Notes below for
        more details.
    Nx : int
        The number of samples in the window.
    fftbins : bool, optional
        If ``True`` (default), create a periodic window, ready to use with `ifftshift`
        and be multiplied by the result of an FFT (see also
        :func:`~scipy.fft.fftfreq`). If ``False``, create a symmetric window, for use
        in filter design. This parameter is ignored, if the window name in the
        `window` parameter has a suffix ``'_periodic'`` or ``'_symmetric'`` appended to
        it (e.g., ``'hann_symmetric'``).
    %(xp_device_snippet)s

    Returns
    -------
    get_window : ndarray
        Returns the created window as a one-dimensional array made of `Nx` samples.

    Raises
    ------
    ValueError
        If the provided parameters do not allow to choose a valid window function
        with valid parameters.

    Notes
    -----
    Note that by default this function returns a periodic window, whereas the wrapped
    window functions return a symmetric window by default. This is caused by the
    `fftbins` parameter having the inverse meaning of the `sym` parameter of the
    wrapped function, which are both ``True`` by default.

    .. currentmodule:: scipy.signal.windows

    The following list shows the wrapped window functions with the respective settings
    of the `window` parameter followed by a short description. Aliases for alternative
    window names are given in parenthesis.

    `barthann` / ``'barthann'``:
        Modified Bartlett-Hann window (aliases: ``'brthan', 'bth'``)
    `bartlett` / ``'bartlett'``:
        Bartlett window (aliases: ``'bart', 'brt'``)
    `blackman`/ ``'blackman'``:
        Blackman window (aliases: ``'black', 'blk'``)
    `blackmanharris` / ``'blackmanharris'``:
        4-term Blackman-Harris window (aliases: ``'blackharr', 'bkh'``)
    `bohman` / ``'bohman'``:
        Bohman window (aliases: ``'bman', 'bmn'``)
    `boxcar` / ``'boxcar'``:
        Rectangular window (aliases: ``'box', 'ones', 'rect', 'rectangular'``)
    `chebwin` / ``('chebwin', at)``:
        Dolph-Chebyshev window with `at` dB attenuation (aliases: ``'cheb'``)
    `cosine` / ``'cosine'``:
        Cosine window (aliases: ``'halfcosine'``)
    `dpss` / ``('dpss', NW)``:
        First window of discrete prolate spheroidal sequence with standardized half
        bandwidth ``NW`` and "approximate" norm.
    `exponential` / ``'exponential'`` / ``('exponential', center, tau)``:
        Exponential / Poisson window centered at ``center`` (default: ``None``) with
        deacy ``tau`` (default: ``1``) (aliases: ``'poisson'``)
    `flattop`/ ``'flattop'``:
        Flat top window (aliases: ``'flat', 'flt'``)
    `gaussian` / ``('gaussian', std)``:
        Gaussian with standard deviation `std` (aliases: ``'gauss', 'gss'``)
    `general_cosine` / ``('general cosine', a)``:
        Generic weighted sum of cosine terms with weighting coefficients ``a``
        (aliases: ``'general_cosine'``)
    `general_gaussian` / ``('general gaussian', p, sig)``:
        Generalized Gaussian with shape parameter ``p`` and standard deviation ``sig``
        (aliases: ``'general_gaussian', 'general gauss', 'general_gauss', 'ggs'``)
    `general_hamming` / ``('general hamming', alpha)``:
        Generalized Hamming window with coefficent ``alpha``
        (aliases: ``'general_hamming'``)
    `hamming` / ``'hamming'``:
        Hamming window (aliases: ``'hamm', 'ham'``)
    `hann` / ``'hann'``:
        Hann window (aliases: ``'han'``)
    `kaiser` / ``('kaiser', beta)``:
        Kaiser window with shape parameter ``beta`` (aliases: ``'ksr'``)
    `kaiser_bessel_derived` / ``('kaiser bessel derived', beta)``:
        Kaiser-Bessel derived window with shape parameter ``beta``
        (aliases: ``'kaiser_bessel_derived', 'kbd'``)
    `lanczos` / ``'lanczos'``:
        Lanczos / sinc window (aliases: ``'sinc'``)
    `nuttall`/ ``'nuttall'``:
        Minimum 4-term Blackman-Harris window according to Nuttall
        (aliases: ``'nutl', 'nut'``)
    `parzen` / ``'parzen'``:
        Parzen window (aliases: ``'parz', 'par'``)
    `taylor` / ``'taylor'`` / (``'taylor', nbar, sll, norm)``:
        Taylor window with ``nbar`` adjascent sidelobes (default: ``4``)), ``sll`` dB
        suppression level (default: ``30``) and boolean value ``norm``
        (default: ``True``) (aliases: ``taylorwin``)
    `triang`/ ``'triangle'``:
        Triangle window (aliases: ``'triang', 'tri'``)
    `tukey` / ``'tukey'`` / ``('tukey', alpha)``:
        Tukey window with shape parameter ``alpha`` (default: ``0.5``)
        (aliases: ``'tuk'``)


    Examples
    --------
    This example shows different usages of the `window` parameter:

    >>> from scipy.signal import get_window
    >>> get_window('triang', 7)
    array([ 0.125,  0.375,  0.625,  0.875,  0.875,  0.625,  0.375])
    >>> get_window(('exponential', None, 1.), 9)
    array([ 0.011109  ,  0.03019738,  0.082085  ,  0.22313016,  0.60653066,
            0.60653066,  0.22313016,  0.082085  ,  0.03019738])
    >>> get_window(('kaiser', 4.0), 9)
    array([ 0.08848053,  0.29425961,  0.56437221,  0.82160913,  0.97885093,
            0.97885093,  0.82160913,  0.56437221,  0.29425961])
    >>> get_window(4.0, 9)  # same as previous call
    array([ 0.08848053,  0.29425961,  0.56437221,  0.82160913,  0.97885093,
            0.97885093,  0.82160913,  0.56437221,  0.29425961])

    The following snippet shows different ways to create identical symmetric and
    periodic Bartlett windows:

    >>> from scipy.signal import get_window, windows
    >>> # Symmetric window:
    >>> windows.bartlett(5)  # Parameter `sym` defaults to True
    array([0. , 0.5, 1. , 0.5, 0. ])
    >>> get_window('bartlett', 5, fftbins=False)
    array([0. , 0.5, 1. , 0.5, 0. ])
    >>> get_window('bartlett_symmetric', 5)
    array([0. , 0.5, 1. , 0.5, 0. ])
    >>> # Periodic window:
    >>> windows.bartlett(4, sym=False)
    array([0. , 0.5, 1. , 0.5])
    >>> get_window('bartlett', 4)  # Parameter `fftbins` defaults to True
    array([0. , 0.5, 1. , 0.5])
    >>> get_window('bartlett_periodic', 4)
    array([0. , 0.5, 1. , 0.5])
    >>> # `_periodic' suffix overrides `fftbins` parameter:
    >>> get_window('bartlett_periodic', 4, fftbins=False)
    array([0. , 0.5, 1. , 0.5])

    Note that a periodic window can be created out of a symmetric window by discarding
    the last sample.
    """
    if not (Nx > 0 and isinstance(Nx, numbers.Integral)):
        raise ValueError(f"Parameter {Nx=} is not a positive integer")
    if not isinstance(fftbins, bool):
        raise ValueError(f"Parameter {fftbins=} is not of type bool!")

    if not isinstance(window, str | tuple):
        try: # if parameter window can be converted to a float, return kaiser window:
            beta = float(window)
        except Exception as float_exception:
            err_msg = f"Parameter {window=} must be a tuple, a string or a float!"
            raise ValueError(err_msg) from float_exception
        return kaiser(Nx, beta, not fftbins, xp=xp, device=device)

    if isinstance(window, tuple) and not isinstance(window[0], str):
        raise ValueError(f"First tuple entry of parameter {window=} is not a str!")

    sym = not fftbins
    win_name = window if isinstance(window, str) else window[0]
    if win_name.endswith('_symmetric'):  # overwrite `fftbins` / `sym` if needed
        sym, win_name = True, win_name[:-10]  # remove '_symmetric' from `win_name`
    elif win_name.endswith('_periodic'):
        sym, win_name = False, win_name[:-9]  # remove '_periodic' from `win_name`

    if win_name not in _WIN_FUNCS:
        raise ValueError(f"Invalid window name '{win_name}' in parameter {window=}!")

    func, has_args = _WIN_FUNCS[win_name]
    args = window[1:] if isinstance(window, tuple) else tuple()
    if len(args) > 0 and has_args is False:
        raise ValueError(f"'{win_name}' does not allow parameters, but {window=}!")
    if len(args) == 0 and has_args is True:
        raise ValueError(f"'{win_name}' must have parameters, but {window=}!")
    # has_args == 'OPTIONAL' allows len(args) == 0 as well as len(args) > 0

    if not has_args:
        return func(Nx, sym=sym, xp=xp, device=device)

    # special cases taken from original implementation:
    if func is dpss:
        if len(args) != 1:
            raise ValueError(f"Window {win_name} must have one parameter but {window=}")
        return dpss(Nx, args[0], Kmax=None, sym=sym, xp=xp, device=device)
    if func is general_cosine:
        if not (xp is None and device is None):
            raise ValueError("'general_cosine' does not accept the parameters xp " +
                             "and device not being None!")
        return general_cosine(Nx, *args, sym=sym)

    return func(Nx, *args, sym=sym, xp=xp, device=device)


########## complete the docstrings, on import
_xp_device_snippet = {'xp_device_snippet':
"""\
xp : array_namespace, optional
    Optional array namespace.
    Should be compatible with the array API standard, or supported by array-api-compat.
    Default: ``numpy``
device: any
    optional device specification for output. Should match one of the
    supported device specification in ``xp``.
"""
}


_names = [x for x in __all__ if x != 'general_cosine']
for name in _names:
    window = vars()[name]
    window.__doc__ = doccer.docformat(window.__doc__, _xp_device_snippet)