The pseudo-code for cross-correlation implements the naive method, which is O(M*N), M and N being the lengths of the arrays being correlated. I think it might just have been used to illustrate, but just in case, the proper way of computing the cross-correlation is with the FFT, as that is "linearithmic" in time. It's commonly called "multiplication in the frequency domain" or just frequency domain convolution. Basically, take the FFT of both input sequences, multiply them, and take the inverse FFT of the product.
Original author here. They do this already, the pseudo function 'correlate' operates in the frequency domain. It doesn't negate the fact that you still correlate the signals of different stations pairwise.
If you're curious, the basis of this algorithm is the convolution theorem: https://en.wikipedia.org/wiki/Convolution_theorem
Here is an example implementation of cross-correlation using the FFT, with NumPy:
http://www.eliteraspberries.com/blog/2013/08/application-of-...