# Efficient way to calculate $n$ first elements of cross-correlation using FFT

The cross-correlation between two signals can be calculated using $$\texttt{IFFT}\left(\texttt{FFT}(a) \cdot \overline{\texttt{FFT}(b)}\right)$$ for two binary signals $$a$$ and $$b$$.
How could one speed up the calculation if you want only the first $$n$$ elements of the cross-correlation? In my case $$n$$ is much smaller than the length of the signals, with $$n \approx 5000$$, while the size of the signal might be $$10^9$$

• One way would be to decimate your signals before doing the FFT processing. Do you really need all $10^9$ samples?
– Jdip
Dec 1, 2022 at 14:33
• My problem is that the signal is very noisy and sparse, so using less than $10^9$ samples, I might not be able to find correlations Dec 1, 2022 at 14:35
• You can resort to downsampling if your signal is not wideband. Dec 1, 2022 at 17:47
• @endolith Aaah I see what you mean now. Then yes, I assume he meant the middle $n$ elements. Or circular cross-correlation. I apologize for the confusion, should have realized that's what you were saying the whole time ;)
– Jdip
Dec 1, 2022 at 22:47
• @Jdip Yeah, nevermind, I get it now. "First n elements of the cross-correlation" aren't likely useful unless you're thinking in terms of positive and negative lag and "first" is starting at 0 lag. Dec 1, 2022 at 23:34

For long sequences, the FFT approach really dominates for "fast correlation" (and fast convolution when you don't conjugate). For very small number of samples $$n$$, a direct cross-correlation using dot products would be more efficient, but for the OP's case this is much small than $$n=5000$$. Reducing the FFT block size is the common approach to improve efficiency, but this is at the expense of the resolution bandwidth. Below I determine the cross-over where it would be more efficient to simply compute the cross-correlation directly using dot products to compute the cross-correlation without any loss in resolution bandwidth:

Ultimately the time duration needed is driven by the resolution bandwidth which is $$1/T$$ for a non-windowed duration of time $$T$$. If the correlations desired are due to single tones in the presence of white noise, the noise density relative to the weakest signal level will dictate the minimum length desired for which to reduce the noise below the signal level to be detected. The number of samples over this time is driven by the highest frequency of interest, and therefore sets the sampling rate according to the Nyquist criteria. If the waveform is significantly oversampled, then processing can be significantly reduced by properly decimating the waveform first.

Assuming the OP has already determined that there must be $$10^9$$ samples, we can then explore what the cross-over would be where the direct cross-correlation using repeated dot products (for each time offset) would be more efficient than the FFT computation.

The cross correlation approach described using FFT's is a circular cross-correlation. Given $$N$$ total samples, We could compute a single element directly using $$N$$ real multiplies and $$N$$ real adds (assuming $$a$$ and $$b$$ are real, or 4x more multiples and 2x more adds if $$a$$ and $$b$$ are complex) by using a dot product.

With the FFT approach which computes this dot product for every (circular) offset in the time domain, there are approximately $$2N \log_2(N)$$ real multiplications and $$2N \log_2(N)$$ real additions for every FFT or IFFT computation. So to compute the cross correlation using FFTs, we have three total FFT or IFFT computations, in addition to a product of two of the FFTs requiring $$4N$$ more real multiplications. If we total this up, we get $$3 \times 2N \log_2(N) + 4N$$ real multiplications and $$3 \times 2N \log_2(N)$$ real additions.

From this, the cross-over number of samples $$n$$ below which it would be more efficient to simply use a dot product for computing the cross correlation for the n samples is as follows:

For $$a$$, $$b$$, real:

$$n$$ for same number of multiplications required:

$$n = \frac{6N \log_2(N) + 4N}{N} = 6 \log_2(N)+4$$

$$n$$ for same number of additions required:

$$n = \frac{6N \log_2(N)}{N} = 6 \log_2(N)$$

In the OP's case, $$N= 10^9$$, and for this $$n$$ as limited by multipliers is $$6 \log_2(10^9)+4 \approx 183$$

Meaning for this case, using the FFT for computing the (circular) cross-correlation is more efficient than manually computing a dot product for each offset when $$n>183$$. For the case that $$n \ll 183$$, a dot product for each sample would be increasingly more efficient. If $$a$$ and $$b$$ are complex, then this crossover point reduces by 4, or only 45 samples!

An additional word of caution when correlating over extremely long sequences. At some point, inevitably, the signal of interest will be insufficiently stable and the noise will no longer be stationary, at which point if either of these occur and we continue to correlate beyond this time, the SNR of our desired result will actually degrade rather than improve. There are statistical techniques (such as the Allan Deviation) to test for optimum averaging times for this very reason.

• Can you show your calculations in more detail? $3 \times 2N log_2(N) + 4N$ with $N=10^9$ would be 183e9 multiplications, while the direct method would only be $n(n+1)/2$ = 12e6 multiplications, no? Dec 1, 2022 at 22:45
• Nevermind I get it now. There's no reason you'd want the "first 5000 outputs" in terms of the actual cross-correlation output. OP must mean "the first 5000 outputs" starting at full overlap and going in the positive lag direction, in which case it's 1e9 multiplies for each output ≈ 5000e9 which is obviously higher. Dec 1, 2022 at 23:33

Do a blockwise cross-correlation:

1. Take samples [0, 10000) from signal a.
2. Take samples [2500, 7500) from signal b, zero-pad on both sides.
3. Calculate the cyclic cross-correlation using FFTs, transform length is 10000. Store the result.
4. Repeat steps 1-3 for all further blocks, stride is 5000. (The second block uses samples [5000, 15000) from a, [7500, 12500) from b)
5. Average all block results.

You end up with a cross-correlation of 10000 samples, corresponding to shifts [-5000, +5000), of which only the values for [-2500, +2500] are valid. The others are spoilt from wrap-around of the cyclic block convolutions.

• Is that faster than an FFT based cross correlation with full input signals?
– Jdip
Dec 2, 2022 at 1:44
• @Jdip - Much, much faster. Only linear runtime complexity with respect to input size. (N log K for a chunked convolution with chunk size K, compared to N log N for the full convolution). My algorithm essentially overlap-save adapted for cross-correlation. Dec 2, 2022 at 9:25
• But wouldn't you then need to do these N log K computations about N/K times?
– Jdip
Dec 2, 2022 at 9:29
• Ah right those are K-length chunks.. To be a little more exact, I see 6KlogK+4K with N/K of them, so 6NlogK + 4 (instead of 6NlogN+4 for the full one). So with the OP's specs, about 2.5 times faster. If that gives the same result as the full cross-correlation with K lags output, I'm sold! Will try a little simulation later. I wonder what @danboschen thinks about this.
– Jdip
Dec 2, 2022 at 9:39
• @Jdip - Wikipedia has good articles on overlap-add and overlap-save. These methods are well known and so is the N log K runtime complexity for the entire thing. Normally, we use those methods to convolve a long signal (length N) with a short one (length K, typically a filter kernel) and get a long (length N) result. I used the same technique here with two long inputs and a short output. Dec 2, 2022 at 9:45