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;
}