In MatLab I generate Haar distributed matrices like this:
m = randn(m,m)% all elements are normally distributed
u= qr(m) % make qr decomposition and what you get is Haar measure on "u"
So the mathematical statement is that if m is normally distributed, then u is Haar.
The reason is quite trivial — normal distribution is preserved by unitary transformations.
However writing this I begin to doubt myself about tiny details - depending
how they implement qr algorithm the matrix u is not unique, it can be multiplied by
$\operatorname{diag}(\pm1)$. Neverthelss most probably everything should be correct.
