def METHOD_NAME(N, w, top_one_init=True, real=False):
    r"""Generates a banded unitary matrix.
    Args:
        N (int): number of modes
        w (int): bandwidth
        top_one_init (bool): if True places a `1\times1` interferometer in the top-left-most block of the first matrix in the product
        real (bool): return a random real orthogonal matrix
    Returns:
        array: random :math:`N\times N` unitary with the specified block structure
    """
    if N < w + 1:
        raise ValueError("The bandwidth can be at most one minus the size of the matrix.")
    if N == w + 1:
        return random_interferometer(N, real=real)
    U = sp.linalg.block_diag(*[random_interferometer(1, real=real) for _ in range(N)])
    for _ in range(w):
        U = U @ random_block_interferometer(N, top_one=top_one_init, real=real)
        top_one_init = not top_one_init
    return U