Since the integrand has singularity on $\omega=\pm 𝑘$, so we should concern Cauchy Principal Value.
This, I believe, is where your 'mistake' lies. The Cauchy Principal Value is one way to deal with an integral that is divergent at a set of points, but it is not the way to deal with such an integral.
In this case, the Green function as-is is ill-defined because of the singular points, and we cannot perform the inverse Fourier transform that we want. If we want to make progress, we have to try something else.
One thing we can try is to look at a different integral,
$$\int_{-\infty}^\infty \frac{e^{-i \omega t}}{(\omega + i\eta - k)(\omega+i\eta + k)} \text{d} t,$$
where $\eta$ is a small positive constant. Now the poles have been shifted away from the real axis, and thus this is an integral we can actually solve using the residue theorem. In addition, in the limit $\eta \rightarrow 0$ this function superficially seems to actually approach the integral that we want! If you calculate this integral and then take the limit $\eta \rightarrow 0$ only at the end of the calculation, you will find that you obtain the retarded Green function just as in the lecture note that you link to. If you instead took the limit from below, i.e. let $\eta$ be a negative constant, you would obtain the advanced Green function.
It may seem unintuitive to have to pull these sorts of 'tricks'. To gain some intuitive understanding of the issue let us look at it from a slightly different perspective. The following is by no means meant to be rigorous, but rather an attempt at an intuitive explanation.
The Green function $G(\omega)$ basically tells you what the response of the system is when the 'driving term' is of the form $e^{i \omega t}$.
Let us consider a harmonic oscillator described by the equation
$$ \ddot{x} + \omega_0^2 x = F,$$
where $F$ is a driving force. If $F$ is a harmonic function, $F = F_0 e^{i \omega t}$ then the response will be
$$x(t) = \frac{F_0 e^{i\omega t}}{\omega_0^2 - \omega^2},$$
however this solution is clearly not valid if $\omega = \omega_0$. In that case, we will instead find that $x$ grows without bound as time increases. Consider how this changes if we add a little bit of friction to the system, though. Then the equation of motion becomes
$$\ddot{x} + \eta \dot{x} + \omega_0^2 x = F,$$
and the response to a harmonic driving term $F = F_0 e^{i \omega t}$ becomes
$$x(t) = \frac{F_0 e^{i\omega t}}{\omega_0^2 - \omega^2 + i\omega\eta},$$
which is well-defined for all real $\omega$ and, if we take the limit $\eta \rightarrow 0$, gives the same answer as in the friction-less sytem, except when the frequency of the driving term is at the resonance frequency $\omega_0$.
The problem with the friction-less harmonic oscillator is not really that it is friction-less, but rather the driving term; it is of infinite duration, extending both into the infinite past and the infinite future. Such a driving term is not very physical at all! If you restrict yourself to physical driving terms, e.g. terms that have a finite duration and never diverge, then you wouldn't encounter problems like these. But for such a physical driving term, observed over a finite duration of time, there would not be a difference between the response of the friction-less system, and the system with friction in the limit $\eta \rightarrow 0$. Therefore, we can allow ourselves to replace the "pure", friction-less system that we wanted to describe, by a different system that has a negligible but in principle non-zero friction, if we at the same time restrict ourselves to working only with 'physical' driving terms.
The introduction of an $\eta$ when calculating a retarded (or advanced) Green function can be understood intuitively in a similar manner.
Best Answer
The first step is to split up the integrand using partial fractions, $$\int dz\left(\frac{1}{z-a}\frac{1}{z-x}\right)=\int dz\left(\frac{1}{x-a}\frac{1}{z-x}-\frac{1}{x-a}\frac{1}{z-a}\right).$$ Then the $1/(x-a)$ can be pulled out of the $z$-integration, $$\int dz\left(\frac{1}{z-a}\frac{1}{z-x}\right)=\frac{1}{x-a}\int dz\left(\frac{1}{z-x}-\frac{1}{z-a}\right).$$ Taking the principal value, the first term in the remaining integrand gives zero, since it has just a simple pole on the real line contour of integration, $$\mathcal{P}\int_{-\infty}^{+\infty}\frac{dz}{z-x}=\lim_{\epsilon\rightarrow0}\left( \int_{-\infty}^{x-\epsilon}\frac{dz}{z-x}+\int_{x+\epsilon}^{+\infty}\frac{dz}{z-x}\right)\\=\lim_{\epsilon\rightarrow0}\left[\lim_{t\rightarrow+\infty}\left( \int_{-t+x}^{x-\epsilon}\frac{dz}{z-x}+\int_{x+\epsilon}^{t+x}\frac{dz}{z-x}\right)\right]=0.$$ (This is really just a shifted version of $\mathcal{P}\int dz/z=0$, the basic identity for principal value integrals.) With the first term in the integrand taken care of, that just leaves the second term, which has no singularity on the real axis, so we can evaluate it directly as an improper integral (no principal value required), $$-\int_{-\infty}^{+\infty}\frac{dz}{z-a}=-\lim_{t\rightarrow+\infty}\int_{-t}^{t}\frac{dz}{z-a}=\lim_{t\rightarrow+\infty}\left[-\ln\left(\frac{t-a}{-t-a}\right)\right].$$ As $t\rightarrow+\infty$, the complex number $t-a$, written in polar form, approaches $te^{i0}$, while $-t-a\,(\approx-|t|)$ approaches $te^{i\pi}$. (Because the contour of integration passes above the pole at $a$, the argument of the complex number is decreasing from the lower limit at $-\infty$ to the upper limit at $+\infty$.) Therefore the logarithm reduces to $$\lim_{t\rightarrow+\infty}\left[-\ln\left(\frac{t-a}{-t-a}\right)\right]=-\ln\frac{e^{i0}}{e^{i\pi}}=i\pi.$$ Combining that with the overall factor of $1/(x-a)$ that was pulled out of the integrations, we have the Hilbert transform $$H\frac{1}{x-a}=\frac{1}{\pi}\mathcal{P}\int_{-\infty}^{+\infty} dz\left(\frac{1}{z-a}\frac{1}{z-x}\right)=\frac{i}{x-a}.$$