The estimation of derivative is straightforward:
$$x'(n)~=\frac{x(n+1)-x(n-1)}{2}$$
$$x''(n)~={x(n+1)-2*x({n})+x({n-1})}$$
or if you have a signal sampled at $t_i=i\Delta t$,
it is
$$x'(t_{i})~=\frac{x(t_{i+1})-x(t_{i-1})}{2*\Delta t}$$
$$x''(t_{i})~=\frac{x(t_{i+1})-2*x(t_{i})+x(t_{i-1})}{(\Delta t)^2}$$
What you are interested in may be how to smooth the estimation. And yes you can use some recursive filters such as $y(n) = a \cdot x(n) + (1-a) \cdot y(n-1)$, or just implement your estimation with some simple windows (Hann window for example).
To achieve the high SNR without distorting the signal very much, Savitzky–Golay filter also smooths your data by fitting successive sub-sets of adjacent data points with a low-degree polynomial with linear least squares.
EDIT
Matlab code for an N-th derivative of signal row vector x
dx = x; %'Zeroth' derivative
for n = 1:N % Apply iteratively
dif = diff(dx,1); % First derivative
first = [dif(1) dif];
last = [dif dif(end)];
dx = (first+last)/2;
end