It seems to be true. Due to the atomic decomposition, isit is enough to approximate in $H^1$ any "$(1,\infty)$"-atom, that is, a function $a$ such that $$\mathrm{supp}\, a \subset Q,\quad \|a\|_{L^\infty}\leq |Q|^{-1}, \quad \int_Q a=0$$ for a cube $Q$.
In order to do this, we can cut $Q$ into $N=k^n$ smaller equal cubes $Q_1, Q_2, \ldots, Q_N$ for a large $k$. We put $$ b_N(x)=\sum_{j=1}^N \Big( \frac{1}{|Q_j|}\int_{Q_j}a(y)\,dy\Big) \chi_{Q_j}(x). $$ Put $g_N=a-b_N$. Due to Lebesgue differentiation theorem we have $$ \lim\limits_{N\to \infty}\int |g_N(x)|^2dx = 0. $$
Then for big $N$ the following holds: $$ \mathrm{supp}\, g_N\subset Q, \quad \|g_N\|_{L^2}\leq \varepsilon |Q|^{-1/2}, \quad \int_Q g_N = 0. $$
It means that $\varepsilon^{-1} g_N$ is a $(1,2)$-atom. Hence, $\|g_N\|_{H^1}\leq C\varepsilon$.
The atomic decomposition in $H^1$ is described in the book "Weighted norm inequalities and related topics" by Garcia-Querva and Rubio de Francia. Another reference is "Modern Fourier Analysis" by Grafakos.