This is a modification of the Gram-Schmidt process. As usual, the projection of $x$ onto $a$ is defined as $\frac{\langle x,a\rangle}{\langle a,a\rangle}a$, and the projection of $x$ away from $a$ is $x-\frac{\langle x,a\rangle}{\langle a,a\rangle}a$.
Take any basis vector $a\neq0$. If $\langle a,a\rangle\neq0$, then project the other $n-1$ basis vectors away from $a$, and continue the process in that subspace. If $\langle a,a\rangle=0$ and there's another basis vector $b$ not orthogonal to $a$, then project the other $n-2$ basis vectors away from the span of $a$ and $b$ (which is hyperbolic), and continue the process in that subspace. If $\langle a,a\rangle=0$ and all other basis vectors are orthogonal to $a$, then nothing needs to be done in this step; continue the process in the span of the other basis vectors.
(And any hyperbolic plane produced in the process can be given an orthonormal basis. Given $\langle a,a\rangle=0\neq\langle b,a\rangle$, define $b'=\frac{b}{\langle b,a\rangle}-\frac{\langle b,b\rangle a}{2\langle b,a\rangle^2}$, so that $\langle b',b'\rangle=0$ and $\langle b',a\rangle=1$. Then take $c=a+\tfrac12b'$ and $d=a-\tfrac12b'$ as the new basis, with $\langle c,c\rangle=1=-\langle d,d\rangle$ and $\langle c,d\rangle=0$.)
Note that this works as well with degenerate forms. It always produces an orthogonal basis.