only the first
    4 types are implemented in SciPy.

    **Type I**

    There are several definitions of the DST-I; we use the following for
    ``norm="backward"``. DST-I assumes the input is odd around :math:`n=-1` and
    :math:`n=N`.

    .. math::

        y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(n+1)}{N+1}\right)

    Note that the DST-I is only supported for input size > 1.
    The (unnormalized) DST-I is its own inverse, up to a factor :math:`2(N+1)`.
    The orthonormalized DST-I is exactly its own inverse.

    ``orthogonalize`` has no effect here, as the DST-I matrix is already
    orthogonal up to a scale factor of ``2N``.

    **Type II**

    There are several definitions of the DST-II; we use the following for
    ``norm="backward"``. DST-II assumes the input is odd around :math:`n=-1/2` and
    :math:`n=N-1/2`; the output is odd around :math:`k=-1` and even around :math:`k=N-1`

    .. math::

        y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(k+1)(2n+1)}{2N}\right)

    If ``orthogonalize=True``, ``y[0]`` is divided :math:`\sqrt{2}` which, when
    combined with ``norm="ortho"``, makes the corresponding matrix of
    coefficients orthonormal (``O @ O.T = np.eye(N)``).

    **Type III**

    There are several definitions of the DST-III, we use the following (for
    ``norm="backward"``). DST-III assumes the input is odd around :math:`n=-1` and
    even around :math:`n=N-1`

    .. math::

        y_k = (-1)^k x_{N-1} + 2 \sum_{n=0}^{N-2} x_n \sin\left(
        \frac{\pi(2k+1)(n+1)}{2N}\right)

    If ``orthogonalize=True``, ``x[0]`` is multiplied by :math:`\sqrt{2}`
    which, when combined with ``norm="ortho"``, makes the corresponding matrix
    of coefficients orthonormal (``O @ O.T = np.eye(N)``).

    The (unnormalized) DST-III is the inverse of the (unnormalized) DST-II, up
    to a factor :math:`2N`. The orthonormalized DST-III is exactly the inverse of the
    orthonormalized DST-II.

    **Type IV**

    There are several definitions of the DST-IV, we use the following (for
    ``norm="backward"``). DST-IV assumes the input is odd around :math:`n=-0.5` and
    even around :math:`n=N-0.5`

    .. math::

        y_k = 2 \sum_{n=0}^{N-1} x_n \sin\left(\frac{\pi(2k+1)(2n+1)}{4N}\right)

    ``orthogonalize`` has no effect here, as the DST-IV matrix is already
    orthogonal up to a scale factor of ``2N``.

    The (unnormalized) DST-IV is its own inverse, up to a factor :math:`2N`. The
    orthonormalized DST-IV is exactly its own inverse.

    References
    ----------
    .. [1] Wikipedia, "Discrete sine transform",
           https://en.wikipedia.org/wiki/Discrete_sine_transform

    r