function x=randnND(mu,Sigma,N)
%%  x=randN(mu,Sigma,N)
%% mu, sigma: parameters of the Gaussian
%% N: how many samples should be generated
%% 

   mu = mu(:);
   l = length(mu);
   [V,D] = eig(Sigma); %% Calculate the transformation matrix
   ss = sqrt(diag(D));
   T =  V*diag(ss);
   x = T*randn(l,N) + mu*ones(1,N);