I have a sequence of fixed length $N$ input arrays ($x_0$, $x_1$, ..., $x_{n-1}$), and my goal is to compute output arrays as such:
$$d_{00}[l] = \sum_{i = 0}^{N - l} x_0[i]x_0[i + l]$$ $$d_{01}[l] = \sum_{i = 0}^{N - l} x_0[i]x_1[i + l]$$ $$d_{02}[l] = \sum_{i = 0}^{N - l} x_0[i]x_2[i + l]$$ $$...$$ $$d_{10}[l] = \sum_{i = 0}^{N - 1} x_1[i]x_0[i + l]$$ $$...$$ With $l$ having a maximum value that scales linearly in $N$. That is to say, $d_{ij}$ will contain, in it's $k^{th}$ entry, the sum of all pairs in $x_i$ and $x_j$ that are distance $k$ away.
My understanding - Clearly, since each $d_{ij}$ output requires $l \propto N$ entries, each containing $N$ multiplies, the time complexity of producing each $d_{ij}$ output array is $O(N^2)$. These pairwise multiplies hidden behind a summation seem optimal for a Fourier Transform to improve the time complexity to $O(N \log N)$. However, I understand the standard Fourier Transform convolution does not allow for limiting to constant distance, instead allowing easy sum of pairwise indices in the form $\sum_{i=0}^{c} x[i] y[c - i]$.
How can I use this standard convolution, or any variation of the Fourier transform algorithm to make this problem computationally tractable?