Why correlation functions? While this concept is widely used in physics, it is really puzzling (at least for beginners) that you just have to multiply two functions (or the function by itself) at different values of the parameter and then average over the domain of the function keeping the difference between those parameters:
$$C(x)=\langle f(x'+x)g(x')\rangle$$
Is there any relatively simple illustrative examples that gives one the intuition about correlation functions in physics? 
 A: The correlation function you wrote is a completely general correlation of two quantities,
$$\langle f(X) g(Y)\rangle$$
You just use the symbol $x'$ for $Y$ and the symbol $x+x'$ for $X$.
If the environment - the vacuum or the material - is translationally invariant, it means that its properties don't depend on overall translations. So if you change $X$ and $Y$ by the same amount, e.g. by $z$, the correlation function will not change.
Consequently, you may shift by $z=-Y=-x'$ which means that the new $Y$ will be zero. So
$$\langle f(X) g(Y)\rangle = \langle f(X-Y)g(0)\rangle = \langle f(x)g(0) \rangle$$
As you can see, for translationally symmetric systems, the correlation function only depends on the difference of the coordinates i.e. separation of the arguments of $f$ and $g$, which is equal to $x$ in your case.
So this should have explained the dependence on $x$ and $x'$.
Now, what is a correlator? Classically, it is some average over the probabilistic distribution
$$\langle S \rangle = \int D\phi\,\rho(\phi) S(\phi)$$
This holds for $S$ being the product of several quantities, too. The integral goes over all possible configurations of the physical system and $\rho(\phi)$ is the probability density of the particular configuration $\phi$.
In quantum mechanics, the correlation function is the expectation value in the actual state of the system - usually the ground state and/or a thermal state. For a ground state which is pure, we have
$$\langle \hat{S} \rangle = \langle 0 | \hat{S} | 0 \rangle$$
where the 0-ket-vector is the ground state, while for a thermal state expressed by a density matrix $\rho$, the correlation function is defined as
$$\langle \hat{S} \rangle = \mbox{Tr}\, (\hat{S}\hat{\rho})$$
Well, correlation functions are functions that know about the correlation of the physical quantities $f$ and $g$ at two points. If the correlation is zero, it looks like the two quantities are independent of each other. If the correlation is positive, it looks like the two quantities are likely to have the same sign; the more positive it is, the more they're correlated. They're correlated with the opposite signs if the correlation function is negative.
In quantum field theory, correlation functions of two operators - just like you wrote - is known as the propagator and it is the mathematical expression that replaces all internal lines of Feynman diagrams. It tells you what is the probability amplitude that the corresponding particle propagates from the point $x+x'$ to the point $x'$. It is usually nonzero inside the light cone only and depends on the difference of the coordinates only.
An exception to this is the Feynman Propagator in QED. It is nonzero outside the light cone as well, but invokes anti-particles, which cancel this nonzero contribution outside the light cone, and preserve causality.
Correlation functions involving an arbitrary positive number of operators are known as the Green's functions or $n$-point functions if a product of $n$ quantities is in between the brackets. In some sense, the $n$-point functions know everything about the calculable dynamical quantities describing the physical system. The fact that everything can be expanded into correlation functions is a generalization of the Taylor expansions to the case of infinitely many variables.
In particular, the scattering amplitude for $n$ external particles (the total number, including incoming and outgoing ones) may be calculated from the $n$-point functions. The Feynman diagrams mentioned previously are a method to do this calculation systematically: a complicated correlator may be rewritten into a function of the 2-point functions, the propagators, contracted with the interaction vertices.
There are many words to physically describe a correlation function in various contexts - such as the response functions etc. The idea is that you insert an impurity or a signal into $x'$, that's your $g(x')$, and you study how much the field $f(x+x')$ at point $x+x'$  is affected by the impurity $g(x')$.
A: There is another reason, although it is not intuitive.  A stochastic process is almost completely characterised by its auto-correlation function.  More precisely, if the process is stationary (of course, all these methods only work after a process has been de-trended and all cycles analysed and filtered out first) and Gaussian, and centred,  then it is completely characterised by the auto'correlation function.  This is analogous to the elementary fact that a normal random variable is completely characterised by its mean and standard deviation.
But wait.  There's more.  Even if the process is not Gaussian, it is characterised if you know not just the usual autocorrelation, which is also called the two-point correlation functions, but if you know all the higher auto-correlations as well.  (i.e, three-point, four-point, etc.).  This is analogous to the (difficult) "problem of moments". solved by my (academic PhD. genealogy) ancestor Marcel Riesz and that alcoholic genius of Sweden, Carlman, that says if you know all the moments of a random variable, it is determined up to equivalence.
And, in practice, it is the correlations that are the most accessible to measurement.  Most experiments with particles, including the famous Bell inequality experiments of Aspect, are measurements of correlations.  There is probably some profound philosophical significance to this....
A: a very intuitive example for correlation functions can be seen in laser speckle metrology. 
If you shine light on a surface which is rough compared to the wavelength, the resulting reflected signal will be somehow random. This can also be stated as that you cannot say from one point of a signal how a neighbouring one looks like - they are uncorrelated. Such a field is often referred to as speckle pattern.
This fact can be used. Suppose you take an image $A(x,y)$ of such a random scattered field, a movement of the image $$(x,y)\rightarrow (x+\delta_x, y+\delta_y)  = (x',y')$$ thus $$B(x,y) \approx A(x',y')$$
will be clearly visible and since all information are statistical, one finds that 
$$C(\Delta_x,\Delta_y) = \int B(x,y) A(x + \Delta_x, y + \Delta_y) dx dy $$
will only have a "big" contribution at $(\Delta_x,\Delta_y) \equiv (\delta_x, \delta_y)$ of some peaked form. The width of the peak will be given by some physical properties of the illumination, roughness of the surface etc. - it directly corresponds to the local variation of the field. 
If we had now in the field some periodic variation we could see that $C$ will have several peaks corresponding to the image's (or field's) self-similarity.  
So, analyzing the correlation of a quantity will give you information on how fast it changes and if it is somehow self-similar.
I hope you don't mind that I have chosen an application coming from a more practical viewpoint. 
Sincerely  
Robert
PS.: More can be found in all the very rich works done by Goodman.
A: Excellent question, Kostya. Lubos already gave a detailed answer using general arguments in the language of QFT.
In astrophysics and cosmology, however, there is another, and very simple, reason why we use the correlation functions all the time. It turns out that the mean value of the function $f(\vec{x})$, denoted $\langle f(\vec{x})\rangle$, can often not be predicted by the theoretical model (e.g. hot Big Bang model with inflationary stage early on, cold dark matter at late times, etc... or whatever other model you wish to consider) -  while its correlation $\langle f(\vec{x})f(\vec{y})\rangle$ can be predicted. Here $f$ can refer to any cosmological observable quantity, and $\vec{x}$ and $\vec{y}$ refer to spatial coordinates.
The most common example would be to consider the excess density of dark matter, $f(\vec{x})\equiv \delta\rho(\vec{x})/\rho$, where $\rho$ is the mean density (units of which are kilograms per meter cubed for example) and $\delta\rho(\vec{x})$ is the excess over- or under-density at location $\vec{x}$, and over some region which I will not specify for simplicity of the argument. By definition, the mean of $f$ is zero, so we explicitly indicate that we are not interested in the mean (alternatively, we cannot easily get the mean density of the universe from first principles). But the correlation function, $\langle \delta\rho(\vec{x})\delta\rho(\vec{y})/\rho^2\rangle $ can be related to fundamental parameters of the universe, in particular details of the inflationary epoch, dark matter density, etc. Details of this are involved, and are taught in a graduate course in cosmology. Suffice it to say that theory predicts not the mean of the function (1-point correlation function), but rather its (co)variances (2-point correlation function). 
Intuitively, the two-point correlation function of $\delta\rho/\rho$ is related to the "probability that, given an overdense region of dark matter at location $\vec{x}$, there is an overdense region at location $\vec{y}$", and this probability is determined by the good old law of gravity - and can be predicted from first principles. 
Theory also in principle predicts the 3-point (e.g. $\langle f(\vec{x})f(\vec{y})f(\vec{z})\rangle$, and higher-point correlation functions, but those are both harder to calculate theoretically and measure observationally.  Nevertheless, there is a thriving subfield in particle physics and cosmology of predicting theoretically, and measuring observationally, these so-called higher-order correlation functions. 
One final ingredient in all this is the role of measuring the correlation function. The angular averaging sign, $\langle\cdot\rangle$ implies that we should be averaging over different realizations of the system - that is, the universe - in the same underlying cosmological model. This is clearly impossible, since we have only one universe to measure! Instead, we assume statistical homogeneity (which is the same as translational invariance from Lubos' post). Then, instead of averaging over the different universes, cosmologists average $f(\vec{x})f(\vec{y})$ over different locations ($\vec{x}, \vec{y}$) in our universe which have a fixed distance between the two points $|\vec{x}-\vec{y}|$. This way, using the statistical homogeneity assumption, we can get good measurements of the correlation function of any quantity we desire. 
