Here is a quick way to compute it in the case of an odd sphere (so even number of variables) and $p$ is an odd prime. Note that for a general odd number its easy to reduce to this case by the Chinese reminder theorem and an easy lifting from $p$ to $p^m$. The equation $x^2+y^2=a$ is equivalent to the equation $N_{k[i]/k}(\alpha)=a$ where $k=\mathbb{F}_p$ is the field with $p$-elements.
If $p\equiv 3 mod(4)$ then this equation has exactly $p+1$ solutions if
$a \ne 0$ and one solution if $a=0$ because the Norm map $Nm:k[i]^\times\to k^\times$ is a surjective homomorphism.

If we denote by $f_\ell(a)$ the function that count the solution with $\ell$ variables, then byt clearly we have

$$f_{2\ell}(a)=\sum_{a_1+...+a_\ell = a}f_2(a_1)f_2(a_2)...f_2(a_\ell)$$,
because we can count conditionally on the value of $x_{2i-1}^2+x_{2i}^2$ for every $1\le i\le \ell$ and sum the results. In other words, $f_{2\ell}$ is the $\ell$-fold convolution $f_2*f_2*...*f_2$. Let $I$ denote the constant function with value 1 on $k$, and $\delta_0$ the function which is $1$ at $0$ and $0$ otherwise. The calculation we did for $f_2$ gives exactly that
$f_2=(p+1)I -p\delta_0$. Using the linearity of the convolution and the convolutions $I*I=pI$, $\delta_0*f=f$ for every $f$ we get that
$$((p+1)I - p\delta_0)^{*\ell}=\sum_{k=0}^{\ell-1} {\ell \choose k} (-1)^{\ell - k}(p+1)^kp^{k-1}p^{\ell-k}I+(-1)^\ell p^\ell \delta_0=$$

$$=p^{\ell-1} \sum_{k=1}^{\ell} {\ell \choose k} (-1)^{\ell - k}(p+1)^kI+(-1)^\ell p^\ell\delta_0=p^{\ell-1}(p^\ell-(-1)^\ell)I+(-1)^\ell p^\ell \delta_0$$

and this (if I don't have any computation mistake!) gives you the answer in the case $p\equiv3 mod 4$. A similar argument works for $p\equiv 1 mod 4$ with different $f_2$ (here the algebra $k[i]$ splits into $k\oplus k$)

For even spheres, its slightly more complicated since eventually you have to convolve with the function $f_1$ which is the sum of $I$ and a quadratic character, but this is not a big difference when you notice that if $\chi$ is such a character then $\chi*I=0$ and $(I+\chi)*(I+\chi)=I+2\chi+\chi*\chi=f_2$
so you have all the convolution formulas you need.