  Recovering a signal (function) from finitely many binary or Fourier samples
is one of the core problems in modern medical imaging, and by now there exist a
plethora of methods for recovering a signal from such samples. Examples of
methods, which can utilise wavelet reconstruction, include generalised
sampling, infinite-dimensional compressive sensing, the
parameterised-background data-weak (PBDW) method etc. However, for any of these
methods to be applied in practice, accurate and fast modelling of an $N \times
M$ section of the infinite-dimensional change-of-basis matrix between the
sampling basis (Fourier or Walsh-Hadamard samples) and the wavelet
reconstruction basis is paramount. In this work, we derive an algorithm, which
bypasses the $NM$ storage requirement and the $\mathcal{O}(NM)$ computational
cost of matrix-vector multiplication with this matrix when using Walsh-Hadamard
samples and wavelet reconstruction. The proposed algorithm computes the
matrix-vector multiplication in $\mathcal{O}(N\log N)$ operations and has a
storage requirement of $\mathcal{O}(2^q)$, where $N=2^{dq} M$, (usually $q \in
\{1,2\}$) and $d=1,2$ is the dimension. As matrix-vector multiplications is the
computational bottleneck for iterative algorithms used by the mentioned
reconstruction methods, the proposed algorithm speeds up the reconstruction of
wavelet coefficients from Walsh-Hadamard samples considerably.
