## Cross Correlation Using FFTW3

As part of my thesis project, I needed to get the Cross Correlation of two signal functions as implemented in Matlab, but using anything else available in C/C++. The obvious way was to implement the definition of Discrete Cross Correlation manually, which is an

$O\left({n}^{2}\right)$

algorithm, yet that's a bad idea when you have a formula that puts Cross Correlation in terms of the Fourier transform, since you can then take advantage of the Fast Fourier Transform (which, by Cooley Tukey algorithm is

$O\left(n\mathrm{log}n\right)$

).

Then this post is about how I got Cross Correlation on C using FFTW3, since I had to spend quite some time compiling info to get it. Furthermore, the goal of the following code was to imitate the exact output of Matlab's xcorr function, which uses a slightly different formula from that on the previous link (namely, using the dot product of the first Fourier transform by the conjugate of the second one, that is

${ℱ}^{-1}\left\{ℱ\left\{𝕒\right\}·\overline{ℱ\left\{𝕓\right\}}\right\}$

for input signals

$𝕒$

and

$𝕓$

).

Step 1: Use C's complex numbers. It'll make your life easier when dealing with conjugates (and printing stuff, if you want to do so). To do that, include complex.h header before fftw3.h:

[code lang="c"]
#include
#include
[/code]

Step 2: Create "extended" versions of the input signals. Input signals will be of a given size, say

$N$

, but you'll need to zero-pad those to a length of

$2N-1$

so that they fit on the output signal during the transformation calculations:

[code lang="c"]
fftw_complex * signala_ext = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
fftw_complex * signalb_ext = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
fftw_complex * out_shifted = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
fftw_complex * outa = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
fftw_complex * outb = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
fftw_complex * out = (fftw_complex *) fftw_malloc(sizeof(fftw_complex) * (2 * N - 1));
[/code]

Step 3: Create transformation plans. If you have worked with FFTW3 before, you already know all transformations are executed through plans. Here we need 3 plans: 2 plans for transforming signals

$𝕒$

and

$𝕓$

and one to reverse-transform the dot product of the transformed signals:

[code lang="c"]
fftw_plan pa = fftw_plan_dft_1d(2 * N - 1, signala_ext, outa, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan pb = fftw_plan_dft_1d(2 * N - 1, signalb_ext, outb, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_plan px = fftw_plan_dft_1d(2 * N - 1, out, out_shifted, FFTW_BACKWARD, FFTW_ESTIMATE);
[/code]

Note that, despite the fact you'll probably only have real-valued data, you should use complex transformations: using real versions of these will overcomplicate the whole process... believe me, I tried that! Note also that you may reuse the first plan if you are experienced enough with FFTW3.

Step 4: Transform, take dot products, and transform back. Use C's complex functions to get the conjugates of the values in the transformed signals:

[code lang="c"]
fftw_execute(pa);
fftw_execute(pb);

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

fftw_execute(px);
[/code]

Step 5: Reorder data. The final step is to reorder the resulting array so that the 0-frequency appears at the beginning, since the transformation will place this value in the middle of the signal (this is the equivalent to Matlab's fftwshift function):

[code lang="c"]
for (int i = 0; i < 2 * N - 1; i++)
result[i] = out_shifted[(i + N) % (2 * N - 1)] / (2 * N - 1);
[/code]

And you are done :)

The source code for building this library, if you want to incorporate it in your code under the terms of the GPL, is at my repository.