function Psi = get_Psi(N, M, epi)
    Psi = zeros(N,M*N);
    for n = 0 : N-1
        Cn = floor(rand()*M);
        for q = 0 : N-1
            for p = 0:M-1
                Psi(n+1,q*M+p+1) = exp(1i*2*pi*p/M*Cn+1i*2*pi*q/N*n*(1+Cn*epi)) / sqrt(N);
            end
        end
    end
end