1
$\begingroup$

I want to obtain the convolution of two discretized real functions $f$ and $g$, $$ c(t) = \int_{-\infty}^{+\infty} \mathrm{d}{x} \, f(x) \, g(t-x) \tag{1} $$ via discrete Fourier transform (DFT).

As a concrete example, let $x \in [-6, 6]$ be the $x$-axis on which the functions are defined, with $$ f(x) = \Theta(x) \, e^{-x} ~, \\ g(x) = \Pi(x/2) ~, $$ where $\Theta(x)$ is the Heaviside step-function and $\Pi(x)$ is the rectangular function. Then the exact analytical expression for the convolution will be (see also figure below) $$ C(t) = \begin{cases} 0 , & t \leq -1 \\ 1 - e^{t + 1} , & -1 < t \leq 1 \\ e^{-t} (e - e^{-1}) , & t \geq 1 \end{cases} $$ exact convolution of f and g

The $x$-axis is discretized into $N = 2^K$ points. DFT (and inverse DFT) can be used to obtain the periodic (circular) convolution $c_P$ of $f$ and $g$, $$ c_P[n] = DFT^{-1} ( DFT(f) \cdot DFT(g) ) ~. $$

According to what I have seen in text-books (e.g. section 2.4 of this), the aperiodic (‘infinite’) convolution (corresponding to infinite ‘signals’ for $f$ and $g$) can be obtained by zero-padding $f$ and $g$ arrays at the end by $N-1$ zeros. Then the discretized convolution corresponding to Eq. 1 should be obtained as $$ c[n] = DFT^{-1} ( DFT(f_{padded}) \cdot DFT(g_{padded}) ) $$

I tried to do this in Python using scipy.fftpack (see the code here). But the results of this discrete convolution does not match with the analytic one (see figure); namely, the scale and the position of the discrete convolution is not correct.

exact vs. discrete convolution of f and g

Is this the correct way to perform the discrete (infinite) convolution of $f$ and $g$?


*Note: The correct result is actually obtained by padding $f$ by $N-1$ zeros at the end, re-arranging (center-padding) $g$ as $\hat{g} := (g_{+} ~|~ N-1 ~ \text{zeros} ~|~ g_{-})$, where $g_{\pm}$ denotes the right/left half of $g$, and applying DFT as $$ c[n] = DFT^{-1} ( DFT(f_{padded}) \cdot DFT(\hat{g}) ) ~; $$ finally one should divide elements of $c$ by $\sum_{i} g_i$. But I cannot understand the reason behind this.

$\endgroup$
  • $\begingroup$ What DFT are you using? You probably forgot to divide by $N$ or $\sqrt N$ when you took the inverse DFT. $\endgroup$ – Rodrigo de Azevedo Jun 12 '16 at 9:17
  • 1
    $\begingroup$ You'll also need to scale for the discrete time step $\Delta x$. $\endgroup$ – Brian Borchers Jun 12 '16 at 13:32
2
$\begingroup$

You want to calculate $$c(t)=\int_{-6}^6dx f(x)g(t-x)$$. Your functions are not periodic, but you pad them with zeros so that assuming them to be repeated periodically does not mess up the calculation of the convolution. You decided to pad them with zeros from $6$ to $18$. So you now want to calculate $$c(t)=\int_{-6}^{18}dx f(x)g(t-x)$$ by discrete Fourier transform on the interval from $-6$ to $18$. This means that you need to sample your functions at $N$ equally-spaced points $x_j=j\Delta x - 6$, $j=0,\dots N-1$, where $\Delta x=24/N$. Note the shift by $-6$ in the definition of the $x_j$ that is necessary because your interval starts at $-6$. We can now approximate the integral by a sum: $$c(t)\approx\Delta x\sum_{j=0}^{N-1}f(x_j)g(t-x_j).$$ Notice the factor of $\Delta x$ that you forgot in your calculation and that explains the wrong amplitude you found. Next we consider only a discrete set of $t$s so that $t-x_j$ is also one of the discrete points at which you have sampled your functions. It is convenient to set $t=x_n-6$ for some $n$, because then $t-x_j=x_n-6-x_j=x_{n-j}$ and the continuous convolution you want to calculate is approximated by a discrete convolution of the form that can be calculated by discrete Fourier transform $$c(x_n-6)\approx\Delta x\sum_{j=0}^{N-1}f(x_j)g(x_{n-j}).$$ It is this shift by $-6$ that explains why the graphs do not line up in your figure.

$\endgroup$

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.