Convergence rate for $L^2$ convergence Let $f \in L^2(\mathbb R)$ then it is well-known that 
$$ \widetilde{f}(x):=\sum_{n \in \mathbb Z} \frac{1}{\varepsilon}\int_{[n\varepsilon,(n+1)\varepsilon]} f(s) \ ds 1_{[n\varepsilon,(n+1)\varepsilon)}(x)$$
converges in the $L^2$ sense to $f.$
But even more is true, as we can write 
$$(\widetilde{f}-f)(x):=\sum_{n \in \mathbb Z} \frac{1}{\varepsilon}\int_{[n\varepsilon,(n+1)\varepsilon]} f(s)-f(x) \ ds 1_{[n\varepsilon,(n+1)\varepsilon)}(x).$$
Then, it follows that if $f$ is absolutely continuous that
$$(\widetilde{f}-f)(x):=\sum_{n \in \mathbb Z} \frac{1}{\varepsilon}\int_{[n\varepsilon,(n+1)\varepsilon]} \int_x^s f'(\tau) d \tau \ ds 1_{[n\varepsilon,(n+1)\varepsilon)}(x).$$
Hence, it is tempting to ask whether there are conditions on $f$ being in some Sobolev space such that
$$\left\lVert \widetilde{f}-f \right\rVert_{L^2}=\mathcal O(\varepsilon).$$
Moreover, perhaps this generalizes to higher dimensions, i.e. $\mathbb R^d$ rather than $\mathbb R.$
It seems that what would be sufficient is that $$\sum_{n \in \mathbb Z} \frac{1}{\varepsilon} \sup_{s \in [n \varepsilon, (n+1)\varepsilon)}\left\lvert f'(s)\right\rvert 1_{[n\varepsilon,(n+1)\varepsilon)}(x)$$ is square-integrable. However, I fail to see whether there is some natural Sobolev space in which functions have this property.
 A: $\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}} 
\newcommand{\PP}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}} 
\newcommand{\tf}{\widetilde{f}}$ 
As you noted,
\begin{equation*}
 \tf-f=\frac{1}{\ep}\sum_n a_n 1_{[n\ep,(n+1)\ep)},
\end{equation*}
where 
\begin{equation*}
 a_n(x):=\int_{n\ep}^{(n+1)\ep} \int_x^s f'(\tau) d \tau \ ds. 
\end{equation*}
For $x\in[n\ep,(n+1)\ep)$, 
\begin{equation*}
 |a_n(x)|\le\int_{n\ep}^{(n+1)\ep}\int_{n\ep}^{(n+1)\ep} |f'(\tau)| d \tau\ ds
 =\int_{n\ep}^{(n+1)\ep} |f'(\tau)| d \tau\ \ep\le b_n\ep\sqrt\ep,
\end{equation*}
where 
\begin{equation*}
 b_n:=\sqrt{\int_{n\ep}^{(n+1)\ep} |f'(\tau)|^2 d \tau}. 
\end{equation*}
So, 
\begin{equation*}
|\tf-f|\le\sum_n b_n \sqrt\ep \, 1_{[n\ep,(n+1)\ep)}
\end{equation*}
and 
\begin{equation*}
 \|\tf-f\|_2\le\sqrt{\sum_n b_n^2 \ep^2}=\ep\|f'\|_2\le\ep\|f\|_{W^{1,2}}, 
\end{equation*}
as desired. 

Consider now the general case when $f$ is a function on $\R^d$ with $d\ge2$ and 
\begin{equation*}
 \tf=\frac{1}{\ep^d}\sum_{n\in\Z^d} 1_{\de_n}\int_{\de_n} ds\,f(s),
\end{equation*}
where $\de_n:=\prod_1^d\,(n_i\ep,(n_i+1)\ep)$ for $n=(n_1,\dots,n_d)$. 
Take any real 
\begin{equation*}
 p>d
\end{equation*}
and let $q:=\frac p{p-1}$, so that $\frac1p+\frac1q=1$. 
We have
\begin{equation*}
 \tf-f=\frac{1}{\ep^d}\sum_n a_n 1_{\de_n}, \tag{!}
\end{equation*}
where 
\begin{equation*}
 a_n(x):=\int_{\de_n} ds \int_0^1 dt\, f'(x+t(s-x))\cdot(s-x) 
\end{equation*}
$f':=\nabla f$ and $\cdot$ is the dot product. 
By the change of variables from $s$ to $\tau=x+t(s-x)$, for $x\in\de_n$ we have 
\begin{equation*}
 a_n(x)=\int_0^1 dt\int_{\de_n} \frac{d\tau}{t^d} \, f'(\tau)\cdot\frac{\tau-x}t\,1_{(t_*,\infty)}(t),
\end{equation*}
where 
\begin{align*}
t_*:=t_*(n,x,\tau)&:=\max_1^d\Big(\frac{\tau_i-x_i}{x_i-n_i\ep}\vee\frac{x_i-\tau_i}{(n_i+1)\ep-x_i}\Big) \\ 
&\ge\max_1^d\frac{|\tau_i-x_i|}{\ep}\gg\frac{|\tau-x|}{\ep},
\end{align*}
with the constants in $\gg$ and $\ll$ depending only on $d$ and $p$, and
$|\ \, |$ also denoting the Euclidean norm on $\R^d$. Note that 
\begin{equation*}
 \int_{t_*}^\infty\frac{dt}{t^{d+1}}\ll\frac1{t_*^d}\ll\frac{\ep^d}{|\tau-x|^d}. 
\end{equation*}
It follows that for $x\in\de_n$ 
\begin{equation*}
 |a_n(x)|\ll\ep^d\,\int_{\de_n} \frac{d\tau}{|\tau-x|^{d-1}} \, |f'(\tau)|
 \le\ep^d\,b_n c^{1/q}, \tag{!!}
\end{equation*}
by H\"older's inequality, where 
\begin{equation*}
 b_n:=\Big(\int_{\de_n}d\tau\,|f'(\tau)|^p\Big)^{1/p}, 
\end{equation*}
\begin{equation*}
 c:=\int_{|\tau-x|\le\ep\sqrt d}\frac{d\tau}{|\tau-x|^{(d-1)q}}
 \ll\int_0^{\ep\sqrt d}\frac{dr}{r^{(d-1)q}}\,r^{d-1}
 \ll\ep^{(p-d)/(p-1)}. 
\end{equation*}
Thus, by (!) and (!!), 
\begin{equation*}
 \|\tf-f\|_p\ll\Big(\sum_n b_n^p \ep^{p-d}\ep^d\Big)^{1/p}=\ep\|f'\|_p\le\ep\|f\|_{W^{1,p}}, 
\end{equation*}
as desired. 
