How to find all the local minima/maxima in a range

This can be done using event location within NDSolve. I start off as below (note f is slightly modified from what you have, mostly to rescale it).

GetRLine3[MMStdata_, IO_: 1][x_: x] := 
  ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][
     x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], 
    Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ 
     Range[-4, 4, 0.08]}];

xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f = First[100*D[GetRLine3[{data}, 3][x], x]];

We'll recapture f using NDSolve, and locate the points where the derivative vanishes in the process.

vals = Reap[
    soln = y[x] /. 
      First[NDSolve[{y'[x] == Evaluate[D[f, x]], 
         y[-9.9] == (f /. x -> -9.9)}, y[x], {x, -9.9, 30}, 
        Method -> {"EventLocator", "Event" -> y'[x], 
          "EventAction" :> Sow[{x, y[x]}]}]]][[2, 1]];

Visual check:

Plot[f, {x, -9.9, 30}, 
 Epilog -> {PointSize[Medium], Red, Point[vals]}]

Sounds like a job for Ted's RootSearch package.

Clear[f];
SeedRandom[1];
GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO,
  Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], 
  Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f[x_] = First[D[GetRLine3[{data}, 3][x], x]];

Note that I changed the way yo defined f slightly.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Needs["Ersek`RootSearch`"]
{{a, b}} = InterpolatingFunctionDomain[Head[f[x]]];
criticalPoints = x /. RootSearch[f'[x] == 0, {x, a, b}];
mins = Select[criticalPoints, f''[#] > 0 &];
maxs = Select[criticalPoints, f''[#] < 0 &];

Well, I get an error so I'm not certain how well it worked, but it does look good.

Plot[f[x], {x, -9, 30}, Epilog ->{PointSize[Medium],
  Blue, Point[{#, f[#]} & /@ mins],
  Darker[Red], Point[{#, f[#]} & /@ maxs]
}]]

Some comments

• I didn't actually notice the importance of separating the maxima from the minima right away and have edited accordingly
• Of course, it's feasible that $f'(x)=0$ without a max or min occurring at $x$, although I think that the random nature of your function makes the probability of this occurring zero.
• Not sure we can really prove we've got all the extremes.

Edit in Response to Artes

If you use FindMaximum from two different starting values you will likely get two slightly different approximations to the same root. Here's an example with exactly one max.

g[x_] = Exp[-(Cos[x] - x)^2];
x1 = x /. Last[FindMaximum[g[x], {x, 0.8}]];
x2 = x /. Last[FindMaximum[g[x], {x, 0.7}]];
x1 - x2

8.35447*10^-9

Since you're fitting your data with a cubic spline, and defining $f(x)$ as the derivative of the interpolation function, then if we solve for $f'(x)=0$ for the critical points, the resulting function is a piecewise linear function. We can pull the interpolation x-grid data and use it to construct linear functions which are easy to solve for the roots without resorting to the FindMinimum, etc...

Clear[f];
SeedRandom[1];
GetRLine3[MMStdata_, IO_: 1][x_: x] := 
ListInterpolation[#, InterpolationOrder -> IO,Method -> "Spline"] & /@ 
(({{#[[1]]}, #[[2]]}) & /@ # & /@MMStdata);
data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] +
(Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@Range[-4, 4, 0.08]}];
xLimits = {Min@#1, Max@#1} & @@ Transpose[data];
f = First[GetRLine3[{data}, 3][x]]

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
grid = First@InterpolatingFunctionCoordinates[f];

Show[{Plot[f''[x], {x, -9, 30}, PlotRange -> All],ListPlot[Thread[{grid,  
f''[grid]}], PlotStyle -> Red], 
Graphics[Map[{Dashed, Line[#]} &,Partition[Thread
[{grid, f''[grid]}], 2, 1]]]}]

The dashed curve is the reconstruction from the interpolation grid points shown in red. Next, $\it{eqn}$ defines the linear parameterizations where we have an intersection with the $x-$axis if $0<t<1$ when $y=0$. Since the interpolation function is a polynomial, the minima and maxima alternate and can be pulled from the zeros list in order.

NB: The dashed curve is not exactly $f''(x)$ here, as the red interpolation sample points don't lie exactly on the critical points in the first figure. According to the documentation, taking the derivative of an InterpolatingFunction returns a new InterpolatingFunction, however, it seems the grid values are the same for both. I guess you can always find $f'''(x)=0$ to get those.

pts = Thread[{grid, f''[grid]}];
eqn = Map[#[[1]] (1 - t) + #[[2]] (t) &, Partition[pts, 2, 1]];
zeros = Select[
Flatten[Table[{x, t} /. Solve[eqn[[j]] == {x, 0}, {x, t}], {j, 
   Length[eqn]}], 1], 0 <= #[[2]] <= 1 &][[All, 1]];
zerosM = zeros[[1 ;; -1 ;; 2]];
zerosN = zeros[[2 ;; -1 ;; 2]];

Plot[f'[x], {x, -9, 30}, 
Epilog -> {PointSize[Medium], Red, Point[{#, f'[#]} & /@ zerosM], 
Blue, Point[{#, f'[#]} & /@ zerosN]}, PlotRange -> All]