Why use the complex conjugate for Fourier spectra division?

Question

I saw several codes written where Fourier spectra are divided with the complex conjugate like this:

af = fftn(double(img1));
bf = fftn(double(img2));
cp = af .* conj(bf) ./ abs(af .* conj(bf));

in this script among others.

Is this related to handling complex division? Reading the documentation about the ./ operator, it is stated that it handles complex numbers. So is this wrong?:

af./bf

Answer 1

The expressions af./bf and af.*conj(bf)./abs(bf).^2 are completely equivalent in MATLAB, if that's what you're asking. There is no clear connection, however, between that question and the example you've shown. abs(bf).^2 does not appear in the denominator in your example.

The only reason conj() is being used in the code you've shown is because it is the Fourier dual of time inversion

I.e., f(t)<-->F(k) implies f(-t)<--->conj(F(k)), for real-valued time signals f(t).

This has a specific application to time delay analysis using phase correlation.

Answer 2

You could rewrite this expression avoiding the conjugation as

(af./bf)./abs(af./bf). 

However, the given form of the expression has the advantage that you can desingularize the division by adding a small epsilon to the denominator,

(af.*conj(bf))./(1e-40+abs(af.*conj(bf)))

Answer 3

Consider the following equivalent (within about 1e-15) code:

cpX = exp(1i*(angle(af)-angle(bf)));

You can compute the normalized cross power spectrum as you have shown with the complex conjugate (cp = af .* conj(bf) ./ abs(af .* conj(bf))) or by explicitly subtracting the phase as above.

Considering that the FFT of a shifted impulse is a complex exponential, the cpX equation should give some insight into how "phase correlation" allows you to find a translation between two images. The location of the peak in the inverse FFT gives the translation.