History of convolution Let $f, g\in L^{1}(\mathbb R),$
we may define the convolution of $f$ and $g$ as follows:
$f\ast g(x)= \int_{\mathbb R} f(x-y)g(y) dy, (x\in \mathbb R).$
It is well known that it can be defined on general locally compact group as well.
It the operation convolution (I think)  in analysis (perhaps, in other branch of mathematics as well) is  like one of the most useful operation (perhaps after the four fundamental operations addition, subtraction, multiplication, division) 

MY Question: How old  the  operation convolution  is? In other words, the idea of convolution goes back to whom? Why does it get this much attention in mathematics?

 A: Citing from jeff560.tripod.com:

CONVOLUTION. Expressions that would now be described as “convolutions”
  appear in Laplace’s earliest work on sums of independent random
  variables, “Mémoire sur l’inclinaison moyenne des orbites des comètes,
  sur la figure de la terre, et sur les functions,” Mém. Acad. R. Sci.
  Paris (Savants Étrangers), 7, (1773), 503-540, OC 8, 279-321. A
  succession of French and Russian mathematicians followed Laplace and
  used convolutions without, it seems, evolving a name for them. See A.
  Hald History of Mathematical Statistics from 1750 to 1930 (1998).
In the early 20th century convolutions appeared in the field of
  INTEGRAL EQUATIONS. In his 1906 paper “Grundzüge einer allgemeinen
  Theorie der linearen Integralgleichungen. Vierte Mittelung”
  Nachrichten von d. Königl. Ges. d. Wissensch. zu Göttingen
  (Math.-physik. Kl.) (1906) 157-227 Hilbert used the word Faltung,
  meaning “folding” or “plaiting.” See J. Dieudonné History of
  Functional Analysis (pp. 113-4). The term became standard in the new
  field of FUNCTIONAL ANALYSIS.
In the 1930s an English equivalent was found. In 1932 Aurel Wintner
  used “folding expression” in his “Remarks on the Ergodic Theorem of
  Birkhoff,” Proceedings of the National Academy of Sciences, 18, (3),
  p. 251. Norbert Wiener stuck to Faltung in his 1933 book, The Fourier
  Integral and Certain of its Applications, maintaining that “there is
  no good English word” (p. 45) Actually there was a rarely used and
  rather formal English synonym for “folding” and Wintner used
  “convolution” in his 1934 article, “On Analytic Convolutions of
  Bernoulli Distributions,” American Journal of Mathematics, 56, p. 662.
  This became the usual English term. See the Encyclopaedia of
  Mathematics entry.
[This entry was contributed by John Aldrich. A citation was provided
  by Yaakov Stein.]

About the importance:
The convolution has a nice property, when Fourier transformed:
$$
\mathcal{F}(f*g) = c \, \mathcal{F}(f) \mathcal{F}(g)
$$
for some constant specific to the definition of the Fourier transform used. 
And it is the continous cousin of the Cauchy product.
A: Note: I really like the other answer. This is just a side note to add a little intuition as to why convolution has a useful function in Mathematics.
To see why it gets so much attention, consider the product of two power series
\begin{align}            \sum_{n=0}^{\infty}f(n)z^{n}\sum_{n=0}^{\infty}g(n)z^{n}
  & = \sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}f(k)g(n-k)\right)z^{n} \\
  & = \sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}f(n-k)g(k)\right)z^{n}
\end{align}
The coefficient of $z^{n}$ is the sum of all coefficients with powers that sum to $n$, i.e., $(n-k)+k=n$. The Laplace transform can be viewed as a generalization of a power series
$$
             \mathscr{L}\{f\}(s)=\int_{0}^{\infty}f(t)(e^{-s})^{t}dt.
$$
In fact, the Mellin transform makes this more explicit
$$
                       \int_{0}^{\infty}f(t)z^{t}dt.
$$
The product of two Laplace transforms is another object of the same type after the common powers of $e^{-s}$ are collected by convolution:
\begin{align}
        & \int_{0}^{\infty}f(t)e^{-st}dt\int_{0}^{\infty}g(t)e^{-st}dt \\
   & = \int_{0}^{\infty}\left(\int_{0}^{t}f(u)g(t-u)du\right)e^{-st}dt \\
   & = \int_{0}^{\infty}\left(\int_{0}^{t}f(t-u)g(u)du\right)e^{-st}dt
\end{align}
When dealing with the Fourier transform, all positive and negative powers must be collected:
\begin{align}
    & \int_{-\infty}^{\infty}f(x)e^{-isx}dx\int_{-\infty}^{\infty}g(x)e^{-isx}dx \\ & =\int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty}f(u)g(x-u)du\right)e^{-isx}dx \\
   & = \int_{-\infty}^{\infty}\left(\int_{-\infty}^{\infty}f(x-u)g(u)du\right)e^{-isx}dx
\end{align}
You can see how this might extend to groups.
A: I heartily recommend the part 23, Haar Measure. Convolution,  of the book Bourbaki. Elements of the History of Mathematics, Springer. I have the book but in Spanish. This chapter 23, not numbered as such in the original version, is truly extraordinary and bright. I copy here from Internet  the beginning of this part 23 which unfortunately does not reach the convolution that is treated in the last part of the chapter and begins: 
“The history of the convolution product, is more complex”.

