(Edit: the original claim was much more ambitious).
It is possible to transform the representation into a representation which is slightly less general, without extending the field.
Let $\mu:G\to GL_d(\mathbb{R})$ s.t. $\rho(G)=M\mu(G)M^{-1}$ for $G\in GL_d(\mathbb{C})$. W.l.o.g. $\mu$ is orthogonal, i.e. $hh^\top=I$ for any $h\in\mu(G)$. Thus $$|G|I=\sum_{h\in\mu(G)} hh^\top=\sum_{h\in\mu(G)} hh^*= \sum_{g\in\rho(G)} M^{-1}gM (M^{-1}gM)^*=\sum_{g\in\rho(G)} M^{-1}gM M^*g^*(M^*)^{-1}.$$ Multiplying on the right by $M^*$ and on the left by $M$ one obtains $$|G|MM^*=\sum_{g\in\rho(G)} gM M^*g^*.$$
Up to scaling by a positive real, there is unique, by irredicibility of $\rho$, Hermitean p.s.d. matrix $H$ providing $G$-invariant sesquilinear scalar product, i.e. $gHg^*=H$ for any $g\in\rho(G)$. Therefore $H=MM^*$. Note that $H$ only records the info about $M$ up to a unitary factor. Namely, $M=PU$, with $U$ unitary and $P$ Hermitean positive definite ($PU$ is called left polar decomposition of $M$). Thus $H=MM^*=PUU^*P^*=PP^*=P^2$.
Note that $H$ affords a (unique, in fact), $LDL^*$ decomposition, with $D=diag(D_1,\dots,D_d)$, and $D_k>0$ for all $1\leq k\leq d$.
It is well-known that $L$ and $D$ can be computed by solving a system of linear equations, without extending the field (e.g. as on p.43 of our preprint (=MSc thesis of my student K.Hymabaccus)).
Set $M=LD^{1/2}U$. Then $$\mu(G)=M^{-1}\rho(G)M=U^*D^{-1/2}L^{-1}\rho(G)LD^{1/2}U.$$ Thus we can perform the conjugation by $L$, and so reduce to the case of $M=D^{1/2}U$.