## Cross Correlation Using FFTW3

**Posted:** June 23rd, 2010 **Tags:** C++, FFTW, Math

2 Comments »

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

${\mathcal{F}}^{-1}\left\{\mathcal{F}\left\{\mathbb{a}\right\}\xb7\overline{\mathcal{F}\left\{\mathbb{b}\right\}}\right\}$for input signals

$\mathbb{a}$and

$\mathbb{b}$).

**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

, 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

$\mathbb{b}$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.