# Counting the combinations of $i+j-k=\omega$ in a triple sum for any $\omega$

I have a function defined as the summation shown below, where the results contains multiple orthogonal terms at frequencies of $\omega=i+j-k$.

$s(t)=\sum_{i=1}^N\sum_{j=1}^N\sum_{k=1}^Ncos(i+j-k)t$

I would like to know if there is a direct way to determine, for any $\omega$, how many combinations of $i$, $j$, and $k$ in this triple-sum will have the same frequency $\omega=i+j-k$. Also, any answer that points to a process for generalizing to other similar problems would be welcome.

This is part of a more general problem where the expression is like this:

$s(t)=\sum_{i \in W}\sum_{j \in X}\sum_{k \in Z}(A_iA_jA_k)^2cos(i+j-k)t$

where $W$, $X$ and $Z$ are subsets from the range $[1,N]$.

This problem is related to calculating the third-order response of a system in the frequency domain.

With $1\le i,j,k \le N$ one can count the number of sums $i+j-k=\omega$ using polynomial multiplication. Note that also $-N \le -k \le -1$ which will help avoid the subtraction. The values possible for $\omega$ are from $2-N$ to $2N-1$ inclusive.
To reframe it with nonnegative variables let $x=i-1,y=j-1,z=(-k)-(-N)$ so that now each of $x,y,z$ ranges from $0$ to $N-1$. And the relation $i+j-k=\omega$ becomes $$x+y+z=r=\omega+N-2,$$ where now $r$ ranges from $0$ to $3N-3$.
Then the number of solutions to $x+y+z=r$ is the coefficient of $t^r$ in the expansion of $$(1+t+t^2+\cdots+t^{N-1})^3.$$ Looking on line I didn't find much in the way of a simple "closed form" formula for this coefficient. However there were some fast implementations of polynomial multiplication referred to. If it makes coding any easier the sum being cubed is the same as $$\frac{t^N-1}{t-1}.$$