From a mathematical point of view (it means without any connection to a practical use described in the summary), the authors consider spherical means defined by \[ f:\mathbb R^d \rightarrow \mathbb R,\quad Mf(y,r)=\frac 1 {\omega_{d-1}} \int_{S^{d-1}}f(y+r\xi)d\sigma(\xi) \] (centerpoint \(y\in \mathbb R^d\), radius \(r>0\), \(\sigma\) denotes the surface measure on the sphere and \(\omega_{d-1}=\sigma(S^{d-1})\)) and use iterative methods for their computation. A spectral discretization via trigonometric polynomials in such a way that computations can be done via a nonequispaced fast Fourier transform (FFT), however, the sparse FFT is used in the three-dimensional case. The studied theory is illustrated numerically.