function [C_n, Psi] = get_Psi(N, M, epi)
    if nargin == 2
        epi = 0;
    end
    k = 1 / sqrt(M * N);
    C_n = zeros(M, 1);

    % Determined code 1
    % C_n = repelem(0:M-1, N/M);

    % Determined code 2
    % C_n = repmat(0:M-1, 1, N/M);

    % Determined code 3
    C_n = repmat(0:M-1, 1, N/M);
    C_n = C_n(randperm(length(C_n)));

    Psi = zeros(N,M*N);
    for n = 0 : N-1
        col = zeros(1, M*N);
        % C_n(n+1) = floor(rand()*M);
        for q = 0 : N-1
            for p = 0:M-1
                f1 = p/M*C_n(n+1);
                f2 = q/N*n*(1+C_n(n+1)*epi);
                col(q*M+p+1) = exp(1j*2*pi*(f1 + f2));
            end
        end
        Psi(n+1, :) = col * k;
    end
end