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

Question

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$

Answer 1

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.

Answer 2

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.