There exists a variety of options for this task but let me stress first that this is an extremely complicated and difficult issue that is still subject of current research because analytical continuation is an ill posed problem!
1) The 'analytical' analytical continuation can be performed when the function $f(\mathrm i\omega)$ under consideration is a rational function of $\mathrm i\omega$. So $$f(\mathrm i\omega)=\frac{1}{\mathrm i\omega}$$ can be continued to the complex plane $\mathrm i\omega \rightarrow z\in\mathbb C$ while $$f(\mathrm i\omega)=\frac{e^{\mathrm i\omega\beta}}{\mathrm i\omega}$$ is not a rational function of $\mathrm i\omega$ and making the replacement here is a mistake. Instead one needs to evaluate the exponential first and find $e^{\mathrm i\omega\beta}=\pm 1$ depending on the statistics.
2) Directly inferred from this replacement rule comes the expansion of a function in a finite Laurent series
$$f(z)=\sum^{m_2}_{n=m_1} a_n z^n,\;\; m_1,m_2\in\mathbb Z$$
where the coefficients may be calculated from the numerical values known at $m_2-m_1$ Matsubara energies.
3) One of the oldest methods to do a numerical analytical continuation is the Pade approximation. The function in question is expanded in a continued fraction
$$f(z)=b_0+\frac{a_1z}{1-\frac{a_2z}{1-\frac{a_3z}{1-...}}}.$$
The coefficients can be computed from a Pade table, see http://en.wikipedia.org/wiki/Pad%C3%A9_table
Method 1) is exact and other than for almost trivial calculations of little practical value. 2) and 3) suffer from cutoff effects due to the limited amount of available Matsubara points at which the function value might also have a numerical error as is the case for data from Quantum Monte Carlo calculations. But in fact analytical continuation is very volatile towards cutoff and noise effects. This is where physical considerations have to be accounted for.
To tackle the cutoff one can approximate the tail (large $\mathrm i\omega$ or respectively $z$ expansion) of the function with an analytical form that can often times be computed exactly from the many body problem or general physical necessities, e.g. the 1-Particle Green's function of a fermionic system always has the form $\frac{1}{z}+\frac{a_1}{z^2}+...$. The tail can be used to compute an arbitrary number of expansion coefficients but keep in mind that the interesting low energy spectrum of your system is strongly influenced by small Matsubara energies and less so by the tail so from computing a large number of coefficients from the tail one gains little to nothing.
The treatment of statistical noise is even more delicate than the cutoff and the reason why a lot of people try to avoid calculating on the Matsubara axis altogether.
4) A prominent method for noisy data is the maximum entropy method about which you can read more here http://arxiv.org/pdf/1001.4351v1.pdf where you will also find references to alternative techniques.
