I want to numerically evaluate a convolution that contains the function

$$f(t) = \dfrac{1}{e^{t}-1},$$


$$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.

  • $\begingroup$ The functions involved must be Fourier transformable for you to be able to use an FFT technique. Have you checked if it is? $\endgroup$ – mathreadler Mar 21 '17 at 10:53
  • 1
    $\begingroup$ I just tried putting $1/(e^{t}-1)$ into WolframAlpha and it does give me a solution. $\endgroup$ – Medulla Oblongata Mar 21 '17 at 11:22
  • 1
    $\begingroup$ You can try build the grid so the singularity ends up in the middle of one discrete sample. Then use the Cauchy principal value as an approximation for that sample: $\int_{-\Delta_T/2}^{\Delta_T/2}f(t)dt$, $\Delta_T$ being your time length for a sample. All the other samples you could just sample $f(t)$ by direct computation. $\endgroup$ – mathreadler Mar 21 '17 at 11:55
  • 1
    $\begingroup$ Very fast temporally located events like singularities are usually inconvenient with FFT because they spread all over the spectrum. The Fourier basis functions which you build linear combinations of are global but the singularity is very localized. Forces us to use very many non-zero coefficients. $\endgroup$ – mathreadler Mar 21 '17 at 12:05
  • $\begingroup$ That sounds interesting. Do you mind reposting your comments (with possibly more detail) as an answer below? Thanks $\endgroup$ – Medulla Oblongata Mar 21 '17 at 12:13

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.

  • $\begingroup$ Do you know of any good references on signal processing and that cover these issues involving singularities? I'd like to read more into this. $\endgroup$ – Medulla Oblongata Mar 21 '17 at 12:45
  • $\begingroup$ There are very many books on signal processing, it's a bit difficult to say any book in particular without knowing a little about the background of the student, content of previous courses and things like it. $\endgroup$ – mathreadler Mar 21 '17 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.