Delay/Lag calculation using cross correlation in frequency domain

Question

I need to compute cross correlation of two signals using frequency domain. Hence, I am using the below equation for calculation.

 max(IFFT(FFT(a)*conj(FFT(b))))

Using this equation, how to find the lag? (lag similar to matlab built-in function of xcorr) (I have observed that when cross correlated in frequency domain, the max value is always present in the first index of the output array.) Any basic ideas will be helpful for me to learn and apply it.

Answer 1

It is incorrect that the first sample will always be maximum. The lag is simply determined by the index of the maximum of the absolute value of the result.

Here is a simple example to illustrate the proper operation and result:

x = randn(512,1);
y = circshift(x, 20);    # y is a shifted version of x
out = ifft(fft(y).*conj(fft(x));
plot(abs(y))
[peak, lagp1] = max(abs(out))  # lag is one less than lagp1  

Doing the frequency domain (circular correlation) is great for this application when the lag is sufficiently smaller than the entire data block captured. Unlike a linear convolution, it is not affected by a DC offset that is the same on each signal (the entire result shifts up) and computes the result efficiently. Linear correlation is a good choice for applications when we are searching for instances of a waveform of relatively small time duration within a much larger block of data. For this we can use the same approach by zero padding the shorter data set out to the length of the larger block, or do the correlation in the time domain as an FIR filter. When the data sets are of similar duration in time, there is no need for zero padding or linear correlation.

Answer 2

You can inspect the MATLAB in-built function, try typing edit xcorr.m at the MATLAB command prompt. That should give you an idea on how things should be set up (lags, how to zero-pad your inputs, etc).

Answer 3

You do have the basic idea as to how to implement a frequency-domain cross correlation calculation. Now, as Hilmar stated in their comment, frequency domain multiplication ends up with circular convolution (or correlation if you complex-conjugate one of the sequences). Thus, you'll have to zero-pad your sequences accordingly.

From time-domain cross-correlation calculations we know that the resulting sequence is $L = M + N - 1$ long, where $M$ is the length of the first sequence and $N$ the length of the second. This means that you should zero-pad your sequences accordingly, so that both will end up with that length. This effectively can be done like (in MATLAB/Octave code)

% Create the sequences
x1 = rand(100, 1);
x2 = rand(80, 1);

% Get the lengths
x1Len = length(x1);
x2Len = length(x2);
L = x1Len + x2Len - 1;

% Zero pad
xzp1 = [x1; zeros(L - x1Len, 1)];
xzp2 = [x2; zeros(L - x2Len, 1)];

Now, you can use the calculations you have presented to get the cross correlation function. Please keep in mind that you may have to shift the result. For example, in MATLAB/Octave the negative frequencies are placed after the positive frequencies from $\frac{f_{s}}{2}$ to $f_{s}$, where $f_{s}$ is the sampling frequency. This results in the cross correlation function being circularly shifted by half the length of the whole window. In MATLAB/Octave you can use fftshift() to perform this operation on the result of your ifft() function.

Alternatively, you could calculate the cross correlation at the lags of interest with the following formula

$$ r_{x_{1}, x_{2}} \left(\tau \right) = \sum_{f = -\frac{f_{s}}{2}}^{\frac{f_{s}}{2}} \mathcal{R} \left\{ \mathbf{R}_{x_{1}, x_{2}} e^{-j 2 \pi f \tau} \right\}$$

where $\mathcal{R}$ denotes the real part, $\mathbf{R}_{x_{1}, x_{2}}$ is the cross spectrum, $f$ the frequency (you will use the frequencies corresponding to the central frequency of each bin in your FFT) and $\tau$ the lag of interest. The cross spectrum is given by

$$ \mathbf{R}_{x_{1}, x_{2}} = X_{1} \overline{X_{2}} $$

where $X_{1}$ and $X_{2}$ are the spectra of the time-domain sequences and $\overline{\left[ ~ \cdot ~ \right]}$ denotes complex conjugation.