My thoughts / background information:
It is easy to find an inclusion homomorphism for complex matrices into real matrices: considering the one dimensional case note that multiplying a complex number $z$ by a complex number of unit length is the same as rotation the complex number $z$ by a certain angle. Writing down the rotation matrix in $\mathbb R^2$ makes it obvious that mapping $a + ib$ to $\displaystyle \left(\begin{array}{cc} a & b \\ -b & a \end{array}\right)$ is the desired inclusion homomorphism in one dimension and generalising it to $n$ dimensions is not difficult.
My question:
What is not clear to me is generalising this to quaternionic matrices: given a quaternion $p = a + jb$ (caution! $a$ and $b$ are now in $\mathbb C$) and a unit quaternion $q = \cos {\theta \over 2} + (u_x i + u_y j + u_z k)\sin {\theta \over 2}$ I am aware that $qpq^{-1}$ is a rotated quaternion by angle $\theta$.
But how, from that, can I derive that the homomorphism in $1$ dimension should be
$$ a + bj \mapsto \left(\begin{array}{cc} a & b \\ -\overline{b} & \overline{a} \end{array}\right)$$
How to derive this map?
I think you're mislead by focussing on the "rotation" aspect, instead of just the operation of "multiplication by a quaternion (or complex number)".
For a fixed quaternion $a+bj$, the map $$ (x+yj) \mapsto (x+yj)(a+bj) $$ is a $\mathbf{C}$-linear map from $\mathbf{H}$ to $\mathbf{H}$. Here we are viewing $\mathbf{H}$ as a two-dimensional complex vector space where scalar multiplication is from the left; this is an important detail since if $\lambda \in \mathbf{C}$ and $q \in \mathbf{H}$, then $\lambda q \neq q \lambda$ in general.
Multiplying out the right-hand side and simplifying (using that $jz=\bar{z}j$ for complex numbers $z$), we get $(xa-y\bar{b})+(xb+y\bar{a})j$, which means that the matrix representation of this linear map is $$ (x,y) \mapsto (x,y) \begin{pmatrix} a & b \\ -\bar{b} & \bar{a} \end{pmatrix} . $$ (This looks a little unusual, since normally we let matrices act on column vectors rather than row vectors, but it has to do with the choice of putting scalars on the left.)
Now, composing two such maps shows that $$ (x+yj) \mapsto (x+yj)(a+bj)(c+dj) $$ corresponds to $$ (x,y) \mapsto (x,y) \begin{pmatrix} a & b \\ -\bar{b} & \bar{a} \end{pmatrix} \begin{pmatrix} c & d \\ -\bar{d} & \bar{c} \end{pmatrix} , $$ so the product of the matrices indeed corresponds to the product of the corresponding quaternions.