Time-domain cross-correlation with padded signals [closed]

Question

I have a frequency domain cross-correlation implementation written in C (based on: https://github.com/dMaggot/libxcorr).

It uses the library FFTW3 and this is the gist of it:

//2 signals of size N, and an output array of size 2*N-1
signala = [0, 0, 0, ..., signal_a_orig] //padded to 2*N
signalb = [signal_b_orig, 0, 0, 0, ...] //padded to 2*N

fftw_plan pa = fftw_plan_dft_1d(2 * N - 1, signala, outa, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan pb = fftw_plan_dft_1d(2 * N - 1, signalb, outb, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan px = fftw_plan_dft_1d(2 * N - 1, out, result, FFTW_BACKWARD, FFTW_ESTIMATE);

fftw_execute(pa);
fftw_execute(pb);

fftw_complex scale = 1.0/(2 * N -1);
for (int i = 0; i < 2 * N - 1; ++i)
    out[i] = outa[i] * conj(outb[i]) * scale;

fftw_execute(px);
//result contains the result

Here are the outputs of an example:

N = 10
for (int i = 0.0; i < N; ++i) {
    in1[i] = 6.0*i + (7.5*i + 13.5)*I;
    in2[i] = -5.0*i + (14.3*i)*I;
}

xcorr(in1, in2, out, N);

Here are the outputs printed of the result of this cross-correlation.

1737.450000     -607.500000
3977.100000     -2257.200000
6641.700000     -4825.800000
9654.000000     -8190.000000
12936.750000    -12226.500000
16412.700000    -16812.000000
20004.600000    -21823.200000
23635.200000    -27136.800000
27227.250000    -32629.500000
30703.500000    -38178.000000
25489.800000    -32022.000000
20546.400000    -26056.800000
15950.550000    -20405.700000
11779.500000    -15192.000000
8110.500000     -10539.000000
5020.800000     -6570.000000
2587.650000     -3408.300000
888.300000      -1177.200000
0.000000        0.000000

On the other hand, I have the initial workings of a time-domain autocorrelation piece of code that I took from here: https://github.com/JorenSix/TarsosDSP/blob/master/src/core/be/tarsos/dsp/pitch/McLeodPitchMethod.java

Here's what it looks like in my C code:

//original signals, not padded. size = N = 10
for (int tau = 0; tau < N; ++tau) {
    complex acf = 0 + 0*I;
    for (int i = 0; i < N - tau; ++i) {
         acf += conj(signala[i]) * signalb[i+tau];
    }
    result[N - 1 - tau] = creal(acf) - I*cimag(acf);
}

I got acf += conj(signala[i]) * signalb[i+tau] from this Wikipedia article describing discrete time cross correlation: https://en.wikipedia.org/wiki/Cross-correlation

The part I don't understand is this one: result[N - 1 - tau] = creal(acf) - I*cimag(acf). Honestly, I figured this out by trial and error to make my time domain code give me the same results as the frequency domain code. I can't find any resources online that say why.

Also, this only gives me half a signal compared to the double-padded FFT implementation.

The results are good:

1737.450000     -607.500000
3977.100000     -2257.200000
6641.700000     -4825.800000
9654.000000     -8190.000000
12936.750000    -12226.500000
16412.700000    -16812.000000
20004.600000    -21823.200000
23635.200000    -27136.800000
27227.250000    -32629.500000
30703.500000    -38178.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000
0.000000        0.000000

The first 10 values are identical to the FFT implementation, and, of course, since I don't use padded signals, I'm only producing N values instead of 2*N-1.

Now, I can't figure out how to produce the second half of the values.

Simply bumping up the size limits of the time-domain loops and using the padded signals instead of the original signals doesn't give good results.

Is there a way I can get identical outputs from the FFT and time-domain implementations? Should I just take half of the output into consideration? results[0:N/2] is identical for both implementations.

Thanks in advance.

edit Here's the modified time-domain code based on the best answer:

-       for (int tau = 0; tau < N; ++tau) {                                                                                        
+       for (int tau = 0; tau < N2; ++tau) {                                                                                       
                acf = 0 + 0*I;                                                                                                     
-               for (int i = 0; i < N - tau; ++i) {                                                                                
-                       acf += conj(signala[i]) * signalb[i+tau];                                                                  
+               for (int i = 0; i < N2; ++i) {                                                                                     
+                       acf += signala_ext[(i+tau) % N2] * conj(signalb_ext[i % N2]);                                              
                }                                                                                                                  
-               result[N - 1 - tau] = creal(acf) - I*cimag(acf);                                                                   
+               result[tau] = acf;                                                                                                 
        }              

Answer 1

You can consider the code below, which is written in Python. But, you should get the idea of how to do it in the time domain (but note that it is not very efficient). I have not implmented a possible optimization, since the signals are padded with zeros. So, the actual summation can be done in a shorter fashion.

Note that the approach with fft/ifft by definition performs actually a circulant correlation/convolution. But, by padding the signals with N-1 zeros, you get the same result as if you did a linear convolution. I have also done a circulant correlation in the time domain, by using the modulo-operation. You could do a linear one by not adding the values, if one index is out-of-bounds (or you adjust the summation boundaries).

N = 10
i = np.arange(N)
Z = np.zeros(N-1)
inA = 6*i + (7.5*i+13.5)*1j
inB = -5*i + (14.3*i)*1j

padA = np.hstack([Z, inA])
padB = np.hstack([inB, Z])

N2 = 2*N-1
scale = 1/N2

outA = np.fft.fft(padA)
outB = np.fft.fft(padB)

out = outA * outB.conj() * scale

result = N2*np.fft.ifft(out)

acorr2 = np.zeros(N2, dtype=complex)
acorr3 = np.correlate(inA, inB, 'full')  # Using the standard python correlate function
for tau in range(N2):
    acf = 0j
    for i in range(N2):
        acf = acf + (padA[(i+tau) % N2]) * np.conj(padB[(i) % N2])
    acorr2[tau] = acf  # calculating it by hand
plt.subplot(121)
plt.plot(np.imag(result))
plt.plot(np.imag(acorr2), 'r-x')
plt.plot(np.imag(acorr3), 'g-o')
plt.title('imag part')
plt.subplot(122)
plt.plot(np.real(result))
plt.plot(np.real(acorr2), 'r-x')
plt.plot(np.real(acorr3), 'g-o')
plt.title('real part')
plt.tight_layout()