First forget QFT for a while and think about Classical Field Theory. Consider the Klein-Gordon field more precisely. Its Lagrangian is
$$\mathcal{L}(\phi,\partial_\mu\phi)=\dfrac{1}{2}\partial^\mu\phi \partial_\mu \phi-\dfrac{1}{2}m^2\phi^2$$
In this Lagrangian the variable is $\phi$. Now since $\phi$ is a function defined on spacetime, $\phi$ depends on time in particular and you can compute in a reference frame $\dot{\phi} = \partial_0 \phi$.
One then defines the conjugate momentum, as in Classical Mechanics, to be
$$\pi = \dfrac{\partial \mathcal{L}}{\partial(\partial_0\phi)}.$$
For this field, what do we get? If you work this out you'll find $\pi = \dot{\phi}$.
This is all classical. Then with this in hands you can quantize.
After all, quantizing the field means that you want to turn $\phi,\pi$ into operators obeying
$$[\phi(x),\phi(y)]=[\pi(x),\pi(y)]=0$$
$$[\phi(x),\pi(y)]=i\delta(x-y).$$
Thus you already need $\pi$ to talk about quantization. As in Quantum Mechanics, you need both position and momentum to impose the canonical commutation relations.
By the way, there's a small detail. The commutation relations are taken at equal times. In that case they hold among the Schrodinger picture operators $\phi(\mathbf{x}),\pi(\mathbf{y})$, since they are defined at the same initial time.
So if you want to compute $\pi$ from $\phi$, you can do it classicaly and then impose the canonical commutation relations, or you can do it in the Heisenberg picture and you'll get the same results.
Edit: the mode decomposition can be achieved in Classical Field Theory, the only thing is that the coefficients will be numbers. The equation of motion is
$$(\Box + m^2)\phi=0$$
Take the Fourier transform in the spatial variable so that denoting the Fourier transform by $\hat{\phi}$ you have
$$\partial^2_t \hat{\phi} +(|\mathbf{p}|^2+m^2)\hat{\phi}=0$$
define $\omega_{p}^2=|\mathbf{p}|^2+m^2$ and $p = (\omega_p, \mathbf{p})$. The equation is parametrized by $\mathbf{p}$ and can easily be solved to give
$$\hat{\phi}(\mathbf{p},t)=a_p e^{-i\omega_p t}+b_p e^{i\omega_p t}$$
now apply the reality condition of the Fourier transform
$$\hat{\phi}(-\mathbf{p},t)=\hat{\phi}(\mathbf{p},t)^\ast.$$
You arrive at the condition
$$a_{-\mathbf{p}}e^{-i\omega_p t}+b_{-\mathbf{p}}e^{i\omega_p t}=a_{\mathbf{p}}^\ast e^{i\omega_p t}+b_{\mathbf{p}}^\ast e^{-i\omega_p t}$$
the linear independence of the exponentials then gives $a_{-\mathbf{p}}=b_{\mathbf{p}}^\ast$ and $b_{-\mathbf{p}}=a_{\mathbf{p}}^\ast$. Now you have
$$\hat{\phi}(\mathbf{p},t)=a_pe^{-i\omega_p t}+a_{-p}^\ast e^{i\omega_p t}$$
now apply the Fourier inversion to get
$$\phi(x)=\int \dfrac{d^3 p}{(2\pi)^3}(a_p e^{-i\omega_p t}+a_{-p}^\ast e^{i\omega_p t}) e^{i \mathbf{p}\cdot \mathbf{x}}$$
if you change variables on the second term making $p \to -p$ since $\omega_{-p}=\omega_p$ you get the formula
$$\phi(x)=\int \dfrac{d^3 p}{(2\pi)^3}(a_p e^{-ipx}+a_p^\ast e^{ipx}).$$
The $\sqrt{2\omega_p}$ is then included for convenience to get a Lorentz invariant result (it amounts for a redefinition of $a_p$). The final answer is
$$\phi(x)=\int \dfrac{d^3 p}{(2\pi)^3}\dfrac{1}{\sqrt{2\omega_p}}(a_p e^{-ipx}+a_p^\ast e^{ipx}).$$
Please, understand that this doesn't derive the Fock space representation. This is just a classical calculation that in turns motivates the mode decomposition in terms of the Fock space ladder operators.
By the way, there's a cleaner and more elegant approach with the spacetime Fourier transform that can be found in the question A question on using Fourier decomposition to solve the Klein Gordon equation.
