The explicit formula of the Riesz transform of the box spline B2(x, y) : = B2(x)B2(y) is given, where B2 is the cardinal B‐spline of order 2. By using the wavelet multilevel method, a sampling recovery scheme derived from B2 is established to recover the Riesz transform of the functions in Sobolev space
with s > 1. For any fixed level, our recovery is involved with a finite sum series. Since the Riesz transforms of some functions are continuous but
has numerical singularity at some points, it is necessary to eliminate the numerical singularity. We first establish the shift‐perturbation error estimate of the multilevel sampling approximation, derived from B2, to the functions in
. By the perturbed approximation system, we give a method to eliminate the numerical singularity. Numerical simulations are conducted to test the recovery efficiency.