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.