- How can we measure meson decay constants?
I am not an experimental physicst, but I think that the best way to obtain the decay constant is to study processes like $\pi^+\to \mu^+ \nu$ and extract them from the branching ratio:
$$\rm{Br}(\pi^+\to \mu^+\nu)=\dfrac{G_F^2 m_{\pi^+} m_\mu^2}{8 \pi}\left(1-\dfrac{m_\mu^2}{m_{\pi^+}^2} \right)^2 f_{\pi^+}^2 |V_{ud}|^2 \tau_{\pi^+} ,$$
which is measured nowadays with great precision.
@dmckee answer's suggests that we can also extract the decay constant from the pion form factor, but this method seems less precise, because it is more difficult to measure form factors than decay constants (but maybe I'm wrong...). If you take a look at PDG, you'll see that the process $\pi^+\to \mu^+\nu$ is measured with an incredible precision.
One last comment about decay constants: actually, these quantities can be computed for pions using Lattice QCD methods and the theoretical error bars are comparable to the experimental ones! You can even find very precise computations for more exotic mesons, like $D$, $B$ and $B_s$.
- For your theoretical question:
It depends on the process you are considering! For example, if you have $\pi^+\to \mu^+\nu$, then you must take a second order term. In this term, you need a current $J^\mu_{q}$ related to the annihilation $u \bar{d}\to W^+$ and a leptonic current $J^\mu_\ell$ related to the creation $W^+\to\mu^+\nu$. Then, the time ordered product will only apply to the $W^+W^-$ term and it will give you simply the $W^+$ boson propagator.
From my experience, I would suggest you to integrate-out the vector bosons, because the corrections to the fermi theory are negligible. In this case, you can write an effective Hamiltonian:
$$\mathcal{H}_{\text{eff}}=-\sqrt{2} G_F V_{ub} [\overline{u}\gamma_\mu (1-\gamma_5)d][\bar{\mu}_L \gamma^\mu {\nu_\mu}_L] +\text{h.c.}, $$
and it is much simpler to read the amplitude and to relate it with the decay constant, because the hadronic part factorizes:
$$\mathcal{A}=-i\langle \mu^+,\nu | {H}_{\text{eff}} |\pi^+\rangle =i\sqrt{2} G_F V_{ub} \langle 0 |\overline{u}\gamma_\mu \gamma^5 d|\pi^+\rangle\cdot \bar{u}(p_\nu)[\gamma^\mu(1-\gamma_5)/2 ]v(p_\mu),$$
where $$\langle 0 |\overline{u}\gamma_\mu \gamma^5 d|\pi^+\rangle=-i p_\mu f_{\pi^+}.$$