def find_scale_index(level, x='0', y='1'):
    """
    Gray code order is used here to generate indices for the scales of each decomposition level in a
    wavelet packet transform. Either 'a' and 'd' (approximation and detail) or '0' and '1' are used.
    E.g., for a level=1 decomposition, the possible combination will be 'ad' and 'da'. For level=2,
    it will be 'aa', 'ad', 'da', 'dd'. The indices are also ordered by their corresponding frequency
    bandpasses.

    Parameters
    ----------
    level : int
        Dcomposition level for WPT.
    x : '0', str, optional
        First index in the gray code. 'a' is also generally used.
    y : '1', str, optional
        Second index in the gray code. 'd' is also generally used.

    Returns
    -------
    graycode_order : list
        The list of indices ordered by their corresponding frequencies.
    """

    graycode_order = [x, y]
    for _ in range(level - 1):
        graycode_order = [x + path for path in graycode_order] + \
                         [y + path for path in graycode_order[::-1]]
    return graycode_order