I want to numerically evaluate a convolution that contains the function
$$f(t) = \dfrac{1}{e^{t}-1},$$
like
$$h(t) = f(t) \circ g(t)$$
where $g(t)$ is a smooth well-behaved function that converges to zero at $\pm \infty$.
I know that to perform the integration $\int_{-a}^{a}f(t)dt$ we can simply use the Cauchy principal value. But what is the best way to solve a convolution with a singularity? I want to be able to calculate convolutions numerically in Matlab using a FFT-based technique.
You can build the discretized approximation of your function $$f(t) = \frac{1}{e^t-1}$$
When discretizing the function if using time steps of length $\Delta_T$ : $\{\cdots, -2\Delta_T,-\Delta_T,0,\Delta_T, 2\Delta_T, \cdots\}$. For many functions a trivial way to sample would be to just calculate $f(t)$ for the value directly. For our function which has a singularity this has high likelihood of introducing error wherever the function has too violent behaviour inside of one sample. Instead we can choose to sample with short term integrals:
$$f[k] = \frac{1}{\Delta_T}\int_{\Delta_T(k-1/2)}^{\Delta_T(k+1/2)}f(t)dt$$
This way we can use the Cauchy principal value as replacement for this integral for the one sample that hits right on the singularity: $t\approx 0, (k=0)$. We can choose to also do it for other samples, but since the function is so smooth everywhere else it is likely that the central Riemann sum is quite close to the short time integral value:
$$\Delta_T \cdot f(k\Delta_T ) \approx \int_{\Delta_T(k-1/2)}^{\Delta_T(k+1/2)}f(t)dt$$
This is often depicted in calculus books with middle point of functions on intervals times interval width being almost equal to integral (area under curve).
When considering FFT for calculating convolutions one should bear in mind an important property of the basis functions (complex exponentials): they are all global (no compact support). Using linear combinations of such functions to describe something with a very localized behavior (like singularities are) can be a bad idea = require many non-zero values. And we may also lose precision.