Let $\rho: \mathbb{R}^3\rightarrow\mathbb{R}$ be a radially symmetric function, i.e., $\rho(x,y,z)=\rho(r)$, where $r^2=x^2+y^2+z^2$. Then the Fourier transform can be computed.
$$ \begin{equation} \begin{split} \hat{\rho}(k) &=\int_{\mathbb{R}^3}\mathrm{d}x\mathrm{d}y\mathrm{d}z \exp\left(-\mathrm{i}\vec{k}\cdot\vec{r}\right)\rho(|\vec{r}|)\\ &=\int_{0}^{\pi}\mathrm{d}\theta\int_{0}^{2\pi}\mathrm{d}\varphi \int_{0}^{\infty}\mathrm{d}r\ r^2\sin(\theta)\exp\left(-\mathrm{i}kr\cos(\theta)\right)\rho(r)\\ &=\frac{4\pi}{k}\int_{0}^{\infty}\mathrm{d}r \ r\sin(kr)\rho(r)\\ \end{split} \end{equation} $$
In my case $\rho(r)$ describes a particle density.
Q: How to approximate this integral numerically if one has one array of length $N$ with $\rho$-values and another array of length $N$ with $r$-values (equidistant).** Can I relate this integral to a discrete Fourier transform? Or is there another way to compute the Fourier transform of such an array?
If it wasn't for the $0$ (instead of $-\infty$) in the lower bound of the integral, I could just compute the discrete sine transform of $r\cdot\rho(r)$ (is that correct?).
In the end I want to compute the convolution of $\rho$ with another function $w$ via
$$ (\rho*w) = \mathcal{F}^{-1}\left\{\hat{\rho}\cdot\hat{w}\right\} . $$
Indeed you can use the discrete sine transform. Just extend the integral to $-\infty$ by using $\rho(-r)=\rho(r)$, and divide by $2$. (If $\rho(x,y,z)$ is smooth, $\rho(r)$ only contains even powers of $r$, so this symmetric extension is then also smooth.)
Whether this yields a good approximation of course strongly depends on the range and resolution of the $r$ values compared to the range and structure of $\rho$.