Suppose we are given moments of a measure on the interval [0,1]. Is there some practical way to recover the measure itself? I am particularly interested in the case where the measure density is given by some smooth positive function, and I would like to determine the behavior of this function near x=0.

There are several approaches to practical resolution of the problem. The moment data can typically be interpreted as linear constraints on the density function that is recovered by minimizing a functional. Depending on the application, the functional might be energy, Shannon entropy, or some other applicationspecific quantity. Take a look at my answer to an earlier question. Despite the fact that the earlier question dealt with multivariate measures, the references in the answer should give you a good head start. 


If you only have a finite number of moments Maximum entropy methods are often a good way to proceed, as was mentioned in Budisic's answer. If you have all the moments, thats equivalent to having the characteristic function $$\hat{f}(\xi) = 1+i\xi \mu_1\frac{1}{2}\xi^2\mu_21/(3!)i\xi^3\mu_3+1/(4!)\xi^4\mu_4+...$$ So in principle you could do the inverse Fourier transform and recover the distribution. If doing the inverse transform is too challenging, and since you said you mostly care about the local behavior around 0, you might also try saying: $$f(0) = \int_{\infty}^{\infty} \hat{f}(\xi)\,d\xi$$ $$f'(0) = \int_{\infty}^{\infty} \xi\hat{f}(\xi)\,d\xi$$ Etc. 


An addendum to Kai's answer: the moments $m_n$ of the measure $\mu$ are the coefficients of the Laurent expansion at infinity of the Cauchy transform $$\hat \mu(z)=\int \frac{d\mu(t)}{zt}$$ of $\mu$ (just expand $(zt)^{1}$ there and substitute it in the integral):$$\hat \mu(z)=\sum_{n=0}^\infty \frac{m_n}{z^{n+1}}. $$ You can recover the original measure (at least its absolutely continuous part) as the difference of the boundary values of $\hat \mu$ on $(0,1)$ (Sokhotski's theorem). From the practical point of view, you could approximate $\hat \mu$ from a finite number of moments using the (diagonal, i.e. both numerator and denominator of degree $\leq n$) Padé approximants at infinity. This will be close to optimal, and the approximation converges uniformly in the whole complex plane except the interval $[0,1]$ with a geometric rate (although the closer to the interval you are, the slower the rate of convergence is). Finally, the denominators $Q_n$ of these Padé approximants will be the orthogonal polynomials with respect to $\mu$. 


An answer on your question is given in the following book: [ST] J. A. Shohat and J. D. Tamarkin "The problem of moments". See [ST], page 90 and what follows. In particular, the case of an absolutely continuous measure is considered on page 95. 

