Notation: $x=(t,\boldsymbol x)$; $k=(k_0,\boldsymbol k)$; $kx=k_0t-\boldsymbol k\cdot\boldsymbol x$; $\mathrm dx=\mathrm dt\;\mathrm d^3\boldsymbol x$; etc.
You can in principle perform the Fourier decomposition on both space and time variables, but to do so you'll need several properties of the Dirac's delta funciton:
The first one is: let $\xi\in\mathbb R$; then
$$
\delta(f(\xi))=\sum_{f(\xi_i)=0} \frac{\delta(\xi-\xi_i)}{|f'(\xi_i)|} \tag{1}
$$
where the sum is over every $\xi_i$ such that $f(\xi_i)=0$, ie, over the roots of $f(\xi)$.
The second one is that, given $g(\xi)$ a known function, the distributional solution of $g(\xi)f(\xi)=0$ is $f(\xi)=h(\xi)\delta(g(\xi))$ for an arbitrary function $h(\xi)$. If you believe these, then the Fourier decomposition is as follows:
Let $\phi(x)$ be the solution of
$$
(\partial^2+m^2)\phi(x)=0
$$
Take the Fourier transform of the equation to find
$$
(k^2-m^2)\phi(k)=0 \tag{2}
$$
where
$$
\phi(k)=\int \mathrm dx\ \mathrm e^{ikx} \phi(x)
$$
As $\phi(x)$ is a distribution, the solution of $(2)$ is $\phi(k)=h(k)\delta(k^2-m^2)$ for an arbitrary function $h(k)$. Inverting the Fourier Transform, we find
$$
\phi(x)=\int\mathrm dk\ \mathrm e^{-ikx}h(k)\delta(k^2-m^2)
$$
Next, use $(1)$ to expand the delta over the roots of $k^2-m^2$. These roots are easily found to be $k_0=\pm \omega(\boldsymbol k)$, where $\omega(\boldsymbol k)=+(\boldsymbol k^2+m^2)^{1/2}$. Therefore, it is immediate to get
$$
\phi(x)=\int\mathrm dk\ \mathrm e^{-ikx}h(k)\frac{1}{2\omega}\left[\delta(k_0-\omega)+\delta(k_0+\omega)\right]
$$
and, after integrating over $\mathrm dk_0$ using the deltas, we find
$$
\phi(x)=\int\frac{\mathrm d \boldsymbol k}{2\omega}\ \left[\mathrm e^{-i\omega t} \mathrm e^{i\boldsymbol k\cdot\boldsymbol x}h(\omega,\boldsymbol k)+\mathrm e^{+i\omega t} \mathrm e^{i\boldsymbol k\cdot\boldsymbol x}h(-\omega,\boldsymbol k)\right]
$$
Finally, make the change of variable $\boldsymbol k\to-\boldsymbol k$ in the second term, which yeilds the usual expansion
$$
\phi(x)=\int\frac{\mathrm d \boldsymbol k}{2\omega}\ \left[\mathrm e^{-ikx}a(\boldsymbol a)+\mathrm e^{+ikx}b^\dagger(\boldsymbol k)\right]
$$
where I have defined $a(\boldsymbol k)=h(\omega,\boldsymbol k)$ and $b^\dagger(\boldsymbol k)=h(-\omega,-\boldsymbol k)$.
As you can see, the solution is the same as yours (modulo some irrelevant prefactor that can be reabsorbed into the definition of $h(k)$), though the algebraic procedure to find it is a bit harder.
Best Answer
The important insight is that it's actually the whole combination $$ \frac{d^3 k}{2(2\pi)^3 E_\mathbf k}, \qquad E_\mathbf{k} = \sqrt{\mathbf k^2 + m^2} $$ that forms a Lorentz-invariant measure. To see this, note that if we define $k= (k^0, \mathbf k)$ and use the identity $$ \delta(f(x)) = \sum_{\{x_i:f(x_i) = 0\}} \frac{\delta(x-x_i)}{|f'(x_i)|} $$ then we get $$ \delta(k^2 - m^2)=\frac{\delta(k^0 - \sqrt{\mathbf k^2+m^2})}{2\sqrt{\mathbf k^2+m^2}} + \frac{\delta(k^0 + \sqrt{\mathbf k^2+m^2})}{2\sqrt{\mathbf k^2+m^2}} $$ so the original measure can be rewritten as $$ \frac{d^3 k}{2(2\pi)^3 E_\mathbf k}=\frac{d^3k\,d k^0}{2(2\pi)^3 k^0}\delta(k^0 - \sqrt{\mathbf k^2+m^2}) = \frac{d^4k}{(2\pi)^3}\delta(k^2-m^2)\theta(k^0) $$ which is manifestly Lorentz invariant for proper, orthochronous Lorentz transformations. The rest of your manipulations go through unscathed, and you get the result you want!
Hope that helps!
Cheers!