[proofplan]
The argument has three pieces. First, we recall the MRA structure: $V_j$ is an increasing chain of closed subspaces with $\bigcap V_j = \{0\}$, $\overline{\bigcup V_j} = L^2(\mathbb{R})$, the dilation property $f \in V_j \iff f(2 \cdot) \in V_{j+1}$, and translates of the father wavelet $\varphi(\cdot - k)$ form an orthonormal basis of $V_0$. Define the *detail spaces* $W_j := V_{j+1} \ominus V_j$ as the orthogonal complement of $V_j$ in $V_{j+1}$. The mother wavelet $\psi$ is constructed so that $\{\psi(\cdot - k) : k \in \mathbb{Z}\}$ is an orthonormal basis of $W_0$. Second, dilation transports this to: $\{\psi_{j,k} : k \in \mathbb{Z}\}$ is an ONB of $W_j$ for every $j \in \mathbb{Z}$. Third, the spaces $W_j$ are pairwise orthogonal and $L^2(\mathbb{R}) = \overline{\bigoplus_{j \in \mathbb{Z}} W_j}$, so the union of the level-$j$ orthonormal bases is an orthonormal basis of $L^2(\mathbb{R})$.
[/proofplan]
[step:Recall the MRA axioms and define the detail spaces $W_j$]
A [multiresolution analysis (MRA)](/page/Multiresolution%20Analysis) of $L^2(\mathbb{R})$ is a sequence of closed subspaces $\{V_j\}_{j \in \mathbb{Z}}$ of $L^2(\mathbb{R})$ together with a function $\varphi \in V_0$ — the father wavelet — satisfying
. $V_j \subseteq V_{j+1}$ for all $j \in \mathbb{Z}$ (nesting);
. $\bigcap_{j \in \mathbb{Z}} V_j = \{0\}$ (separation);
. $\overline{\bigcup_{j \in \mathbb{Z}} V_j} = L^2(\mathbb{R})$ (density);
. $f \in V_j$ if and only if the dilation $f(2\, \cdot) \in V_{j+1}$ (scale invariance);
. $\{\varphi(\cdot - k) : k \in \mathbb{Z}\}$ is an orthonormal system in $V_0$ whose closed linear span equals $V_0$.
For each $j \in \mathbb{Z}$, set the *detail space*
\begin{align*}
W_j := V_{j+1} \ominus V_j = \{ g \in V_{j+1} : (g, h)_{L^2} = 0 \text{ for all } h \in V_j\},
\end{align*}
the orthogonal complement of $V_j$ inside $V_{j+1}$. Each $W_j$ is a closed linear subspace of $V_{j+1}$, and we have the orthogonal direct sum
\begin{align*}
V_{j+1} = V_j \oplus W_j
\end{align*}
inside $L^2(\mathbb{R})$. The mother wavelet $\psi \in W_0$ is constructed (by the standard MRA construction via the high-pass filter) so that
. $\{\psi(\cdot - k) : k \in \mathbb{Z}\}$ is an orthonormal basis of $W_0$.
We take property (W) as part of the hypothesis that "$\psi$ is the mother wavelet associated with the MRA"; the construction of $\psi$ from $\varphi$ is the [Mother Wavelet Construction Theorem](/theorems/???), and we cite it as a named result. The hypothesis we verify when invoking it is the existence of the MRA and the associated low-pass filter $m_0$ — both granted.
[/step]
[step:Verify that $\{\psi_{j,k} : k \in \mathbb{Z}\}$ is an orthonormal basis of $W_j$ via dilation]
Define for $j \in \mathbb{Z}$ the dilation operator
\begin{align*}
D_j: L^2(\mathbb{R}) &\to L^2(\mathbb{R}), \\
f &\mapsto 2^{j/2}\, f(2^j \, \cdot).
\end{align*}
A direct change of variables (Lebesgue measure scaling: $d\mathcal{L}^1(2^j x) = 2^j\, d\mathcal{L}^1(x)$) gives
\begin{align*}
\|D_j f\|_{L^2}^2 = \int_{\mathbb{R}} 2^j |f(2^j x)|^2\, d\mathcal{L}^1(x) = \int_{\mathbb{R}} |f(y)|^2\, d\mathcal{L}^1(y) = \|f\|_{L^2}^2,
\end{align*}
so $D_j$ is an isometry. Since $D_{-j}$ is the inverse, $D_j$ is unitary.
By the MRA scale-invariance axiom (M4), $D_j$ maps $V_0$ onto $V_j$ and $V_1$ onto $V_{j+1}$, hence sends $W_0 = V_1 \ominus V_0$ unitarily onto $W_j = V_{j+1} \ominus V_j$. By property (W), $\{\psi(\cdot - k) : k \in \mathbb{Z}\}$ is an orthonormal basis of $W_0$. Applying the unitary $D_j$ — which preserves orthonormal bases —
\begin{align*}
D_j[\psi(\cdot - k)](x) = 2^{j/2}\, \psi(2^j x - k) = \psi_{j,k}(x).
\end{align*}
Therefore $\{\psi_{j,k} : k \in \mathbb{Z}\}$ is an orthonormal basis of $W_j$ for every $j \in \mathbb{Z}$.
[/step]
[step:Show that distinct detail spaces are orthogonal]
We claim that for all $j, j' \in \mathbb{Z}$ with $j \ne j'$,
\begin{align*}
W_j \perp W_{j'}.
\end{align*}
Without loss of generality $j < j'$. Then $W_j \subset V_{j+1} \subseteq V_{j'}$ (using (M1) iteratively), and by definition $W_{j'} = V_{j'+1} \ominus V_{j'}$, so every element of $W_{j'}$ is orthogonal to every element of $V_{j'}$, in particular to every element of $W_j$. The orthogonality $W_j \perp W_{j'}$ for $j \ne j'$ follows.
Combined with the previous step, the family
\begin{align*}
\bigcup_{j \in \mathbb{Z}} \{\psi_{j,k} : k \in \mathbb{Z}\} = \{\psi_{j,k} : j, k \in \mathbb{Z}\}
\end{align*}
consists of orthonormal vectors of $L^2(\mathbb{R})$: orthonormality within each level $W_j$ comes from step 2, orthogonality across levels from this step.
[/step]
[step:Prove completeness using the density and separation axioms (M2), (M3)]
It remains to show that the orthonormal system $\{\psi_{j,k} : j, k \in \mathbb{Z}\}$ is total in $L^2(\mathbb{R})$ — equivalently, that
\begin{align*}
L^2(\mathbb{R}) = \overline{\bigoplus_{j \in \mathbb{Z}} W_j},
\end{align*}
where the bar denotes closure and $\bigoplus$ is the algebraic orthogonal sum.
By a telescoping argument, for any $j_0 < j_1$ in $\mathbb{Z}$,
\begin{align*}
V_{j_1} = V_{j_0} \oplus W_{j_0} \oplus W_{j_0+1} \oplus \cdots \oplus W_{j_1 - 1},
\end{align*}
where each $\oplus$ is the orthogonal direct sum from step 1. This holds by repeated application of $V_{j+1} = V_j \oplus W_j$.
Letting $j_1 \to \infty$ first and then $j_0 \to -\infty$, we proceed in two stages.
Stage A. Density (M3): $\overline{\bigcup_{j \in \mathbb{Z}} V_j} = L^2(\mathbb{R})$. Let $P_{V_j}: L^2(\mathbb{R}) \to L^2(\mathbb{R})$ denote the orthogonal projection onto the closed subspace $V_j$. Density (M3) is equivalent to
\begin{align*}
\lim_{j \to \infty} P_{V_j} f = f \qquad \text{in } L^2(\mathbb{R}) \text{ for every } f \in L^2(\mathbb{R}).
\end{align*}
We verify this equivalence: the limit-projection statement says that the union $\bigcup V_j$ is dense in the orthogonal complement of $\{0\}$, which is all of $L^2$, so it is exactly (M3).
Stage B. Separation (M2): $\bigcap_{j \in \mathbb{Z}} V_j = \{0\}$, equivalent to
\begin{align*}
\lim_{j \to -\infty} P_{V_j} f = 0 \qquad \text{in } L^2(\mathbb{R}) \text{ for every } f \in L^2(\mathbb{R}).
\end{align*}
The equivalence: $P_{V_j} f$ is monotonically decreasing in norm as $j \to -\infty$ (since $V_{j-1} \subseteq V_j$ implies $\|P_{V_{j-1}}f\|_{L^2} \le \|P_{V_j}f\|_{L^2}$), so the strong limit exists and lies in $\bigcap_j V_j = \{0\}$ by (M2).
Combining: for any $f \in L^2(\mathbb{R})$ and any $j_0 < j_1$ in $\mathbb{Z}$, by the orthogonal decomposition $V_{j_1} = V_{j_0} \oplus \bigoplus_{j_0 \le j < j_1} W_j$ and Pythagoras,
\begin{align*}
P_{V_{j_1}} f = P_{V_{j_0}} f + \sum_{j = j_0}^{j_1 - 1} P_{W_j} f, \qquad \|P_{V_{j_1}}f - P_{V_{j_0}}f\|_{L^2}^2 = \sum_{j = j_0}^{j_1 - 1} \|P_{W_j} f\|_{L^2}^2.
\end{align*}
Letting $j_1 \to \infty$ (so $P_{V_{j_1}} f \to f$ in $L^2$ by Stage A) and $j_0 \to -\infty$ (so $P_{V_{j_0}}f \to 0$ in $L^2$ by Stage B),
\begin{align*}
\|f\|_{L^2}^2 = \lim_{j_0 \to -\infty} \lim_{j_1 \to \infty} \|P_{V_{j_1}}f - P_{V_{j_0}}f\|_{L^2}^2 = \sum_{j \in \mathbb{Z}} \|P_{W_j}f\|_{L^2}^2,
\end{align*}
and analogously $f = \sum_{j \in \mathbb{Z}} P_{W_j}f$ with $L^2$ convergence. The series converges absolutely in $L^2$ norm because the partial sums are Cauchy: for $-N_2 \le -N_1 < N_1 \le N_2$,
\begin{align*}
\Big\| \sum_{|j| > N_1, j \le N_2}\!\!\!\!\!\!\!\!P_{W_j}f \Big\|_{L^2}^2 = \!\!\sum_{|j| > N_1, j \le N_2}\!\!\!\!\!\!\|P_{W_j}f\|_{L^2}^2 \to 0
\end{align*}
as $N_1 \to \infty$ (by absolute convergence of the displayed Parseval-style sum).
Hence $L^2(\mathbb{R}) = \overline{\bigoplus_{j \in \mathbb{Z}} W_j}$. Since $\{\psi_{j,k} : k \in \mathbb{Z}\}$ is an ONB of $W_j$ (step 2) and the $W_j$ are mutually orthogonal (step 3), we conclude that $\{\psi_{j,k} : j, k \in \mathbb{Z}\}$ is an orthonormal basis of $L^2(\mathbb{R})$.
[/step]
[step:Derive the wavelet expansion and Parseval identity]
Let $f \in L^2(\mathbb{R})$. Since $\{\psi_{j,k}\}$ is an ONB, the abstract Hilbert-space ONB expansion yields
\begin{align*}
f = \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} \langle f, \psi_{j,k}\rangle_{L^2}\, \psi_{j,k}
\end{align*}
with convergence in $L^2(\mathbb{R})$, where $\langle f, g\rangle_{L^2} = \int_{\mathbb{R}} f\bar g\, d\mathcal{L}^1$ denotes the $L^2$ inner product. Parseval's identity is the corresponding norm equality
\begin{align*}
\|f\|_{L^2}^2 = \sum_{j \in \mathbb{Z}} \sum_{k \in \mathbb{Z}} \big| \langle f, \psi_{j,k}\rangle_{L^2}\big|^2,
\end{align*}
which is again the standard ONB Parseval. This completes the proof.
[/step]