$\newcommand\R{\mathbb R}\newcommand{\al}{\alpha}\newcommand{\de}{\delta}\newcommand{\J}{\mathcal J}
\newcommand{\be}{\beta}\newcommand{\ep}{\varepsilon}$Yes, such a construction exists.
Indeed, take any real $\ep>0$. For a real $C\ge0$ and $\al\in(0,1]$, let us say that a function $g$ is $(C,\al)$-Hölder on a set $S\subseteq[0,1]$ if $|g(y)-g(x)|\le C|y-x|^\al$ for all $x,y$ in $S$.
Without loss of generality (wlog), the function $f$ is $(1,\al)$-Hölder on $[0,1]$.
Let $Z:=\{z\in I:=[0,1]\colon f(z)=0\}$ and $N:=I\setminus Z$. Wlog, $\{0,1\}\in Z$ (otherwise, extend $f$ appropriately to an interval $[A,B]\supset I$ so that $f(A)=0=f(B)$ and then shrink the interval $[A,B]$ to $I$). So, $N=\bigcup_{J\in\J}J$ for some (countable) set $\J$ of pairwise disjoint nonempty open subintervals of $I$.
It is enough to construct a smooth function $f_\ep$ such that for each $J\in\J$
\begin{equation*}
\|f-f_\ep\|_{C^\be(J)}\le100\ep \tag{10}\label{10}
\end{equation*}
and
\begin{equation*}
|f_\ep|\le|f|\text{ on } J. \tag{20}\label{20}
\end{equation*}
To begin such a construction, take any real $\de>0$. Let $\J_\de$ denote the set of all intervals $J\in\J$ of length $>\de$. Of course, the set $\J_\de$ is finite. Let
\begin{equation*}
f_\ep(x):=0\text{ for }x\in I\setminus\bigcup_{J\in\J_\de}J.
\end{equation*}
If $\de$ is small (which will be henceforth assumed), then on $I\setminus\bigcup_{J\in\J_\de}J$ the function $f-f_\ep=f$ is small and $(1,\al)$-Hölder. So, $f-f_\ep=f$ is $(\de_1,\be)$-Hölder on $I\setminus\bigcup_{J\in\J_\de}J$ for a small $\de_1$. This follows because $|y-x|^\al<<|y-x|^\be$ if $\be\in(0,\al)$ (as given) and $|y-x|<<1$, whereas $|y-x|^\al\asymp|y-x|^\be\asymp1$ if $|y-x|\asymp1$. (We write $E<<F$ if $E=o(F)$, $E\ll F$ if $E=O(F)$, and $E\asymp F$ if $E\ll F\ll E$.)
We shall build the function $f_\ep$ separately on each interval $J=(a,b)\in\J_\de$, in such a manner that
\begin{equation*}
\text{$f_\ep=0$ near the endpoints of $J$. } \tag{30}\label{30}
\end{equation*}
Then the smoothness of $f_\ep$ on each interval $J\in\J$ will be enough for the smoothness of $f_\ep$ on the entire interval $I$.
Take indeed any interval $J\in\J_\de$.
Note that (i) $f(a)=0=f(b)$ and (ii) either $f>0$ on $(a,b)$ or $f<0$ on $(a,b)$. Wlog, $f>0$ (on $(a,b)$). Moreover, $f$ is continuous. So, $M:=f(x_*)\ge f(x)$ for some $x_*\in J$ and all $x\in J$.
For any $c\in(0,M]$, the points
\begin{equation*}
x_+(c):=\min\{x\in[x_*,b)\colon f(x)\le c\},\quad x_-(c):=\max\{x\in(a,x_*]\colon f(x)\le c\}
\end{equation*}
are well defined. Moreover,
\begin{equation*}
f\ge c\text{ on }J_c:=[x_-(c),x_+(c)],\quad f(x_\pm(c))=c. \tag{40}\label{40}
\end{equation*}
Let
\begin{equation*}
g:=(f-c)1_{J_c}. \tag{45}\label{45}
\end{equation*}
Then $g$ is $(1,\al)$-Hölder on $\R$ and hence so is
\begin{equation*}
g_\eta:=g*K_\eta, \tag{47}\label{47}
\end{equation*}
where, for each $\eta\in(0,\de)$, we let $K_\eta$ be any smooth nonnegative function supported on the interval $[-\eta,\eta]$ such that $\int K_\eta=1$.
Moreover, for $x\in J_c=[x_-(c),x_+(c)]$,
\begin{align}
g_\eta(x)&=\int_{-\eta}^\eta du\,K_\eta(u)(f(x-u)-c)1(x-u\in J_c) \notag \\
&\le\eta^\al+\int_{-\eta}^\eta du\,K_\eta(u)(f(x)-c)1(x-u\in J_c) \notag \\
&\le\eta^\al+(f(x)-c)\le f(x) \notag
\end{align}
assuming that $\eta^\al\le c$, which will be indeed assumed henceforth.
Further, for $x\in[x_+(c),x_+(c)+\eta]$,
\begin{align}
g_\eta(x)&=\int_{-\eta}^\eta du\,K_\eta(u)(f(x-u)-c)1(x-u\in J_c) \notag \\
&=\int_{-\eta}^\eta du\,K_\eta(u)(f(x-u)-f(x_+(c))1(x-u\in J_c) \notag \\
&\le\int_{-\eta}^\eta du\,K_\eta(u)(x_+(c)-(x-u))^\al 1(x-u\in [x_-(c),x_+(c)]) \notag \\
&\le\int_{-\eta}^\eta du\,K_\eta(u)\eta^\al 1(x-u\in [x_-(c),x_+(c)])
\le\eta^\al. \notag
\end{align}
On the other hand, again for $x\in[x_+(c),x_+(c)+\eta]$, we have $f(x)\ge f(x_+(c))-(x-x_+(c))^\al\ge c-\eta^\al\ge\eta^\al\ge g_\eta(x)$ assuming that
\begin{equation*}
2\eta^\al\le c, \tag{50}\label{50}
\end{equation*}
which will be indeed assumed henceforth.
Similarly, $g_\eta(x)\le\eta^\al$ for $x\in[x_-(c)-\eta,x_-(c)]$ given \eqref{50}.
Also, $g_\eta(x)=0\le f(x)$ for $x\in J\setminus J_c$.
We conclude that
\begin{equation}
g_\eta\le f\text{ on }J.
\end{equation}
So, letting
\begin{equation*}
f_\ep:=g_\eta\text{ on }J=(a,b), \tag{60}\label{60}
\end{equation*}
we see that $f_\ep$ is smooth on $J=(a,b)$, and conditions \eqref{20} and \eqref{30} hold.
Moreover, $g_\eta$ is $(1,\al)$-Hölder on $\R$ and hence $f_\ep$ is $(1,\al)$-Hölder on $J$.
So, it remains to check that $f_\ep$ is uniformly close to $f$ on $J$, also uniformly in $J\in\J_\de$.
By \eqref{45}, on $J_c$ the function $g$ is uniformly close to $f$ if $c$ is small enough, which will be henceforth assumed.
Also, on $J\setminus J_c$ we have $0\le f\le\max((b-x_+(c))^\al,(x_-(c)-a)^\al)$, which is small (since $c$ was assumed to be small), and hence $f$ is close to $g$ (which is $0$ on $J\setminus J_c$). So, $g$ is uniformly close to $f$ on $J$, uniformly in $J\in\J_\de$.
Finally, because $\eta$ is small, and in view of \eqref{60} and \eqref{47}, we conclude that indeed $f_\ep$ is uniformly close to $f$ on $J$, uniformly in $J\in\J_\de$. $\quad\Box$
Best Answer
Yes, this works, and the only ingredient we need is the estimate $\int_r^{\infty} e^{-t^2}\, dt\lesssim e^{-r^2}$.
We then have (for example) \begin{align*} \int_r^{\infty} |f_r(x)|\, dx &=\frac{1}{\sqrt{\pi}}\int_{-r}^r dy\int_r^{\infty} dx\, e^{-(x-y)^2} =\frac{1}{\sqrt{\pi}}\int_{-r}^r dy\int_{r-y}^{\infty}dt\, e^{-t^2} \\ & \lesssim \int_{-r}^r e^{-(r-y)^2}\, dy \lesssim 1 . \end{align*} Similarly, $$ \int_{-r}^r |1-f_r(x)|\, dx = \frac{1}{\sqrt{\pi}}\int_{-r}^r dx \left( \int_{-\infty}^{\infty}dy\, e^{-(x-y)^2} - \int_{-r}^r dy\, e^{-(x-y)^2} \right) , $$ and now we're back in the same situation as above.