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 n2

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

Onlogn

).

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𝕒·𝕓¯

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.