Since visoft requested i post the solution, here goes. It has been a while so not all the details are still fresh so feel free to ask anything if things are unclear!
Let's say we have an image $I$ and apply a filter, which is given by convolution by some $F$. The pixel values of $I$ convolved with $F$ are all the output of the inner product of $I$ with a translated (and less importantly: mirrored) version of $F$.
What i missed at the time was that each pixel should be seen as a coefficient of the new representation of the image. So instead of the 16 subbands i mentioned, the number of coefficients of our new representation is 16 times the number of pixels in the image.
Some notation: for pixel number $(a,b)$, let us write $F_{(a,b)}$ for $F$ translated $(a,b)$ steps (and mirrored) and $(I * F)(a,b)$ for the $(a,b)$-th pixel of $I$ convolved with $F$.
Also, let us call all the translated versions of $F$ our basisfunctions. One can check from the definition of convolution that if we convolve our filter output again with $F$, this results in multiplying every basisfunction $F_{(a,b)}$ by the coefficient $(I * F)(a,b)$ and adding them up. This results in
\begin{align*}
I * F * F &= \sum_{(a,b)} (I * F)(a,b) \cdot F_{(a,b)}\\
&= \sum_{(a,b)} \langle I, F_{(a,b)}\rangle \cdot F_{(a,b)}
\end{align*}
where the last equality follows by definition of convolution.
From this, one easily recognizes something that comes up using an orthogonal basis. Now the collection of $F_{(a,b)}$ for all $(a,b)$ is called a tight frame if the above mathematical expression is equal to $I$ itself! It is thus very similar to an orthogonal basis, except that one does not need the set $\{F_{(a,b)}\}_{(a,b)}$ to be linearly independent (if you are very sharp you will note that in this particular case it most likely will be linearly independent, but see the following).
In the paper, actually about 16 of these $F's$ are used. Let's number them $F_1$ up to $F_n$. This is called a tight frame if the following holds:
$$
I = \sum_n \sum_{(a,b)} \langle I, F_{n,(a,b)} \rangle \cdot F_{n,(a,b)}
$$
This expression is equal to (check this):
$$
I = \sum_n (I * F_n * F_n).
$$
Since convolution is equivalent to multiplication in the fourier domain, one can check this by going to the fourier domain, squaring all the $F_n$ pointwise and adding them up. That's simple to do. If that adds to an image that is equal to 1 in every pixel, we are dealing with a tight frame.
So concluding, i missed two things: convolving with one $F$ gives as much coefficients as there are pixels in the image, and convolving twice actually multiplies the basisfunctions with the coefficients and adds them up.
As I said, I wrote this quickly and it has been a while, so if anything is unclear feel free to ask again!
