Definitions
We will perform this in general dimension $d$, setting d=4 only when required; this makes some steps clearer. Let us begin with a definition of what it we want to calculate; note that we use the mostly negative metric signature ($+-\ldots-$ here).
\begin{align}
\Delta_F(x) &= \int \frac{d^d k}{(2\pi)^d} \frac{1}{k^2-m^2+i\epsilon}e^{-i k x}
\end{align}
Our goal is to make the $k^0$ integral equal to a contour integral, which we evaluate via the residue theorem, rather than manually evaluating it. We then do the easy (d-2) spatial integrals, leaving the final one in a form that we recognise as the integral representation of a special function, which for $d=4$ is $K_1(z)$.
We will also assume that $x^2 \ge 0$; the $x^2 <0$ case is similar (see Huang's QFT, section 2.9).
Contour integration
We will close the real $k^0$ integration from $(-\infty, \infty)$ via an infinitely large semicircle. Ideally we want the contribution of these semicircles to be zero, such that
$$\int_{-\infty}^\infty dk^0 \, f(k) = \int_{-\infty}^\infty dk^0 \, f(k) + \int_{\text{semicircle}} dk^0 \, f(k) = \int_{\gamma} dk^0 \, f(k) = \text{ residues at poles}.$$
In our case, this will only happen if we have $\mathrm{Re}[-ik_0 x_0] < 0$ for the $k_0$ points on this semicircle, so that the exponent does not blow up; this means that we demand that $\mathrm{Im}[k_0]x_0 <0$. This means that:
- If $x^0$ > 0, we need to close where $\mathrm{Im}[k_0] <0$, i.e. the lower half-plane.
- If $x^0$ < 0, we need to close where $\mathrm{Im}[k_0] > 0$, i.e. the upper half-plane.
In both cases, the semicircular arc contributes zero to the integral, by construction.
We need to choose different contours depending on the sign of $x^0$; we will implement this by the step function $\theta$. The two contours completed by semicircles above and below the line are $\gamma^+$ and $\gamma^-$, and they enclose the poles at $k=-\omega$ and $\omega$ respectively, where $\omega = \sqrt{|\vec{k}|^2+m^2-i\epsilon}$.
\begin{align}
\Delta_F (x)&= \theta(x^0) \int_{\mathbb{R}^{d-1}}\frac{d^{d-1} k}{(2\pi)^{d}} \int_{\gamma^-} d k^0 \frac{e^{-ikx}}{k^2-m^2+i\epsilon}+\theta(-x^0) \int_{\mathbb{R}^{d-1}}\frac{d^{d-1} k}{(2\pi)^{d}} \int_{\gamma^+} d k^0\frac{e^{-ikx}}{k^2-m^2+i\epsilon}\\
&= \Delta^{+}_F(x) + \Delta^{-}_F(x)\\
\Rightarrow \Delta^{\pm}_F (x)&=\theta(\pm x^0) \int_{\mathbb{R}^{d-1}}\frac{d^{d-1} k}{(2\pi)^d} \int_{\gamma^\mp } d k^0 \frac{1}{2\omega}\left[\frac{1}{k^0-\omega}-\frac{1}{k^0+\omega} \right]e^{-ikx}
\end{align}
This is a trivial exercise in contour integration for $k^0$.
\begin{align}
\Delta^{\pm}_F (x)&=-2\pi i \theta(\pm x^0) \int_{\mathbb{R}^{d-1}} \frac{d^{d-1}k}{(2\pi)^{d}}\frac{1}{2\omega}e^{\mp i \omega x^0}e^{i\vec{k} \cdot \vec{x}}\\
&=-i \theta(\pm x^0) \int_{\mathbb{R}^{d-1}} \frac{d^{d-1} k}{(2\pi)^{d-1}}\frac{1}{2\omega}e^{- i \omega |x^0|}e^{i\vec{k} \cdot \vec{x}}\\
\end{align}
Adding the two back together, since the expressions are identical the $\theta$s disappear:
\begin{align}
\Delta_F (x)&=-\frac{i}{2}\int_{\mathbb{R}^{d-1}} \frac{d^{d-1} k}{(2\pi)^{d-1}}\frac{1}{\omega}e^{- i \omega |x^0|}e^{i\vec{k} \cdot \vec{x}}
\end{align}
Rewriting the spatial integrations
Now let's cheat, and use the fact that the answer must be Lorentz symmetric, and therefore only a function of $x^2$. Assuming that $x^2 >0$, we can therefore boost to a frame where $\vec{x}=0$, evaluate the integral, and boost back later. We will use spherical coordinates in $(d-1)=3$ dimensions, with $l=|\vec{k}|$, where we note that
$$
\int d \Omega_{d-2} = \Omega_{d-2} = \frac{2 \pi^{\frac{d-1}{2}}}{\Gamma(\frac{d-1}{2})}\, \text{= surface area of a $(d-2)$-sphere}
$$
This does not actually converge, so let's shift $|x^0| \to |x^0| -i\delta$, $\delta \rightarrow 0$; this to ensure convergence of the integral, by making the integrand vanish for large $l$. We are left with:
\begin{align}
\Delta_F(x)&=-\frac{i}{2} \frac{1}{(2\pi)^{d-1}}\left(\int d \Omega_{d-2}\right) \lim_{\delta \rightarrow 0} \int_{0}^\infty dl \frac{l^{d-2}}{\sqrt{l^2+m^2}}e^{-i \sqrt{l^2+m^2}(|x^0| - i \delta)}\\
&=-\frac{1}{2} \frac{\Omega_{d-2}}{(2\pi)^{d-1}} (d-3) \lim_{\delta \rightarrow 0} \frac{1}{|x^0|-i\delta}\int_{\mathbb{R}^+} dl \, l^{d-4} e^{- i \sqrt{l^2+m^2}(|x^0| - i \delta)}
\end{align}
Where to reach the second line we just integrated by parts.
Setting $d=1+3$
We now set $d=4$, and use the usual $\Omega_2 = 4\pi$. Now making a change of variable $\frac{|\vec{k}|}{m}=\sinh \alpha$ we obtain:
\begin{align}
\Delta_F (x)= -\frac{1}{2} \frac{m}{2 \pi^2} \lim_{\delta \rightarrow 0}\frac{1}{|x^0|-i\delta}\int_{\mathbb{R}^+}d\alpha \cosh \alpha\,e^{- i m \cosh \alpha\,(|x^0| - i \delta)}
\end{align}
This is a Bessel function in disguise, since if the following integral exists, we have
\begin{align}
\int_{\mathbb{R}^+}d\alpha \cosh \alpha\,e^{- z \cosh \alpha} = K_1(z)
\end{align}
To ensure the integral converges, we need to close the contour of integration by using the trick of "go to the infinity, make a small angle upward and then go back to the origin".
We finally boost back to a frame where $\vec{x} \neq \vec{0}$; by Lorentz invariance of the final answer, all we need to do is modify $|x^0| \rightarrow \sqrt{(x^0)^2-|\vec{x}|^2}\equiv \sqrt{x^2}$. One finds:
\begin{equation}
\Delta_F(x)= -\frac{m}{4\pi^2} \lim_{\delta \rightarrow 0} \frac{1}{\sqrt{x^2}-i\delta} K_1(im (\sqrt{x^2}-i\delta))
\end{equation}
As described nicely here, this can then be rewritten without the $\delta$ limit using the principal-value prescription (See Weinberg QFT I, around equation 3.1.22); in doing we obtain the delta function on the light-cone mentioned in the comments.
Coefficient cross-check, via the massless limit.
We might have lost some factors of $i$, $-1$, 2, or $\pi$ in the above process, so it is always worth a cross-check. As usual, we compare our result to a limit where we know (or know how to find) the answer.
Taking the limit as $m\to 0^+$, we obtain the following
\begin{equation}
\Delta^{\text{massless}}_F(x)= \frac{i}{4\pi^2} \frac{1}{(\sqrt{x^2}-i\delta)^2} = \frac{i}{4\pi^2} \frac{1}{x^2-i\delta}
\end{equation}
which we can cross-check, say, with Peskin & Schroder equation 19.40, which calculates the position space free massless fermion propagator with our conventions, using in the process exactly the result we want to check in the process.
\begin{align}
\int \frac{d^4 k}{(2\pi)^4} \frac{i \not{k}}{k^2}e^{-i k (y-z)} &= -\not{\partial_y}\left( \int \frac{d^4 k}{(2\pi)^4} \frac{1}{k^2}e^{-i k (y-z)}\right)\\
& = -\not{\partial_y} \left(\frac{i}{4\pi} \frac{1}{(y-z)^2}\right)\\
& = -\frac{i}{2\pi^2} \frac{\gamma^\mu (y-z)_\mu}{(y-z)^4}
\end{align}
