Quick and practical, but doesn't properly take into account repeated median values, e.g. [1, 2, 2, 2, 2] would have the first 3 quartiles equal to 2, but this code would give [1, 2, NaN].

Nice code, thanks! Just for fun, slight to medium speed improvement (sorry my k = your N):

%instead of:
%[Y{k:-1:1}] = ndgrid(X) ;
%Y = reshape(cat(k+1,Y{:}),[],k) ;
%fill Y directly:
nX = numel(X);
Y(nX^k, k) = X(1); % Initialize to give shape and data type.
s(1:k-1) = {ones(1, nX)};
for i = 1:k
Y( ((i-1)*nX^k + 1) : (i*nX^k) ) = X(s{1:i-1}, :, s{i:k-1});
X = permute(X, [1:i-1, i+1, i]);
end
% This has inverse order for the combinations as you had, easy to change.

Theoretical note, in case someone (as I did) wants to get a limited number of random combinations, then this method fails for N > 51 since then eps(2^52) = 1, i.e. can't represent the integers accurately any more with double type. Of course, here we run out of memory with such a large N anyway.

Also, with bitget, I was using a loop over bits instead of a loop over combinations. I haven't checked, but I assumed it would be faster. I'm also going to check with d2b instead of bitget. (found on FEX, and modified to work on arrays).

Actually, the bug was for any power of 2 (1, 2, 4, ...). Ali's modification works, but a more elegant fix is:
c = floor(log2(x)) + 1;
No need for any other modification from original code (no special case for 1, no need to remove preceding zeros).