Solving a 3D (almost radial) convolution with FFT

I have a 3D integral that is almost a radial convolution of the form $$\int d^{3}k'h(\mathbf{k'})g(|\mathbf{k-k'}|)$$ and I am looking for a fast and efficient algorithm (e.g. FFT) to solve it numerically. In the following I give you a brief description of the functions I am using and the problems I encountered.

The functions and the radial case

Consider the Fourier transform of a spherical top hat function: $$\tilde{W}_{i}(k,R_{i},R_{i+1})=\frac{1}{V_{i}}\int W_{i}(r)e^{i\mathbf{k\cdot r}}d^{3}r\\=\frac{4\pi}{V_{i}k}\int_{R_{i}}^{R_{i+1}}\sin(kr)rdr\equiv\tilde{W}_{i}(k)$$ where $R_{i}$ are some radiuses. If $f(k)$ is a general function (it comes from an interpolation, so there is no analytical formula) the convolution integral $$h_{ij}(k)=\int d^{3}k'f(|\mathbf{k'}|)\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)$$ is radial, thus it can be solved with a one-dimensional FFT algorithm.

The $k$-derivative and the convolution $d_{ij}^{(2)}(\mathbf{k})$

Considering now the $k$-derivative of one of the two window function, some terms in the integral are no longer radial: $$\int d^{3}k'f(|\mathbf{k'}|)\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial k} = \\ = \int d^{3}k'f(k')\frac{\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)}{|\mathbf{k}-\mathbf{k'}|}\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial|\mathbf{k}-\mathbf{k'}|}\left(k-\frac{\mathbf{k}\cdot\mathbf{k}'}{k}\right) \\ \equiv kd_{ij}^{(1)}(k)-\frac{\mathbf{k}}{k}\cdot d_{ij}^{(2)}(\mathbf{k}).$$ Particularly the integral $d_{ij}^{(2)}$ is still a convolution but it isn't radial: $$d_{ij}^{(2)}(\mathbf{k})=\int d^{3}k'f(|\mathbf{k'}|)\mathbf{k'}\frac{\tilde{W_{i}}\left(|\mathbf{k}-\mathbf{k'}|\right)}{|\mathbf{k}-\mathbf{k'}|}\frac{\partial\tilde{W}_{j}\left(|\mathbf{k}-\mathbf{k'}|\right)}{\partial|\mathbf{k}-\mathbf{k'}|}$$

Problems

Now, solving the integral by using common integration procedures is not really doable, because the functions $\tilde{W}_{i}(k)$ oscillate really rapidly so the numerical integral takes a lot of time and the result is never really reliable and it depends a lot on the chosen precision.

On the other hand a three dimensional FFT algorithm is really memory consuming. In fact $\tilde{W_{i}}(k)$ is a rapidly oscillating function and it needs a lot of points to be well sampled (around 5000/7000 points for each dimension) thus in three dimensions the matrices become really big. But for example I don't know anything about possible methods to parallelize the process.

I also checked other possibilities like Hankel transforms, Bessel functions for the angular part or Legendre polynomials, but without any success. If it was possible to rewrite it in the form of a 2-dimensional convolution, perhaps it would also be enough.