Ok, I found the real solution for my problem. I am posting it even if the correct answer to my (misplaced) question is the one selected.

Basically, my question explains how to estimate the mean knowing the covariance, and the answer how to estimate the covariance knowing the mean. But my actual problem was estimating with both parameters unknown.

I found the answer on [Wikipedia][1] with the derivation explained [here][2]. The multivariate normal's conjugated prior is the Normal-inverse-Wishart, that is basically a distribution over multivariate Normals.

The prior parameters that need to be specified are $\mathbf{\mu}_0$ to define the mean, $\mathbf{\Psi}$ to define the covariance, and two scalar values $\kappa_0$ and $\nu_0$ that I would say define how confident we are on the estimation of the first two parameters respectively.

The updated distribution after observing $n$ samples of a $p$-variate Normal has the form 

$$\mathbf{P}(\boldsymbol\mu, \mathbf{\Sigma} | \mathbf{X}) \sim \mathcal{NIW} \left( \frac{\kappa_0\boldsymbol\mu_0+n\mathbf{\bar{x}}}{\kappa_0+n} ,\, \kappa_0+n,\, \nu_0+n ,\,  \boldsymbol\Psi + \mathbf{C} + \frac{\kappa_0 n}{\kappa_0+n}(\mathbf{\bar{x}}-\boldsymbol\mu_0)(\mathbf{\bar{x}}-\boldsymbol\mu_0)^T \right)$$

where 

$$\mathbf{\bar{x}} = {1 \over n} \sum_{i=0}^{n} \mathbf{x_i} $$

$$\mathbf{C} = \sum_{i=1}^n (\mathbf{x_i} - \mathbf{\bar{x}}) (\mathbf{x_i} - \mathbf{\bar{x}})^T $$

so my desired estimated parameters are

$$ E (\boldsymbol\mu | \mathbf{X}) = {{\kappa_0\boldsymbol\mu_0+n\mathbf{\bar{x}}} \over{\kappa_0+n} } $$
$$ E (\mathbf{\Sigma} | \mathbf{X}) = \frac{ \boldsymbol\Psi + \mathbf{C} + \frac{\kappa_0 n}{\kappa_0+n}(\mathbf{\bar{x}}-\boldsymbol\mu_0)(\mathbf{\bar{x}}-\boldsymbol\mu_0)^T }{ \nu_0 + n - p - 1}$$

  [1]: http://en.wikipedia.org/wiki/Conjugate_prior
  [2]: http://www.cs.ubc.ca/~murphyk/Papers/bayesGauss.pdf