I did not thought much about that issue, so my answer may not be the best one, but it has the advantage of being quite simple to understand, it is exact even if you restrict yourself to a discrete set of universal gates, and can be generalized to more complicated relations. On the other side, this is not a unitary process because you need to repeat it a polynomial number of times until you get your expected states (see maybe Gianni version for a ).

Let $k = \lceil\log N\rceil$ be the smallest number of bits used to represent $N$. Then, you create the state
$$\sum_{x =0}^{2^k-1} |x\rangle$$
by applying $k$ Hadamard gates on $k$ states initialized to $|0\rangle$. Then, if we define $f(x) = 1$ iff $x < N$ (here you can put any relation you like depending on the states you are interested in, you just want to have a constant fraction of $x$ such that $f(x)=1$ to ensure the algorithm is efficient, which is the case here because at least half of the $x$ have this property by definition of $k$). Then, you add one auxiliary qubit, and you apply the unitary associated to $f$ (which is $U_f |x\rangle|b\rangle =|x\rangle|b \oplus f(x)\rangle$) in order to obtain the state:
$$\sum_{x=0}^{2^k-1}|x\rangle |f(x)\rangle$$
Then, you measure this second register. If you obtain a $1$, the state you obtain is:
$$\sum_{x \in \{x \mid f(x)=1\}} |x\rangle = \sum_{x = 0}^{N-1} |x\rangle$$
(so you have the expected state) Otherwise, you start again from scratch (the probability to obtain the expected state depends on the fraction of $x$ such that $f(x) = 1$. In our case, this is lower bounded by $1/2$ by definition of $k$ so you need to repeat it only a few times before you succeed)