14
$\begingroup$

I am working in digital image restoration field. I have read all things about convolution, that, for an LTI system, if we know its impulse response, then we can find its output by just using convolution between input and impulse response.

Can anyone tell me that what is the main mathematical philosophy behind it? Your experience with will tell more to me than just internet surfing about it.

$\endgroup$
13
$\begingroup$

The Idea of Convolution

My favorite exposition of the topic is in one of Brad Osgood's lectures on the Fourier Transform. The discussion of convolution begins around 36:00, but the whole lecture has additional context that's worth watching.

The basic idea is that, when you define something like the Fourier Transform, rather than working directly with the definition all the time, it's useful to derive higher-level properties that simplify calculations. For example, one such property is that the transform of the sum of two functions is equal to the sum of the transforms, i.e.

$$ F\{f + g\} = F\{f\} + F\{g\}. $$

That means if you have a function with an unknown transform, and it can be decomposed as a sum of functions with known transforms, you basically get the answer for free.

Now, since we have an identity for the sum of two transforms, it's a natural question to ask what the identity is for the product of two transforms, i.e.

$$ F\{f\}F\{g\} = \space ?. $$

It turns out that when you calculate the answer, convolution is what appears. The entire derivation is given in the video, and since your question is mostly conceptual, I won't recapitulate it here.

The implication of approaching convolution in this way is that it's an intrinsic part of the way the Laplace Transform (of which the Fourier Tranform is a special case) turns linear constant-coefficient ordinary differential equations (LCCODE's) into algebraic equations. That fact that such a transform is available to make LCCODE's analytically tractable is a large part of the reason why they're studied in signal processing. For example, to quote Oppenheim and Schafer:

Because they are relatively easy to characterize mathematically and because they can be designed to perform useful signal processing functions, the class of linear shift-invariant systems will be studied extensively.

So one answer to the question, is that if you're using transform methods to analyze and/or synthesize LTI systems, sooner or later, convolution will arise (either implicitly or explicitly). Note that this approach to introducing convolution is very standard in the context of differential equations. For example, see this MIT lecture by Arthur Mattuck. Most presentations either present the convolution integral without comment, then derive its properties (i.e. pull it out of a hat), or hem and haw about the strange form of the integral, talk about flipping and dragging, time-reversal, etc., etc.

The reason I like Prof. Osgood's approach is that it avoids all that tsouris, as well as providing, in my opinion, deep insight into how mathematicians probably arrived at the idea in the first place. And I quote:

I said, "Is there a way of combining F and G in the time domain, so that in the frequency domain the spectrums multiply, the Fourier transforms multiply?" And the answer is, yes, there is, by this complicated integral. It's not so obvious. You wouldn't get out of bed in the morning and write this down, and expect that this was going to solve that problem. How do we get it? You said, suppose the problem is solved, see what has to happen, and then we'd have to recognize when it's time to declare victory. And it's time to declare victory.

Now, being an obnoxious mathematician, you cover your tracks and say, "Well, I'm simply going to define the convolution of two functions by this formula."

LTI Systems

In most DSP texts, convolution is usually introduced in a different way (that avoids any reference to transform methods). By expressing an arbitrary input signal $x(n)$ as a sum of scaled and shifted unit impulses,

$$ x(n) = \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}, \tag{1} $$

where

$$ \delta(n) = \begin{cases} 0, & n \ne 0 \\ 1, & n = 0, \end{cases} \tag{2} $$

the defining properties of linear time-invariant systems lead directly a convolution sum involving the impulse response $h(n) = L[\space \delta(n) \space]$. If the system defined by an LTI operator $L$ is expressed as $y(n) = L[\space x(n)\space]$, then by applying the repective properties, namely linearity

$$ \overbrace{L[\space ax_1(n) + bx_2(n)\space ]}^\text{Transform of the sum of scaled inputs} = \overbrace{aL[\space x_1(n)\space ] + bL[\space x_2(n)\space ]}^\text{Sum of scaled transforms} , \tag{3} $$

and time/shift invariance

$$ L[\space x(n) \space] = y(n) \space \xrightarrow{\text{implies}} L[\space x(n-k) \space] = y(n-k), \tag{4} $$

the system can be rewritten as

$$ y(n) = \overbrace{L\left[ \sum_{k=-\infty}^\infty{x(k)\delta(n - k)}\right]}^\text{Tranform of the sum of scaled inputs} = \overbrace{\sum_{k=-\infty}^\infty{x(k)L[\delta(n - k)]}}^\text{Sum of scaled transforms} = \overbrace{\sum_{k=-\infty}^\infty{x(k)h(n-k) }. }^\text{Convolution with the impulse response} $$

That's a very standard way to present convolution, and it's a perfectly elegant and useful way to go about it. Similar derivations can be found in Oppenheim and Schafer, Proakis and Manolakis, Rabiner and Gold, and I'm sure many others. Some deeper insight [that goes further than the standard introductions] is given by Dilip in his excellent answer here.

Note, however, that this derivation is somewhat of a magic trick. Taking another look at how the signal is decomposed in $\left(1\right)$, we can see that it's already in the form of a convolution. If

$$ \overbrace{(f * g)(n)}^\text{f convolved with g} = \sum_{k=-\infty}^\infty{f(k)g(n-k)}, $$

then $\left(1\right)$ is just $x * \delta$. Because the delta function is the identity element for convolution, saying any signal can be expressed in that form is a lot like saying any number $n$ can be expressed as $n + 0$ or $n \times 1$. Now, choosing to describe signals that way is brilliant because it leads directly to the idea of an impulse response--it's just that the idea of convolution is already "baked in" to the decomposition of the signal.

From this perspective, convolution is intrinsically related to the idea of a delta function (i.e. it's a binary operation that has the delta function as its identity element). Even without considering its relation to convolution, the description of the signal depends crucially on the idea of the delta function. So the question then becomes, where did we get the idea for the delta function in the first place? As far as I can tell, it goes at least as far back as Fourier's paper on the Analytical Theory of Heat, where it appears implicitly. One source for further information is this paper on Origin and History of Convolution by Alejandro Domínguez.

Now, those are the two of the main approaches to the idea in the context of linear systems theory. One favors analytical insight, and the other favors numerical solution. I think both are useful for a full picture of the importance of convolution. However, in the discrete case, neglecting linear systems entirely, there's a sense in which convolution is a much older idea.

Polynomial Multiplication

One good presentation of the idea that discrete convolution is just polynomial multiplication is given by Gilbert Strang in this lecture starting around 5:46. From that perspective, the idea goes all the way back to the introduction of positional number systems (which represent numbers implicitly as polynomials). Because the Z-transform represents signals as polynomials in z, convolution will arise in that context as well--even if the Z-transform is defined formally as a delay operator without recourse to complex analysis and/or as a special case of the Laplace Transform.

$\endgroup$
  • $\begingroup$ thanks sir for your invaluable guidance, you have just shown me the right way to follow. Your this help has taught me that how to be a good human begin for others.... :) $\endgroup$ – Mayank Tiwari Jun 27 '13 at 10:11
  • $\begingroup$ How does this great coincidence explain that you need to do the convolution in he case asked? I believe that in every domain, there is an operation that turns into convolution when you revert the arguments into the time domain. Might be we need to do muiltiplication in the time domain to get the response? Why should we multiply the frequencies instead of time sweeps? $\endgroup$ – Val Jun 27 '13 at 11:36
  • 1
    $\begingroup$ Considering that the OP had already asked a question about the role of impulses in relation to LTI systems, I read the question as him using that as an example to motivate a question about where convolution comes from--not necessarily its role in calculating an impulse response per se. Is that what you're asking? $\endgroup$ – datageist Jun 27 '13 at 11:52
  • 1
    $\begingroup$ Saying that we need convolution because it is identical to fourier multiplication sounds nonsense to me in case we do not know why we need the fourier multiplication. When we are given the impulse response, this means time domain and convolution rather than any black magic in fourier basis. I do not think that reference to that question may clarify this. In any case, it is not good to give "localized answers" to general, fundamental (i.e phylosophical) questions. Q&A must be useful for future visitors. $\endgroup$ – Val Jun 27 '13 at 12:12
  • $\begingroup$ Val's comment above is right on target. I wonder how linear systems worked before Fourier transforms were invented/discovered. How on earth did a non-sentient inanimate object discover such a complicated formula? $\endgroup$ – Dilip Sarwate Jun 27 '13 at 13:33
6
$\begingroup$

I once gave the answer in Wikipedia convolution discussion page, which asked basically the same question: Why the time inversion?. The philosophy is that you apply a single pulse at time 0 to your filter and record its response at time 0,1,2,3,4,…. Basically, response will look like a function, h(t). You may plot it. If pulse was n times taller/higher, response pulses will be proportionally taller (this is because linear filter is always assumed). Now, all DSP (and not only DSP) is about what happens when you apply the filter to your signal? You know the impulse response. Your signal (especially digital) is nothing more than a series of pulses of height x(t). It has height/value $x$ at time $t$. Linear systems are cool that you can sum the outputs for each such input pulse to get the response function y(t) for input function x(t). You know that output pulse y(t=10) depends on immediate input x(10), which contributes h(0)*x(10). But there is also contribution, x(9)*h(1), to the output from the previous pulse, x(9), and contributions from even earlier input values. You see, as you add contributions from earlier inputs, one time argument decreases while another increases. You MAC all the contributions into y(10) = h(0)*x(10) + h(1)*x(9) + h(2)*x(8) + …, which is a convolution.

You can think of functions y(t), h(t) and x(t) as vectors. Matrices are operators in the linear algebra. They take input vector (a series of numbers) and produce output vector (another series of numbers). In this case, the y is a product of convolution matrix with vector x,

$\vec y = \begin{bmatrix}y_0\\y_1\\y_2 \\ \vdots \end{bmatrix} = \begin{bmatrix}h_0 & 0 & 0 & \cdots \\ h_1 & h_0 & 0 & \cdots \\ h_2 & h_1 & h_0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix}x_0\\x_1\\x_2 \\ \vdots \end{bmatrix} = H \vec x$

Now, because convolution is a Toeplitz matrix, it has Fourier eigenbasis and, therefore, convolution operator (linear operators are represented by matrices, but matrix also depends on the basis) is a nice diagonal matrix in the fourier domain,

$\vec Y = \begin{bmatrix}Y_0\\Y_1\\Y_2 \\ \vdots \end{bmatrix} = \begin{bmatrix} \lambda_0 & 0 & 0 & \cdots \\ 0 & \lambda_1 & 0 & \cdots \\0 & 0 & \lambda_2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix}X_0\\X_1\\X_2 \\ \vdots \end{bmatrix} = diag(H)\vec X$

Note, much more zeroes and, thus, much simpler computation. This result is know as "convolution theorem" and, as first response answered, it is much simpler in the fourier domain. But, this is phylosophy behind the "convolution theorem", fourier basis and linear operators rather than ubiquitous need for convolution.

Normally, you do convolution because you have your input signal, impulse response and require output in time domain. You may transform into fourier space to optimize computation. But, it is not practical for simple filters, as I've seen in the DSPGuide. If your filter looks like $y[currentTime] = k_1 x[time-1] + k_2 x(time-2) + b * y[time-1]$, it makes no sense to fourier transform. You just do n multiplications, for computing every y. It is also natural for real-time. In real-time you compute only one y at a time. You may think of Fourier transform if you have your signal x recorded and you need to compute the whole vector y at once. This would need NxN MAC operations and Fourier can help to reduce them to N log(N).

$\endgroup$
  • $\begingroup$ A couple notes: how would you extend this description for the continuous-time case (which obviously came before the discrete-time case)? Also, there are many real-time applications that use Fourier-transform-based methods for fast convolution. Saying that outputs are always calculated one at a time for real-time applications just isn't true. $\endgroup$ – Jason R Jun 27 '13 at 13:04
  • $\begingroup$ With that said, nice job pointing out the fact that the Toeplitz structure of the convolution matrix implies that it admits a diagonal representation under a Fourier basis. $\endgroup$ – Jason R Jun 27 '13 at 13:06
  • $\begingroup$ Yes, may be you use fourier transfrom in the real-time. I am far from being DSP expert. I just expressed the idea (which I got from my scarce practice and reading DSPGuide). Anyway, I want to highlight that fourier has nothing to do with the phylosophy of convolution. Might be I need to remove all the fourier-related discussion, since it is distracting. Convolution is natural in the time domain and is needed without any Fourier, no matter how cool the Fourier is. $\endgroup$ – Val Jun 27 '13 at 13:14
  • $\begingroup$ The fact that continous-time case came before historically, does not demand us that we should learn in the same order. I think that it is easier to understand the philosophy of many things, including convolution, starting with the discrete case (and we are in DSP.se to take this approach). In continous-case a series of MAC operations is turned into integration, as we make our pulses shorter and shorter. BTW, integration itself $\int f(x)dx$ is understood as a limit case of the discrete summation: $\int f(x)dx = \lim_{dx \to 0} (\sum f(x) dx) $. So, it cannot come before discrete summation. $\endgroup$ – Val Jun 27 '13 at 13:47
  • $\begingroup$ @JasonR In the continuous setting, you'd replace the Toeplitz matrix by a shift-invariant operator. You then can show that the Fourier basis functions diagonalize this operator. $\endgroup$ – lp251 Jun 28 '13 at 3:30
3
$\begingroup$

Although the previous answers have been really good, I would like to add my viewpoint about convolution where I just make it easier to visualize due to the figures.

One wonders if there is any method through which an output signal of a system can be determined for a given input signal. Convolution is the answer to that question, provided that the system is linear and time-invariant (LTI).

Assume that we have an arbitrary signal $s[n]$. Then, $s[n]$ can be decomposed into a scaled sum of shifted unit impulses through the following reasoning. Multiply $s[n]$ with a unit impulse shifted by $m$ samples as $\delta[n-m]$. Since $\delta[n-m]$ is equal to 0 everywhere except at $n=m$, this would multiply all values of $s[n]$ by 0 when $n$ is not equal to $m$ and by 1 when $n$ is equal to $m$. So the resulting sequence will have an impulse at $n=m$ with its value equal to $s[m]$. This process is clearly illustrated in Figure below.

$\hspace{2cm}$enter image description here

This can be mathematically written as \begin{equation*} s[n] \delta [n-m] = s[m] \delta [n-m] \end{equation*} Repeating the same procedure with a different delay $m'$ gives \begin{equation*} s[n] \delta [n-m'] = s[m'] \delta [n-m'] \end{equation*}

The value $s[m']$ is extracted at this instant. Therefore, if this multiplication is repeated over all possible delays $-\infty < m < \infty$, and all produced signals are summed together, the result will be the sequence $s[n]$ itself. \begin{align} s[n]&= \cdots + s[-2]\delta[n+2] + s[-1]\delta[n+1] + \nonumber \\ &\quad\quad\quad\quad\quad s[0]\delta[n] + s[1]\delta[n-1] + s[2]\delta[n-2] + \cdots \nonumber \\ &= \sum \limits _{m = -\infty} ^{\infty} s[m] \delta[n-m] \label{eqIntroductionSeqImpulses} \end{align}

In summary, the above equation states that $s[n]$ can be written as a summation of scaled unit impulses, where each unit impulse $\delta[n-m]$ has an amplitude $s[m]$. An example of such a summation is shown in Figure below.

$\hspace{3cm}$enter image description here

Consider what happens when it is given as an input to an LTI system with an impulse response $h[n]$.

enter image description here

This leads to an input-output sequence as

enter image description here

During the above procedure, we have worked out the famous convolution equation that describes the output $r[n]$ for an input $s[n]$ to an LTI system with impulse response $h[n]$.

Convolution is a very logical and simple process but many DSP learners can find it confusing due to the way it is explained. We will describe a conventional method and another more intuitive approach.

Conventional Method


Most textbooks after defining the convolution equation suggest its implementation through the following steps. For every individual time shift $n$,

[Flip] Arranging the equation as $r[n] = \sum _{m = -\infty} ^{\infty} s[m] h[-m+n]$, consider the impulse response as a function of variable $m$, flip $h[m]$ about $m = 0$ to obtain $h[-m]$.

[Shift] To obtain $h[-m+n]$ for time shift $n$, shift $h[-m]$ by $n$ units to the right for positive $n$ and left for negative $n$.

[Multiply] Point-wise multiply the sequence $s[m]$ by sequence $h[-m+n]$ to obtain a product sequence $s[m]\cdot h[-m+n]$.

[Sum] Sum all the values of the above product sequence to obtain the convolution output at time $n$.

[Repeat] Repeat the above steps for every possible value of $n$.

An example of convolution between two signals $s[n] = [2\hspace{1mm}-\hspace{-1mm}1\hspace{2mm} 1]$ and $h[n] = [-1\hspace{2mm} 1\hspace{2mm} 2]$ is shown in Figure below, where the result $r[n]$ is shown for each $n$.

Note a change in signal representation above. The actual signals $s[n]$ and $h[n]$ are a function of time index $n$ but the convolution equation denotes both of these signals with time index $m$. On the other hand, $n$ is used to represent the time shift of $h[-m]$ before multiplying it with $s[m]$ point-wise. The output $r[n]$ is a function of time index $n$, which was that shift applied to $h[-m]$.

enter image description here

Next, we turn to the more intuitive method where flipping a signal is not required.

Intuitive Method


There is another method to understand convolution. In fact, it is built on the derivation of convolution equation, i.e., find the output $r[n]$ as \begin{align} r[n] ~ &= ~ \cdots + \nonumber \\ &\qquad\quad s[-2]\cdot h[n+2] ~+ \nonumber \\ &\qquad\qquad \quad\quad\quad s[-1]\cdot h[n+1] ~+ \nonumber \\ &\qquad\qquad \qquad\qquad\quad \quad\quad \quad s[0]\cdot h[n] ~+ \nonumber \\ &\qquad\qquad \quad \qquad\qquad\qquad \quad \quad\quad~ s[1]\cdot h[n-1] ~+ \nonumber \\ &\qquad\qquad \qquad\qquad\qquad\qquad\qquad \quad\quad\quad \quad~ s[2]\cdot h[n-2] ~+ \nonumber \\ &\qquad\qquad \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\quad\quad\quad \cdots \label{eqIntroductionConvolution3} \end{align} Let us solve the same example as in the above Figure, where $s[n] = [2\hspace{1mm}-~\hspace{-2mm}1\hspace{2mm} 1]$ and $h[n] = [-1\hspace{2mm} 1\hspace{2mm} 2]$. This is shown in Table below.

$\hspace{3cm}$enter image description here

Such a method is illustrated in Figure below. From an implementation point of view, there is no difference between both methods.

$\hspace{3cm}$enter image description here

To sum up, convolution tells us how an LTI system behaves in response to a particular input and thanks to intuitive method above, we can say that convolution is also multiplication in time domain (and flipping the signal is not necessary), except the fact that this time domain multiplication involves memory. To further understand at a much deeper level where flipping comes from, and what happens in frequency domain, you can download a sample section from my book here.

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