Let $\delta_n:=\frac{2n}{(2n+1)n(n-1)}$, and $g_n$ defined by
$$g_n(x)=\begin{cases}
n&\mbox{ on }\left(\frac kn-\frac{\delta_n}2,\frac kn+\frac{\delta_n}2\right), 1\leq k\leq n-1\\\
\mbox{ linear }&\mbox{ on }\left(\frac kn-\frac{\delta_n}2-\frac{\delta_n}{2n},\frac kn-\frac{\delta_n}2\right)\\\
\mbox{ linear }&\mbox{ on }\left(\frac kn+\frac{\delta_n}2,\frac kn+\frac{\delta_n}2+\frac{\delta_n}{2n}\right)\\\
0&\mbox{ elsewhere}.
\end{cases}$$
We have that the measure of the support of $g_n$ is $(n-1)(1+1/n)\delta_n$ which converges to $0$, and $\int_{[0,1]}g_ndm=1$, except miscomputation. We can write
$$\left|\int_{[0,1]}g_nfdm-\frac 1n\sum_{k=1}^nf\left(\frac kn\right)\right|\leq 2\delta_n \lVert f\rVert_{\infty}+\operatorname{mod}(f,\delta_n)\frac{n^2}{\delta_n},$$
where $\operatorname{mod}(f,\delta):=\sup\{|f(x)-f(y)|,x,y\in I, |x-y|\leq \delta \}$.
So the three conditions are full-filled by $m$ and $m$ are of course not singular.
Finally, after more extensive searching, I found this post Question about Lebesgue point in Fourier Transform (Big Rudin chapter 9) which asks the same question as I did
and then provides an answer. The answer depends on material from chapter 11 of Rudin Real and Complex Analysis (which is two chapters past the one in which the original exercise appears). So I decided to distill the part of the relevant proof in Rudin (Theorem 11.20)
and the work from the above quoted post into a self-contained answer that, in principle, could have been worked out by anyone who studied up to the point where the exercise appears.
[Equation numbers refer to both this post and my original post.]
Let $f\in L^1$ and let $x$ be a Lebesgue point of $f$. As in Rudin chapter 9,
let $m$ be Lebesgue measure on $R$ divided by $\sqrt{2\pi}$. Define
$$k_x(y)=\chi_{[x-1,x+1]}(y)f(x)\qquad(y\in R).$$
Then $k_x\in L^1$ and by (1)
\begin{equation}
\begin{split}
(k_x*h_\lambda)(x)
&=\int_{-\infty}^\infty k_x(x-y)h_\lambda(y)\,dm(y)
=\int_{-\infty}^\infty\chi_{[x-1,x+1]}(x-y)f(x)h_\lambda(y)\,dm(y)\\
&=\int_{-1}^1 f(x)h_\lambda(y)\,dm(y)
=f(x)\int_{-1}^1\sqrt{\frac{2}{\pi}}\frac{\lambda}{\lambda^2+y^2}\,dm(y)\\
&=-f(x)\frac{1}{\sqrt{2\pi}}\sqrt{\frac{2}{\pi}}
\int_{\cot^{-1}(-\lambda^{-1})}^{\cot^{-1}(\lambda^{-1})}dt
\qquad\text{(having substituted $t=\cot^{-1}(y/\lambda)$)}\\
&=-\frac{f(x)}{\pi}[\cot^{-1}(\lambda^{-1})-\cot^{-1}(-\lambda^{-1})]
=-\frac{f(x)}{\pi}[\cot^{-1}(\lambda^{-1})-(\pi-\cot^{-1}(\lambda^{-1}))]\\
&=f(x)[1-\frac{2}{\pi}\cot^{-1}(\lambda^{-1})]
\end{split}
\end{equation}
so
\begin{equation}
\lim_{\lambda\to 0}\,(k_x*h_\lambda)(x)
=f(x)-\frac{2f(x)}{\pi}\lim_{\lambda\to 0}\cot^{-1}(\lambda^{-1})=f(x).
\end{equation}
Let $f_x=f-k_x\in L^1$. Then
\begin{equation}
\begin{split}
\lim_{\lambda\to 0}\,(f*h_\lambda)(x)
&=\lim_{\lambda\to 0}\,((f_x+k_x)*h_\lambda)(x)
=\lim_{\lambda\to 0}\,(f_x*h_\lambda)(x)
+\lim_{\lambda\to 0}\,(k_x*h_\lambda)(x)\\
&=\lim_{\lambda\to 0}\,(f_x*h_\lambda)(x)+f(x).
\end{split}
\end{equation}
To complete the proof, we will show that $\lim_{\lambda\to 0}\,(f_x*h_\lambda)(x)=0$.
Let $\epsilon>0$ be given. Since $x$ is a Lebesgue point of $f$,
\begin{equation}
\begin{split}
\lim_{\lambda\to 0}\frac{1}{2\lambda}
\int_{x-\lambda}^{x+\lambda}|f_x(y)|\,dm(y)
&=\lim_{\lambda\to 0}\frac{1}{2\lambda}
\int_{x-\lambda}^{x+\lambda}|f(y)-k_x(y)|\,dm(y)\\
&=\lim_{\substack{\lambda\to 0\\\lambda\in(0,1)}}
\frac{1}{2\lambda}\int_{x-\lambda}^{x+\lambda}|f(y)-f(x)|\,dm(y)=0,
\end{split}
\end{equation}
since, in the last integral, restricting $\lambda$ to $(0,1)$ restricts
$y$ to $[x-1,x+1]$ where $k_x(y)=f(x)$. Therefore, there exists $\delta>0$
such that for every $\lambda\in(0,\delta]$
\begin{equation}\tag{6}\label{6}
\frac{1}{2\lambda}\int_{x-\lambda}^{x+\lambda}|f_x(y)|\,dm(y)
<\frac{\epsilon}{\sqrt{2\pi}}.
\end{equation}
Let $\{\lambda_n\}\subseteq(0,1)$ be such that $\lim_{n\to\infty}\,\lambda_n=0$
and consider the sequence $\{s_n\}$ where
$$s_n(y)=h_{\lambda_n}(x-y)\chi_{[x-\delta,x+\delta]^c}(y)|f_x(y)|\qquad(y\in R).$$
Since $\lambda_n\in(0,1)$ we have that
$\lambda_n^2+(x-y)^2(1-\lambda_n)>0$ so that
$\lambda_n^2+(x-y)^2>\lambda_n(x-y)^2$ whence
$s_n(y)\leq|f_x(y)|/\delta^2\in L^1$.
Now $\lim_{n\to\infty}s_n(y)=0$ for all $y\in R$,
so by the Dominated Convergence Theorem,
$$\lim_{n\to\infty}\int_{-\infty}^\infty s_n(y)\,dm(y)=0.$$
Therefore
\begin{equation}\tag{7}\label{7}
\lim_{\lambda\to 0}\int_{[x-\delta,x+\delta]^c}|f_x(y)|h_\lambda(x-y)\,dm(y)=0.
\end{equation}
Our next goal is to prove that
\begin{equation}\tag{8}\label{8}
\sup_{\lambda>0}\int_{x-\delta}^{x+\delta}|f_x(y)|h_\lambda(x-y)\,dm(y)
\leq\sup_{\lambda\in(0,\delta]}\frac{\sqrt{2\pi}}{2\lambda}
\int_{x-\lambda}^{x+\lambda}|f_x(y)|\,dm(y)\equiv M.
\end{equation}
(This will be an adaptation of the second half of Rudin's proof of Theorem 11.20.)
Let $\lambda>0$. Given $n\geq 0$, let $h_{2^n+1}^{(n)}=0$ and for each
$j=1,2,\dots,2^n$ let $I_j^{(n)}=[x-j\delta/2^n,x+j\delta/2^n]$,
$h_j^{(n)}=h_\lambda(j\delta/2^n)$, and define
$$s_n=\sum_{j=1}^{2^n}(h_j^{(n)}-h_{j+1}^{(n)})\chi_{I_j^{(n)}}.$$
Then each $s_j$ is a simple measurable function, and for all $y\in R$
\begin{equation}
0\leq s_1(y)\leq s_2(y)\leq\cdots\leq\chi_{[x-\delta,x+\delta]}(y)h_\lambda(x-y)
\text{ and}\\
\lim_{n\to\infty}\,s_n(y)=\chi_{[x-\delta,x+\delta]}(y)h_\lambda(x-y).
\end{equation}
In words, we are constructing an increasing sequence of step approximations
from below to
$\chi_{[x-\delta,x+\delta]}(y)h_\lambda(x-y)$. We rely on the fact that
$h_\lambda(t)$ is even and strictly decreases from a maximum at $t=0$ as $|t|$
increases.
We then have that for all $y\in R$
\begin{equation}
0\leq|f_x(y)|s_1(y)\leq|f_x(y)|s_2(y)\leq\cdots\text{ and}\\
\lim_{n\to\infty}|f_x(y)|s_n(y)
=|f_x(y)|\chi_{[x-\delta,x+\delta]}(y)h_\lambda(x-y).
\end{equation}
By the Monotone Convergence Theorem applied twice,
\begin{equation}
\begin{split}
\int_{x-\delta}^{x+\delta}|f_x(y)|h_\lambda(x-y)\,dm(y)
&=\lim_{n\to\infty}\int_{-\infty}^\infty|f_x(y)|s_n(y)\,dm(y)\\
&=\lim_{n\to\infty}\sum_{j=1}^{2^n}(h_j^{(n)}-h_{j+1}^{(n)})
\int_{x-j\delta/2^n}^{x+j\delta/2^n}|f_x(y)|\,dm(y)\\
&\leq\lim_{n\to\infty}\sum_{j=1}^{2^n}(h_j^{(n)}-h_{j+1}^{(n)})
\frac{2j\delta}{2^n\sqrt{2\pi}}M\\
&=M\lim_{n\to\infty}\int_{-\infty}^\infty s_n(y)\,dm(y)
=M\int_{x-\delta}^{x+\delta}h_\lambda(x-y)\,dm(y)\\
&\leq M\int_{-\infty}^\infty h_\lambda(x-y)\,dm(y)=M,
\end{split}
\end{equation}
from which \eqref{8} follows. By \eqref{6}, $M\leq\epsilon$, so
\begin{equation}
\limsup_{\lambda\to 0}\int_{x-\delta}^{x+\delta}|f_x(y)|h_\lambda(x-y)\,dm(y)
\leq\epsilon.
\end{equation}
Since $\epsilon>0$ was arbitrary,
$$\limsup_{\lambda\to 0}\int_{x-\delta}^{x+\delta}|f_x(y)|h_\lambda(x-y)\,dm(y)=0,$$
whence
\begin{equation}\tag{9}\label{9}
\lim_{\lambda\to 0}\int_{x-\delta}^{x+\delta}|f_x(y)|h_\lambda(x-y)\,dm(y)=0.
\end{equation}
From \eqref{7} and \eqref{9}, we see that
\begin{equation}
\lim_{\lambda\to 0}\,(|f_x|*h_\lambda)(x)=
\lim_{\lambda\to 0}\int_{-\infty}^\infty|f_x(y)|h_\lambda(x-y)\,dm(y)=0.
\end{equation}
Finally,
\begin{equation}
\begin{split}
|(f_x*h_\lambda)(x)|
&=\Biggl|\int_{-\infty}^\infty f_x(y)h_\lambda(x-y)\,dm(y)\Biggr|\\
&\leq\int_{-\infty}^\infty|f_x(y)|h_\lambda(x-y)\,dm(y)
=(|f_x|*h_\lambda)(x),
\end{split}
\end{equation}
so $\lim_{\lambda\to 0}\,(f_x*h_\lambda)(x)=0$, as was required to prove.
Best Answer
This counter example that shows that the answer is no in general:
$$\lim_n\int f g_n\,dm = \int f\,dm\tag{1}\label{one}$$ One can appeal to uniform continuity to get for any $\varepsilon >0$, a $\delta>0$ such that $|f(x)-f(y)|<\varepsilon$ whenever $|x-y|<\delta$. Taking $n$ sufficiently large (e.g. $n>\frac{1}{\delta}$)
$$ \Big(f\big(\frac{k}{n}\big) -\varepsilon\Big)g_n\mathbb{1}_{\big[\tfrac{k}{n},\tfrac{k+1}{n}\big]}\leq f g_n\mathbb{1}_{\big[\tfrac{k}{n},\tfrac{k+1}{n}\big]}\leq \Big(\varepsilon + f\big(\frac{k}{n}\big)\Big)g_n\mathbb{1}_{\big[\tfrac{k}{n},\tfrac{k+1}{n}\big]} $$ Integration gives $$ \Big( f\big(\frac{k}{n}\big) -\varepsilon\Big)\frac{1}{n}\leq \int_{\big[\tfrac{k}{n},\tfrac{k+1}{n}\big]} f g_n\,dm \leq \Big(\varepsilon + f\big(\frac{k}{n}\big)\Big)\frac{1}{n} $$ Adding over all $0\leq k<n$ results in $$ \Big|\int_I fg_n\,dm -\sum^{n-1}_{k=0}f\big(\frac{k}{n}\big)\frac{1}{n}\Big|\leq\varepsilon $$ As $f$ is continuous, $\sum^{n-1}_{k=0}f\big(\frac{k}{n}\big)\frac{1}{n}\xrightarrow{n\rightarrow\infty}\int_If\,dm$ and the claim follows.