Series Sum Operator Properties and Identities
That equation is almost correct. Tim gave the correct solution in a comment, which is just a scaled version of the given formula. I'd like to show how you can quite easily see this.
First, note that
$$\sum_{m=-N}^N\phi(n-m)=(r\star\phi)(n)\tag{1}$$
where $\star$ denotes convolution, and $r(n)$ is a sequence of ones in the range $n\in[-N,N]$. Second, summing a sequence is also a convolution with the sequence $r(n)$ evaluated at $n=0$. So we can write
$$\begin{align}\sum_{n=-N}^N\sum_{m=-N}^N\phi(n-m)&amp;=(r\star r\star\phi)(n)\Big|_{n=0}\\&amp;=(s\star\phi)(n)\Big|_{n=0}\tag{2}\end{align}$$
where $s(n)$ is the convolution of two rectangular sequences, which is a triangular sequence:
$$s(n)=(2N+1)\left[1-\frac{|n|}{2N+1}\right],\qquad n\in[-2N,2N]\tag{3}$$
Combining $(2)$ and $(3)$ gives the desired result, which equals the formula in your question, scaled by $2N+1$.