Update
Since Version 12, Mathematica now incorporates a range of (underrated IMHO) regularisation methods to Fit
and FindFit
.
As explained in this question, you can do a non-parametric fit to your data using B-Splines, and differentiate this fit:
pts = Table[{x, Sin[2 Pi x] + RandomReal[{-.15, .15}]}, {x, 0,
1, .0125}];
kfun[n_, d_] :=
Join[ConstantArray[0, d], Range[0, 1, 1/(n - d)],
ConstantArray[1, d]];
uparam[pts_] := N[Range[0, 1, 1/(Length[pts] - 1)]];
mbasis[pts_, n_, d_] :=
With[{param = uparam[pts]},
Table[BSplineBasis[{d, kfun[n, d]}, j - 1, param[[i]]], {i,
Length[param]}, {j, n}]];
Clear[ctrlpts];
ctrlpts[lambda_: 0] :=
With[{mat = mbasis[pts, 25, 3],
reg = SparseArray[{{i_, i_} ->
2., {i_, j_} /; Abs[i - j] == 1 -> -1.}, {25, 25}, 0.]},
LinearSolve[Transpose[mat].mat + 10^(lambda) Transpose[reg].reg,
Transpose[mat].(Last /@ pts)]];
Show[ListPlot[pts, AxesLabel -> {x, y}],
ListLinePlot[{First /@ pts, mbasis[pts, 25, 3].ctrlpts[0.25]} //
Transpose, PlotStyle -> Red]]

Note that the fit does not go through all the points as you requested.
Here we consider an explicit penalty function. The idea here is that
we find the best (spline) weights subject to a prior corresponding to a roughness penalty (which allows us to tune how smooth the spine function should be, which involves adding a tunable cost to unsmooth spline).
The fit is now
controlled by the relative weight of the penalty (given as an argument to
ctrlpts). We can differentiate it:
df[x_] = BSplineFunction[ctrlpts[1], SplineDegree -> 3]'[x];
Plot[{df[x], 2 Pi Cos[2 Pi x]}, {x, 0.05, 0.95}]

There are known methods (such as cross validation) to estimate automatically what the proper amount of smoothing should be, depending on
what it is you want to estimate (the function, its derivative, its second derivative etc.).
Note that the behaviour of your basis function at the edge of the requested interval needs to be addressed depending on what a proper boundary should be.
For instance, the Fourier filtering method presented by others is formally
equivalent to this B-Sline fit, while assuming periodic boundary condition
and a particular choice of Wiener filter.