While chaohuang's answer gave you basically everything you needed, I'd like to elaborate a little on the last part, and on generalized functions (a.k.a. distributions) in general; I apologize for the unavoidable pun.
I realize that in physics, you might treat generalized functions, like the delta distribution, in a "classical" way similar to regular functions. The truth is, however, that distributions are really only defined in the context of how they act on a set of (nicely behaved) test functions. These test functions are usually either Schwartz space functions, or compactly supported functions. Occasionally we also consider the (much larger) test space of $C^\infty$ functions. So the reason why $\hat{H}(\xi)$ involves $\text{PV}\left(\dfrac{1}{\xi}\right)$, as opposed to just $\dfrac{1}{\xi}$, is because of how it integrates against a test function in Schwartz space.
To make
\begin{align}
\hat{H}(\xi) &= \int_{-\infty}^{\infty}H(x)e^{-ix\xi}\;dx\\
&= \lim_{\epsilon\to0^+}\;\int_{0}^{\infty}e^{-\epsilon x}e^{-ix\xi}\;dx\\
&= \lim_{\epsilon\to0^+}\;(\epsilon+i\xi)^{-1}\\
&= \lim_{\epsilon\to0^+}\;-i(\xi-i\epsilon)^{-1}\\
&= -i(\xi-i0)^{-1}
\end{align}
rigorous, you need to see how it acts as a linear functional on Schwartz space:
\begin{align}
\hat{H}[\varphi] &= \int_{-\infty}^{\infty}\hat{H}(\xi)\varphi(\xi)\;d\xi\\
&= \lim_{\epsilon\to0^+}\;\int_{-\infty}^{\infty}-i(\xi-i\epsilon)^{-1}\varphi(\xi)\;d\xi\\
&= -i\lim_{\epsilon\to0^+}\;\int_{-\infty}^{\infty}\left(\frac{\xi+i\epsilon}{\xi^2+\epsilon^2}\right)\varphi(\xi)\;d\xi\\
&= -i\left(\lim_{\epsilon\to0^+}\;\int_{-\infty}^{\infty}\left(\frac{\xi}{\xi^2+\epsilon^2}\right)\varphi(\xi)\;d\xi + i\lim_{\epsilon\to0^+}\;\int_{-\infty}^{\infty}\left(\frac{\epsilon}{\xi^2+\epsilon^2}\right)\varphi(\xi)\;d\xi\right)\\
&=-i\lim_{\epsilon\to0^+}A(\epsilon)+\lim_{\epsilon\to0^+}B(\epsilon)
\end{align}
where
\begin{align}
A(\epsilon) &= \int_{-\infty}^{\infty}\left(\frac{\xi}{\xi^2+\epsilon^2}\right)\varphi(\xi)\;d\xi\\
&= \int_{-\infty}^{\infty}\frac{d}{d\xi}\left(\frac{1}{2}\ln(\xi^2+\epsilon^2)\right)\varphi(\xi)\;d\xi\\
&= -\frac{1}{2}\int_{-\infty}^{\infty}\ln(\xi^2+\epsilon^2)\varphi'(\xi)\;d\xi\\
\end{align}
and
\begin{align}
B(\epsilon) &= \int_{-\infty}^{\infty}\left(\frac{\epsilon}{\xi^2+\epsilon^2}\right)\varphi(\xi)\;d\xi\\
&= \int_{-\infty}^{\infty}\frac{d}{d\xi}\left(\arctan\left(\frac{\xi}{\epsilon}\right)\right)\varphi(\xi)\;d\xi\\
&= -\int_{-\infty}^{\infty}\arctan\left(\frac{\xi}{\epsilon}\right)\varphi'(\xi)\;d\xi\\
\end{align}
with
$$\lim_{\epsilon\to0^+}\;A(\epsilon) = -\int_{-\infty}^{\infty}\ln(|\xi|)\varphi'(\xi)\;d\xi$$
and
$$\lim_{\epsilon\to0^+}\;B(\epsilon) = -\frac{\pi}{2}\int_{-\infty}^{\infty}\text{sgn}(\xi) \varphi'(\xi)\;d\xi$$
by various Lebesgue integration theorems, and the fact that $\varphi$ and its derivatives are rapidly decaying. In particular, whereas the area under $1/\xi$ is infinite to either side of the singularity (making the integration of $\varphi/\xi$ ill-defined around $\xi=0$), the function $\ln|\xi|$ has finite integral near $\xi=0$, and so $\lim_{\epsilon\to0^+}\;A(\epsilon)$ is well-defined.
An easy computation then gives us that
$$\lim_{\epsilon\to0^+}\;A(\epsilon) = \lim_{\epsilon\to0^+}\; \int_{\mathbb{R}\setminus(-\epsilon,\epsilon)}\frac{1}{\xi}\varphi(\xi)\;d\xi =: \text{PV}\left(\frac{1}{\xi}\right)[\varphi]$$
by definition of the Cauchy principal value distribution. Similarly, we find that
$$\lim_{\epsilon\to0^+}\;B(\epsilon) = 2\left(\frac{\pi}{2}\varphi(0)\right):= \pi\delta[\varphi]$$
by definition of the delta distribution.
So in conclusion, we have finally arrived at the fact that
$$\hat{H}(\xi) = -i\text{PV}\left(\frac{1}{\xi}\right) + \pi\delta$$
Another solution
The distribution $\mathrm{pv} \frac{1}{x}$ satisfies $x \, \mathrm{pv} \frac{1}{x} = 1.$ Therefore,
$$
2\pi \, \delta(\xi)
= \mathcal{F} \{ 1 \}
= \mathcal{F} \{ x \, \mathrm{pv} \frac{1}{x} \}
= i \frac{d}{d\xi} \mathcal{F} \{ \mathrm{pv} \frac{1}{x} \}
$$
Thus,
$
\mathcal{F} \{ \mathrm{pv} \frac{1}{x} \}
= -i \pi \, \operatorname{sign}(\xi) + C
$
for some constant $C$.
But $\mathrm{pv} \frac{1}{x}$ is odd so its Fourier transform must also be odd, and since $-i \pi \, \operatorname{sign}(\xi)$ is odd while $C$ is even, we must have $C=0.$
Best Answer
We can't use Fubini directly as noted in the OP. Let $S_{\varepsilon,R}:=\{t\in\Bbb R, \varepsilon<|t|<R\}$. Using Fubini's theorem and a rewriting of the inner integral, we come up with the equality $$\int_{\Bbb R\times\Bbb R}\frac{e^{-isx}}x\varphi(s)\chi_{S_{\varepsilon,R}}(x)dxds=-\int_{\Bbb R}\varphi(s)\int_{s\varepsilon}^{sR}\frac{\sin u}ududs.$$ We can find $M$ such that for all $t$, $\left|\int_0^t\frac{\sin u}u du\right|<M$. This allows us to use the dominated convergence theorem in order to take the limit $R\to +\infty$ in the displayed equality. This gives $$\int_{\{|x|>\varepsilon\}}\int_{\Bbb R}\frac{e^{-isx}}x\varphi(s)dxds=-\int_{\Bbb R}\varphi(s)\int_{s\varepsilon}^{+\infty}\frac{\sin u}ududs.$$ Then use dominated convergence theorem again to take the limit with respect to $\varepsilon$.