Here's a simple algorithm that was proposed in Dixon's 1970 paper: http://www.ams.org/journals/mcom/1970-24-111/S0025-5718-1970-0280611-6/S0025-5718-1970-0280611-6.pdf
It's split into 2 parts:
1. Construct a non-scalar commuting matrix $H$ (if possible)
Dixon presents two variants for part 1. I'll just mention the simpler one here. Let $\rho:G \to \text{GL}(n,\mathbb{C})$ be a (unitary) representation. For $r,s = 1,2,\dots,n$, define
$$
H_{rs} = \begin{cases}
E_{rr} &\text{if } r = s \\
E_{rs} + E_{sr} &\text{if } r > s \\
i(E_{rs} - E_{sr}) &\text{if } r < s,
\end{cases}
$$
where $E_{rs}$ is the $n \times n$ matrix with 1 in the $(r,s)^{th}$ entry and 0 everywhere else. Then $\{H_{rs}\}_{r,s=1}^n$ forms a Hermitian basis for the $n \times n$ matrices.
Now for each $r,s$, compute the sum
$$
H = \frac{1}{|G|} \sum_{g \in G} \,\, \rho(g)^* \, H_{rs} \, \rho(g)
$$
Observe that $H$ commutes with all $\rho(g)$.
If $\rho$ is irreducible, then $H$ is a scalar matrix for all $r,s$. Otherwise, there will be some $r,s$ such that $H$ is non-scalar (because $\{H_{rs}\}_{r,s=1}^n$ forms a basis). In this case, proceed to part 2:
2. Use the eigenspaces of $H$ to decompose $\rho$
Let $H = UJU^*$ be the Jordan decomposition of $H$, where $U$ is unitary. Then
$$
U^* \rho(g) U
$$
will have the same block-diagonal form for all $g \in G$, yielding a decomposition of $\rho$.
Of course this assumes that the $n$ and $|G|$ are small enough that one can easily run through all the $H_{rs}$ and compute the Jordan decomposition. There are more sophisticated and efficient methods, but I like the simplicity of this one.
I've written an implementation of this algorithm in Sage.
