In this paper, we first propose different combination methods to compute the Cauchy principal value integrals of oscillatory Bessel functions. By special transformations, the considered integrals are converted to finite integrals and infinite integrals. Then, the finite integrals can be calculated through the Filon-type method, the Clenshaw–Curtis–Filon method and the Clenshaw–Curtis–Filon-type method, respectively. We compute the infinite integral through the numerical steepest descent method. Moreover, the error analysis with respect to frequency ω is given through theoretical analysis. Eventually, we present several numerical experiments which are in accord with our analysis. Particularly, the accuracy can be improved by either using more nodes or adding more derivatives interpolation at endpoints. The accuracy will increase drastically with the growth of frequency ω if both the number of nodes and interpolated multiplicity are fixed.