See if this construction works (I think it does).
Let $\{x_k\}$ be an enumeration of the points where $f$ has a jump discontinuity. For each $k$, let $\delta_k = f(x_k) - \frac12\bigl(f(x_k+) + f(x_k-)\bigr)$, so that $\delta_k$ is the amount by which $f(x_k)$ deviates from the midpoint of the jump. Notice that $\sum|\delta_k|$ converges, because that sum must be less than $f(1)-f(0)$ (the total variation of $f$ over the interval).
For $n=1,\ 2,\ 3,\ \ldots$, define $$g_n(x) = \frac n2\int_{x-(1/n)}^{x+(1/n)}f(t)\,dt$$ and $$h_n(x) = \sum_{\{k\;:\;x-(1/n) \leqslant x_k \leqslant x+(1/n)\}}\bigl(1- n|x-x_k|\bigr)\delta_k$$, and let $f_n(x) = g_n(x) + h_n(x).$ (To cope with integrals that extend beyond the ends of the unit interval, define $f(t) = f(0)$ when $t<0$, and $f(t) = f(1)$ when $t>1$.)
Then $g_n$ is continuous (because it is an integral), $h_n$ is continuous (because it is a uniformly convergent sum of continuous functions), and therefore $f_n$ is continuous. The big question is whether $f_n$ converges pointwise to $f$. The reason I think it does is that $g_n(x)$ represents the mean value of $f$ in the interval $[x-(1/n), \,x+(1/n)]$. As $n\to\infty$, that will converge to $f(x)$ at all points of continuity. At a jump, it will converge to the midpoint of the jump. As for $h_n(x)$, if $f$ is continuous at $x$ then for $n$ large enough there will be no large jumps in the interval $[x-(1/n), \,x+(1/n)]$, and therefore $h_n(x)$ will be small. On the other hand, if $x=x_k$ for some $k$ then, again for $n$ large enough, the sum of all the other jumps in the interval $[x_k-(1/n), \,x_k+(1/n)]$ will be much smaller than the jump at $x_k$, and so $h_n(x_k)$ will be close to $\delta_k$. Thus $f_n(x_k)$ will be close to $\frac12\bigl(f(x_k+) + f(x_k-)\bigr) + \delta_k = f(x_k).$
I don't have the energy to put all the $\varepsilon$s and $\delta$s into that argument to make it stand up properly, but I hope that it points in the right direction.