How to use the kymatio.scattering1d.backend.subsample_fourier function in kymatio

S = {}

    # pad to a dyadic size and make it complex
    U0 = pad(x, pad_left=pad_left, pad_right=pad_right, to_complex=True)
    # compute the Fourier transform
    U0_hat = fft1d_c2c(U0)
    if vectorize:
        # initialize the cursor
        cc = [0] + list(size_scattering[:-1])  # current coordinate
        cc[1] = cc[0] + cc[1]
        if max_order == 2:
            cc[2] = cc[1] + cc[2]
    # Get S0
    k0 = max(J - oversampling, 0)
    if average:
        S0_J_hat = subsample_fourier(U0_hat * phi[0], 2**k0)
        S0_J = unpad(real(ifft1d_c2c(S0_J_hat)),
                     ind_start[k0], ind_end[k0])
    else:
        S0_J = x
    if vectorize:
        S[:, cc[0], :] = S0_J.squeeze(dim=1)
        cc[0] += 1
    else:
        S[()] = S0_J
    # First order:
    for n1 in range(len(psi1)):
        # Convolution + downsampling
        j1 = psi1[n1]['j']
        k1 = max(j1 - oversampling, 0)
        assert psi1[n1]['xi'] < 0.5 / (2**k1)
        U1_hat = subsample_fourier(U0_hat * psi1[n1][0], 2**k1)

S[()] = S0_J
    # First order:
    for n1 in range(len(psi1)):
        # Convolution + downsampling
        j1 = psi1[n1]['j']
        k1 = max(j1 - oversampling, 0)
        assert psi1[n1]['xi'] < 0.5 / (2**k1)
        U1_hat = subsample_fourier(U0_hat * psi1[n1][0], 2**k1)
        # Take the modulus
        U1 = modulus_complex(ifft1d_c2c(U1_hat))
        if average or max_order > 1:
            U1_hat = fft1d_c2c(U1)
        if average:
            # Convolve with phi_J
            k1_J = max(J - k1 - oversampling, 0)
            S1_J_hat = subsample_fourier(U1_hat * phi[k1], 2**k1_J)
            S1_J = unpad(real(ifft1d_c2c(S1_J_hat)),
                         ind_start[k1_J + k1], ind_end[k1_J + k1])
        else:
            # just take the real value and unpad
            S1_J = unpad(real(U1), ind_start[k1], ind_end[k1])
        if vectorize:
            S[:, cc[1], :] = S1_J.squeeze(dim=1)
            cc[1] += 1
        else:
            S[(n1,)] = S1_J
        if max_order == 2:
            # 2nd order
            for n2 in range(len(psi2)):
                j2 = psi2[n2]['j']
                if j2 > j1:
                    assert psi2[n2]['xi'] < psi1[n1]['xi']

# just take the real value and unpad
            S1_J = unpad(real(U1), ind_start[k1], ind_end[k1])
        if vectorize:
            S[:, cc[1], :] = S1_J.squeeze(dim=1)
            cc[1] += 1
        else:
            S[(n1,)] = S1_J
        if max_order == 2:
            # 2nd order
            for n2 in range(len(psi2)):
                j2 = psi2[n2]['j']
                if j2 > j1:
                    assert psi2[n2]['xi'] < psi1[n1]['xi']
                    # convolution + downsampling
                    k2 = max(j2 - k1 - oversampling, 0)
                    U2_hat = subsample_fourier(U1_hat * psi2[n2][k1],
                                               2**k2)
                    # take the modulus and go back in Fourier
                    U2 = modulus_complex(ifft1d_c2c(U2_hat))
                    if average:
                        U2_hat = fft1d_c2c(U2)
                        # Convolve with phi_J
                        k2_J = max(J - k2 - k1 - oversampling, 0)
                        S2_J_hat = subsample_fourier(U2_hat * phi[k1 + k2],
                                                     2**k2_J)
                        S2_J = unpad(real(ifft1d_c2c(S2_J_hat)),
                                     ind_start[k1 + k2 + k2_J],
                                     ind_end[k1 + k2 + k2_J])
                    else:
                        # just take the real value and unpad
                        S2_J = unpad(
                            real(U2), ind_start[k1 + k2], ind_end[k1 + k2])

# 2nd order
            for n2 in range(len(psi2)):
                j2 = psi2[n2]['j']
                if j2 > j1:
                    assert psi2[n2]['xi'] < psi1[n1]['xi']
                    # convolution + downsampling
                    k2 = max(j2 - k1 - oversampling, 0)
                    U2_hat = subsample_fourier(U1_hat * psi2[n2][k1],
                                               2**k2)
                    # take the modulus and go back in Fourier
                    U2 = modulus_complex(ifft1d_c2c(U2_hat))
                    if average:
                        U2_hat = fft1d_c2c(U2)
                        # Convolve with phi_J
                        k2_J = max(J - k2 - k1 - oversampling, 0)
                        S2_J_hat = subsample_fourier(U2_hat * phi[k1 + k2],
                                                     2**k2_J)
                        S2_J = unpad(real(ifft1d_c2c(S2_J_hat)),
                                     ind_start[k1 + k2 + k2_J],
                                     ind_end[k1 + k2 + k2_J])
                    else:
                        # just take the real value and unpad
                        S2_J = unpad(
                            real(U2), ind_start[k1 + k2], ind_end[k1 + k2])
                    if vectorize:
                        S[:, cc[2], :] = S2_J.squeeze(dim=1)
                        cc[2] += 1
                    else:
                        S[n1, n2] = S2_J

    return S

S0_J = unpad(real(ifft1d_c2c(S0_J_hat)),
                     ind_start[k0], ind_end[k0])
    else:
        S0_J = x
    if vectorize:
        S[:, cc[0], :] = S0_J.squeeze(dim=1)
        cc[0] += 1
    else:
        S[()] = S0_J
    # First order:
    for n1 in range(len(psi1)):
        # Convolution + downsampling
        j1 = psi1[n1]['j']
        k1 = max(j1 - oversampling, 0)
        assert psi1[n1]['xi'] < 0.5 / (2**k1)
        U1_hat = subsample_fourier(U0_hat * psi1[n1][0], 2**k1)
        # Take the modulus
        U1 = modulus_complex(ifft1d_c2c(U1_hat))
        if average or max_order > 1:
            U1_hat = fft1d_c2c(U1)
        if average:
            # Convolve with phi_J
            k1_J = max(J - k1 - oversampling, 0)
            S1_J_hat = subsample_fourier(U1_hat * phi[k1], 2**k1_J)
            S1_J = unpad(real(ifft1d_c2c(S1_J_hat)),
                         ind_start[k1_J + k1], ind_end[k1_J + k1])
        else:
            # just take the real value and unpad
            S1_J = unpad(real(U1), ind_start[k1], ind_end[k1])
        if vectorize:
            S[:, cc[1], :] = S1_J.squeeze(dim=1)
            cc[1] += 1