# How to get x-axis when using FFT to cross-correlate

Suppose there are two functions $$f(x)$$ and $$g(x)$$. Each function is nonzero only in the specified interval, $$[m_f, M_f]$$ and $$[m_g, M_g]$$. What I want to calculate is the cross-correlation between two functions: $$(f\star g)(t) = \int_{-\infty}^{\infty} f^\ast(x)g(x+t)dx$$ which can be obtained by applying convolution theorem: $$(f\star g)(t) = \mathcal{F}^{-1}\left\{{\mathcal{F}\left\{f\right\}}^\ast\cdot\mathcal{F}\left\{g\right\}\right\}=\mathcal{F}^{-1}\left\{\hat{f}^\ast\cdot\hat{g}\right\}$$ where $$\hat{.}$$ denotes forward transformation.

Now, to write this in FFT, I prepared four sets of arrays: x_f, x_g, y_f, y_g. x_f and x_g denote the $$x$$ domain of respective functions, and they are 1. zero-pad and 2. sampled in the same sampling. Basically, I took whichever was greater between $$M_f-m_f$$ and $$M_g-m_g$$, defined it to be $$L/2$$, made it as the interval axes. So, x_f[0] is $$m_f$$, and x_f[-1] is $$m_f+L\geq M_f$$. The same goes for x_g. y_f and y_g are defined to be $$f(x)$$ and $$g(x)$$ at each x_f and x_g.

If we take y_c = ifft(fft(y_f).conjugate()*fft(y_g)), we get a somewhat valid result. However, it is not easy to map y_c corresponds to what x_c so that y_c is equal to $$(f\star g)(x)$$ at x_c. How do you calculate this?

Also, what happens if there are three terms? i.e., $$(f\star (g\star h))(t)$$?

Let's say we have two discrete functions $$f[n]$$ and $$g[n]$$ that have finite support on $$[m_f, M_f]\ m_f,M_f \in \mathbb{Z}$$ and $$[m_g, M_g] \ m_g,M_g \in \mathbb{Z}$$. Than the lengths of the functions are $$N_f = M_f-m_f+1$$ and $$N_g = M_g-m_g+1$$

The length of the cross correlation $$h[n]$$ (or convolution) is simply

$$N_h = N_f+N_g-1$$

We treat the cross correlation as "convolve with the time flipped version of g[n]", i.e.

$$h[n] = f[n]*g[-n]$$

where $$*$$ is the convolution operator.

Since we are time flipping around $$n=0$$ The support of the time-flipped version is $$[-M_g,-m_g]$$, i.e. the limits change signs and flip positions. For example $$[3,8]$$ turns into $$[-8,-3]$$

The FFT implements the DFT (Discrete Fourier Transform). Multiplication in the DFT domain implements circular (not linear) convolution. Hence we need to zero pad both sequences to the length $$N_h$$ to avoid time domain aliasing. Note, that you need to zero-pad $$g[n]$$ AFTER time flipping it.

For convolution the support intervals just add, i.e. $$h[n]$$ will be non-zero on $$[m_h,M_h] = [m_f-M_g,M_f-m_g]$$

Programming languages like Matlab don't have a built-in way for managing the x-interval of an array so you have to keep track of it manually. Below is a code example that demonstrate the cross correlation between two different triangular waves.

%% create two tri-angular waves and cross correlate them
% wave 1: width 9, starting at n = 3;
f = conv(ones(5,1),ones(5,1)); f = f/max(f);
mf = [3,11];
% wave 2: width 13, starting at n = -5
g = conv(ones(7,1),ones(7,1)); g = g/max(g);
mg = [-5,7];

%% do the crosscorrelation
Nh = mf(2)+mg(2)-mf(1)-mg(1) + 1;       % length of cross correlation
h2 = ifft(fft(f,Nh).*fft(flip(g),Nh));  % F-domain convolution with padding
mh = mf + flip(-mg);                    % support interval of result

%% now plot it all
clf;
plot(mf(1):mf(2),f,'LineWidth',2);
hold on;
plot(mg(1):mg(2),g,'LineWidth',2);

plot(mh(1):mh(2),h2/max(h2),'LineWidth',2);
grid on
set(gca,'xlim',[-5 18]);
xlabel('Time in samples');
legend('f[n]','g[n]','h[n]');


Since the peak of $$g[n]$$ precedes the peak of $$f[n]$$ by 6 samples, the maximum of the cross correlation will be at $$n=6$$.