function Psi = get_Psi_KR_product(N, M, epi)
    d_n = floor(rand(N, 1)*M) / M;
    zeta_n = 1 + d_n .* epi; % epi = B / f_c

    R = zeros(N, M);
    D = zeros(N, N);
    for n = 0: N-1
        for m = 0: M-1
            R(n+1, m+1) = exp(-1j * 2 * pi * m * d_n(n+1));
        end
    end

    for n = 0: N-1
        for l = 0: N-1
            D(n+1, l+1) = exp(-1j * 2 * pi * l * n * zeta_n(n+1) / N);
        end
    end

    Psi = KR_product(R', D')';
end

function kr = KR_product(F, G)
    nR_F = size(F, 1);
    nR_G = size(G, 1);
    mul = ones(nR_G, 1);
    FF = kron(F, mul);
    GG = repmat(G, nR_F, 1);
    kr = FF .* GG;
end