Modified Discrete Cosine Transform (MDCT)

    Returns the Modified Discrete Cosine Transform with fixed
    window size L of the signal x.

    The window is based on a sine window.

    Parameters
    ----------
    x : ndarray, shape (N,)
        The signal
    L : int
        The window length

    Returns
    -------
    y : ndarray, shape (L/2, 2 * N / L)
        The MDCT coefficients

    See also
    --------
    imdct

def mdct(x, L):
    """Modified Discrete Cosine Transform (MDCT)

    Returns the Modified Discrete Cosine Transform with fixed
    window size L of the signal x.

    The window is based on a sine window.

    Parameters
    ----------
    x : ndarray, shape (N,)
        The signal
    L : int
        The window length

    Returns
    -------
    y : ndarray, shape (L/2, 2 * N / L)
        The MDCT coefficients

    See also
    --------
    imdct
    """
    x = np.asarray(x, dtype=np.float)

    N = x.size

    # Number of frequency channels
    K = L // 2

    # Test length
    if N % K != 0:
        raise RuntimeError('Input length must be a multiple of the half of '
                           'the window size')

    # Pad edges with zeros
    xx = np.zeros(L // 4 + N + L // 4)
    xx[L // 4:-L // 4] = x
    x = xx
    del xx

    # Number of frames
    P = N // K
    if P < 2:
        raise ValueError('Signal too short')

    # Framing
    x = _framing(x, L)

    # Windowing
    aL = np.arange(L, dtype=np.float)
    w_long = np.sin((np.pi / L) * (aL + 0.5))
    w_edge_L = w_long.copy()
    w_edge_L[:L // 4] = 0.
    w_edge_L[L // 4:L // 2] = 1.
    w_edge_R = w_long.copy()
    w_edge_R[L // 2:L // 2 + L // 4] = 1.
    w_edge_R[L // 2 + L // 4:] = 0.
    x[:, 0] *= w_edge_L
    x[:, 1:-1] *= w_long[:, None]
    x[:, -1] *= w_edge_R

    # Pre-twiddle
    x = x.astype(np.complex)
    x *= np.exp((-1j * np.pi / L) * aL)[:, None]

    # FFT
    y = fft(x, axis=0)

    # Post-twiddle
    y = y[:L // 2, :]
    y *= np.exp((-1j * np.pi * (L // 2 + 1.) / L)
                * (0.5 + aL[:L // 2]))[:, None]

    # Real part and scaling
    y = math.sqrt(2. / K) * np.real(y)
    return y