How to get/approximate the derivative of noisy time series? I have a set of  Langevin equations given by
$${\mathbf{\dot{x}}} = \mathbf{-Q \,x} + \mathbf{\eta} \tag{1}$$ 
where $\eta$ is white Gaussian noise and $Q$ is not a function of $x$.
Using Euler's method for SDE, I generated a time series of $\bf{x}$, which (as expected) is noisy due to $\eta$. Using the time series of $\bf{x}$, generating $\langle(\bf{-Qx} + \bf{\eta})(\bf{-Qx} + \bf{\eta)^T}\rangle$ gave me a finite matrix.
However, for the case where I don't know $\bf{Q}$ and $\eta$, I need to recover $\bf\dot{x}$ from the noisy time series $\bf x$. How do I get $\langle\bf\dot{x}\dot{x}^T\rangle$ from the generated $\bf x$? I hope this can be done and would give me the same value for  $\langle(\bf{-Qx} + \bf{\eta})(\bf{-Qx} + \bf{\eta)^T}\rangle$.
Equation 1 has the solution 
$${x} (t) = e^{-Qt}x(0) +\int_{0}^{t}dt'e^{-(t-t)'Q}\eta(t')$$
but now I want to recover $\bf{\dot{x}}$ from the time series of $\bf{x}$. Any insight on this problem is highly appreciated. Thank you.
 A: For the general question of estimating a derivative from a noisy time series, there exists a fairly large literature; the best tool is probably determined by which modeling assumptions you want to make about the data.
For example, a Savitzky-Golay filter is a technique for smoothing data that also provides an analytic, closed-form estimate for the derivative.
A: to recover $\dot{\mathbf{x}}$ from $\mathbf{x}$, maybe one way is to  first get an estimation of $Q$ and $\eta$ (treat $Q$ as a random variable, independent of all $\eta$). Assume an estimation of $Q$ is $\hat{Q}$, and by simple algebra, $<\dot{\mathbf{x}} \dot{\mathbf{x}}^T> = <(\hat{Q}\mathbf{x})(\hat{Q}\mathbf{x})^T> + var[\eta]$.
Probably the simplest method is to apply linear regression -- your output data is $\dot{\mathbf{x}} \approx \frac{\mathbf{x}(t+\Delta t) - \mathbf{x}(t)}{\Delta t}$, input data is $\mathbf{x}$, and your error is $\eta$.
A: Note that
$\newcommand\mean[1]{\left\langle #1 \right\rangle}
\mean{(-Qx+\eta)(-Qx+\eta)^\top}=\mean{Qxx^\top Q^\top}+\mean{\eta \eta^\top}=\mean{Qxx^\top Q^\top}+B\, \delta(0)$, where $B$ is the covariance matrix of the noise (identity matrix if they are normalized uncorrelated white noise). 
The reason that you are having a hard time fitting with the time series is that the second term is proportional to the Dirac $\delta$-function, while the first term is finite. As you take smaller and smaller $\Delta t$s for numerically calculating the derivative the second term grows larger with a factor $(\Delta t)^{-1}$.
As suggested by Xige Yang, you need to apply a linear regression on $\frac{x(t+\Delta t)-x(t)}{\Delta t}$ to find $Q$. Then, $B$ is given by the covariance of the error multiplied by $\Delta t$.
