Away from the singular point the standard quadrature methods should be fine, so the problem is really about computing an integral of the form
$$ \mathrm{P.V.}\int_{-h}^{h} \frac{f(x)}{x}\,dx. $$
I would suggest trying the following things:
If you know a closed-form Taylor series for $f(x)=G(\omega+x)$, the you can choose $h$ to small enough for the series approximation to be accurate, and calculate the principal value in closed form for each series term.
Alternatively, you can try to interpolate $f$ numerically, using something like polyfit, which will give you the Taylor series numerically. Depending on how well-behaved your function is, it is difficult to say from the start whether any of this will definitely work.
Still the simplest way (although numerically unstable) is to write the integral as
$$ \int_0^h \frac{f(x)-f(-x)}{x}\,dx, $$
and apply the standard quadrature algorithm to this. The issue is that the calculation $f(x)-f(-x)$ will lead to inaccuracies on the order of $\int_0^h (\epsilon_{\mathrm{mach}}/x)\,dx$, which may or may not be tolerable. If you have a way to simplify the expression $f(x)-f(-x)$ to make it numerically stable, this is probably the best option.
There are ways to derive quadrature nodes and weights for computing principal values of such integrals (just search for them, even GSL has them), but I couldn't find this in matlab.
These are numerical methods for principal value integrals on the real line. You can also try to compute the integral as
$$ \mathrm{P.V.}\int_{-h}^{h}\frac{f(x)}{x}\,dx = \int_\gamma \frac{f(x)}{x}\,dx + \pi \mathrm{i}f(0), $$
where $\gamma$ is a contour going from $-h$ to $h$ passing above the real line. If you can evaluate $G$ for complex $\omega$, then even standard quadrature routines will be able to compute the integral over $\gamma$.
