How does subpixel image shifting using DFT really work?

I am trying to assess the quality of several image interpolation methods for an application that involves generating subpixel-shifted images. I thought I could compare the results of a subpixel shift using all these interpolation variants with some perfectly shifted image, but it's probably not possible to obtain it (what would be the need for interpolation then?).

I was thinking about using DFT + shifting in the frequency domain, and I am unsure as to how it actually works in comparison to explicitly interpolating the image (using bilinear, bicubic, etc...). I'm sure it can't possibly generate a perfectly shifted image, but I can't put my finger on it. Is subpixel shifting with DFT equivalent to applying interpolation and if so, which one? What is the bias of the pixel values in images obtained using this method? Thanks!

EDIT: After thinking the matter over, I figured since the FFT is an approximation (even more so the DFT) of the original function in terms of harmonics (sine functions), that it would amount to some sort of trigonometric interpolation. I recall a "Fourier series interpolation" formula for discrete data which was a trigonometric interpolation, but not sure if it's connected.

-
The fast Fourier transform (FFT) is an algorithm for the discrete Fourier transform. The DFT is not an approximation of the original function in terms of harmonics, but rather a projection of a signal onto a complex exponential orthogonal basis. – Bryan Feb 19 '13 at 19:57
Okay, but the signal itself is a sampled and quantized approximation of some intensity distribution, and the DFT is limited in respect to harmonic content in comparison to that theoretical distribution. You can get the exact signal back from IDFT, but there will be some bias if you do stuff (like shifting) to it before IDFTing it back. Or am I missing something? – neuviemeporte Feb 21 '13 at 6:54
The DFT indeed takes discretized inputs, but is not limited to quantizated inputs. What the signal is doesn't matter. As you've pointed out, you can get the exact signal back. However, I'm not sure what you mean by "shifting". The properties of shifting in the frequency domain are well-known (complex frequency translation in the time domain). If your desire is to shift in the "time" domain, than you need to think about the DFT dual of that. – Bryan Feb 21 '13 at 13:46
I mean that if I perform some operation on the DFT of a signal (like in my case - subpixel shift of an image in the "pixel domain" using the Fourier shift theorem), that the IDFT will return interpolated results as explained by @hotpaw2's answer. This interpolation is imperfect because the signal is not bandlimited and the DFT itself was calculated from a finite set of quantized (0-255) samples. – neuviemeporte Feb 21 '13 at 20:00

A DFT/FFT, plus added zero-padding in the frequency domain, then a longer IDFT/IFFT, returns interpolated points. These points will be interpolated using a periodic Sinc kernel, which is a perfect interpolation for original data that is strictly band-limited to below half the original sample rate. However, the data will be acted on as if it were circularly wrapped, which may produce odd results at the edges of some images. So you may want to pad the edges of the original source with a nice filler or framing color before interpolation.

If you upsample by 2X (zero-pad the FFT to double the length before the IFFT), then you can do a half-pixel shift using the interpolated points. 3X for a third pixel shift, etc. For shifting, you can throw away the original points plus any excess interpolated points to get the size desired.

-
@hotpaw2: the interpolating kernel for the DFT is not a sinc() of infinite extent, in fact the DFT is a discrete, finite transform. Interpolation by the DFT is equivalent to convolution with the Dirichlet kernel, also called periodic sinc() by some authors: en.wikipedia.org/wiki/Dirichlet_kernel – Arrigo Feb 19 '13 at 19:39
@Arrigo : Agree. Edited answer to fix. – hotpaw2 Feb 20 '13 at 9:01
@hotpaw2: when I pad the FFT to twice the size, the IFFT will yield a reconstruction of twice the size. Not sure what to do with the surplus? Thanks – neuviemeporte Feb 21 '13 at 8:35
Throw away the surplus points you don't need. In a 2X upsample, every other one is shifted, alternating with the reconstructed original points. In a 3X upsample, you get 2 shifted points (by 1/3 and by 2/3) alternating with the originals. Etc. The more you upsample, the more you throw away. – hotpaw2 Feb 21 '13 at 8:51

There are several key insights you need in order to understand how DFT allows you to shift an image.

First, Fourier's theorum: It's probably easier to look at the continuous (i.e., analog) case first. Imagine you have some function, call it g(t). For simplicity, let's say that g(t) is an analog audio recording, so it's a one-dimensional function, which is continuous, and represents the instantaneous pressure as a function of time.

Now, g(t) is one way we can represent our audio recording. Another is G(f). G(f) is the Fourier transform of g(t). So, G(f) == FT(g(t)). G(f) has all the same information as g(t), but it represents that information in the frequency domain instead of the time domain. There are some nit-picky details about Fourier Transforms, which I won't mention.

You can sort of think of G(f) as the "distribution of frequencies" contained in g(t). So, if g(t) is a sine-wave (i.e., a pure tone), then G(f) will be zero everywhere, except at the frequency of that tone. This is probably a good point to mention that G(f) is in general a complex function--that is to say that it returns complex numbers, which can be thought of having a real and imaginary component or a magnitude and phase.

One small digression here: Since g(t) is continuous (both in domain and range), G(f) is also continuous. So, how can G(f) be zero everywhere except for the tone frequency? Well, FT(sin(wt)) = $\delta(w)$. Where $\delta$ is the Dirac delta function.

Ok, so now we've got continuous FT's under our belt.

Here's the second insight: A Discrete Fourier Transform is to a Fourier Transform as a sampled signal is to an analog signal. In this case, the "discrete" refers to quantization of the function's domain (time or frequency), not it's range. (The sampled digital signal you get from your sound card is quantized in both domain and range.)

The digital byte-stream you get from your sound card contains "samples" of the original continuous (analog) signal from the microphone. If we take the DFT of our sampled g(t), we still get a G(f). G(f), remember, is just a different way of representing the information contained in g(t). If we obeyed Nyquist's theorum, the sampled signal g(t) contains all the "intelligence" of the original continuous signal, so our discrete G(f) must contain all the information from our original continuous signal. Parenthetically, G(f) is still a complex function.

This is where the magic of sub-pixel shifting comes in, but in this case I'm going to write about shifting the audio signal in time by less than a sample, since it's the same thing.

Remember how G(f) is a complex function? It needs to be complex in order to represent frequencies that aren't zero at t=0. Remember, sin(0)=0, so sin(2*0)=0, etc. But what if we started recording a quarter of the way through a cycle of the tone? That's where the phase part of G(f) comes from. In this case, the phase would be 90 degrees or pi/2 radians, depending on your preference for representing a quarter of a cycle. So G(tone_frequency) = 0 + i or $e^{{i \pi}\over{2}}$.

That means we can shift our audio recording in time (by any amount we choose, including a fraction of a sample time) simply by modifying the phase of G(t). Actually, that statement is perhaps a bit too casual. For an un-quantized, sampled signal the phase can be adjusted arbitrarily (this is part of the reason I made the distinction between quantization of domain and range earlier). However, for a quantized sampled signal (our byte-stream of audio, for example) the quantization step size (i.e., number of bits) determines the resolution with which we can adjust the phase. When we Inverse Fourier Transform G(f) (or DIFT it, for this sampled signal), the new set of samples g'(t) = DIFT(G(F)) will all be shifted in time by the amount we pick.

Applying this to your pixels simply means using a 2-dimensional FT instead of the 1-dimensional FT discussed here.

-