Repository: scipy/scipy
File: scipy/signal/_waveforms.py (Python)

Code before:
    --------
    A 5 Hz waveform sampled at 500 Hz for 1 second:

    >>> import numpy as np
    >>> from scipy import signal
    >>> import matplotlib.pyplot as plt
    >>> t = np.linspace(0, 1, 500, endpoint=False)
    >>> plt.plot(t, signal.square(2 * np.pi * 5 * t))
    >>> plt.ylim(-2, 2)

    A pulse-width modulated sine wave:

    >>> plt.figure()
    >>> sig = np.sin(2 * np.pi * t)
    >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
    >>> plt.subplot(2, 1, 1)
    >>> plt.plot(t, sig)
    >>> plt.subplot(2, 1, 2)
    >>> plt.plot(t, pwm)
    >>> plt.ylim(-1.5, 1.5)

    """
    xp = array_namespace(t, duty)
    t, w = xp_promote(t, duty, xp=xp)

    y = xp.zeros(t.shape, dtype=xp.float64)

    # width must be between 0 and 1 inclusive
    mask1 = (w > 1) | (w < 0)
    y = xpx.at(y, mask1).set(xp.nan)

    # on the interval 0 to duty*2*pi function is 1
    tmod = t % (2 * xp.pi)
    mask2 = ~mask1 & (tmod < w * 2 * pi)
    y = xpx.at(y, mask2).set(1)

    # on the interval duty*2*pi to 2*pi function is
    #  (pi*(w+1)-tmod) / (pi*(1-w))
    mask3 = ~mask1 & ~mask2
    y = xpx.at(y, mask3).set(-1)
    return y


def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
               retenv=False):
    """
    Return a Gaussian modulated sinusoid:

        ``exp(-a t^2) exp(1j*2*pi*fc*t).``

    If `retquad` is True, then return the real and imaginary parts

Review comment (suggestion):
Reviewing the history, I see that the original intention was to respect dtypes (gh-22252), but it mistakenly always returned `float64`, and there was an explicit decision to always return `float64` (gh-22520). I don't agree with it, but I'm not going to shake things up at this time.

Code after:
    --------
    A 5 Hz waveform sampled at 500 Hz for 1 second:

    >>> import numpy as np
    >>> from scipy import signal
    >>> import matplotlib.pyplot as plt
    >>> t = np.linspace(0, 1, 500, endpoint=False)
    >>> plt.plot(t, signal.square(2 * np.pi * 5 * t))
    >>> plt.ylim(-2, 2)

    A pulse-width modulated sine wave:

    >>> plt.figure()
    >>> sig = np.sin(2 * np.pi * t)
    >>> pwm = signal.square(2 * np.pi * 30 * t, duty=(sig + 1)/2)
    >>> plt.subplot(2, 1, 1)
    >>> plt.plot(t, sig)
    >>> plt.subplot(2, 1, 2)
    >>> plt.plot(t, pwm)
    >>> plt.ylim(-1.5, 1.5)

    """
    xp = array_namespace(t, duty)
    t, w = xp_promote(t, duty, xp=xp, force_floating=True, broadcast=True)

    y = xp.zeros(t.shape, dtype=t.dtype)

    # width must be between 0 and 1 inclusive
    mask1 = (w > 1) | (w < 0)
    y = xpx.at(y, mask1).set(xp.nan)

    # on the interval 0 to duty*2*pi function is 1
    tmod = t % (2 * xp.pi)
    mask2 = ~mask1 & (tmod < w*2*xp.pi)
    y = xpx.at(y, mask2).set(1)

    # on the interval duty*2*pi to 2*pi function is -1
    mask3 = ~(mask1 | mask2)
    y = xpx.at(y, mask3).set(-1)
    return y


def gausspulse(t, fc=1000, bw=0.5, bwr=-6, tpr=-60, retquad=False,
               retenv=False):
    """
    Return a Gaussian modulated sinusoid:

        ``exp(-a t^2) exp(1j*2*pi*fc*t).``

    If `retquad` is True, then return the real and imaginary parts
    (in-phase and quadrature).