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

Question

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)$?

Answer 1

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$.