Why am I not getting the intended coefficients in this 2D Fourier demonstration?

Question

I am trying to demonstrate how the 2D Fourier decomposition of an image works with Matlab and a very, very simple example.

I create a 4 x 4 pixel image with cosine basis vectors as such:

%Basis vectors with coefficients:

f = @(x, y,a0,a1,a2,a3,b0,b1,b2,b3,c0,c1,c2,c3,d0,d1,d2,d3)... 
a0*cos(0*x + 0*y) + a1*cos(0*x + 1*y) +...
a2*cos(0*x + 2*y) + a3*cos(0*x + 3*y) + ...
b0*cos(1*x + 0*y) + b1*cos(1*x + 1*y) + ...
b2*cos(1*x + 2*y) + b3*cos(1*x + 3*y) + ...
c0*cos(2*x + 1*y) + c1*cos(2*x + 1*y) + ...
c2*cos(2*x + 2*y) + c3*cos(2*x + 3*y) + ...
d0*cos(3*x + 1*y) + d1*cos(3*x + 1*y) + ...
d2*cos(3*x + 2*y) + d3*cos(3*x + 3*y)

%Matrix of image space:
[X,Y] = meshgrid(0:3);

%Creating the image from Fourier space by assigning coefficients:

F = f(X,Y,1,1,-1,1,1,1,1,6,3,1,1,1,7,1,1,1);
Fourier = fft2(F);

F =
 27 -5.53409 3.93452 -10.2602
 -0.177994 -12.1365 8.43196 1.46219
 3.50819 6.86722 2.91194 -3.6156
 -18.1397 2.3525 -6.23362 5.35198

Fourier =

 (5.72275,0) (3.14567,1.38924) (36.7478,0) (3.14567,-1.38924)
 (5.46847,-14.2485) (33.0684,2.46052) (43.5603,-51.0061) (11.87,-9.05285)
 (43.9012,0) (44.1778,-31.8071) (63.0469,0) (44.1778,31.8071)
 (5.46847,14.2485) (11.87,9.05285) (43.5603,51.0061) (33.0684,-2.46052)

The image is:

Now I naively want to get back the coefficients with Fourier = fft2(F), but it is not happening - to begin with, there are imaginary components, which I presume are part of a phase shift I did not introduce (I only used cosines in the creation of the original image).

What are my mistakes and misconceptions undermining this demonstration?

I see that the basis vectors can actually be expressed as complex numbers, e.g.

$$\cos(3x + 1x) = \frac 1 2 e^{-3 i x - i y} + \frac 1 2 e^{3 i x + i y}.$$

Yet, predictably, expressing the function with complex terms does not change much - below is the function expressed with complex terms for access:

f = @(x, y,a1,a2,a3,a4,b1,b2,b3,b4,c1,c2,c3,c4,d1,d2,d3,d4)... 
a1*(1/2*exp(-0*i*x-0*i*y) + 1/2*exp(0*i*x+0*i*y)) + a2*(1/2*exp(-0*i*x-1*i*y) + 1/2*exp(0*i*x+1*i*y)) +...
a3*(1/2*exp(-0*i*x-2*i*y) + 1/2*exp(0*i*x+2*i*y)) + a4*(1/2*exp(-0*i*x-3*i*y) + 1/2*exp(0*i*x+3*i*y)) + ...
b1*(1/2*exp(-1*i*x-0*i*y) + 1/2*exp(1*i*x+0*i*y)) + b2*(1/2*exp(-1*i*x-1*i*y) + 1/2*exp(1*i*x+1*i*y)) + ...
b3*(1/2*exp(-1*i*x-2*i*y) + 1/2*exp(1*i*x+2*i*y)) + b4*(1/2*exp(-1*i*x-3*i*y) + 1/2*exp(1*i*x+3*i*y)) + ...
c1*(1/2*exp(-2*i*x-0*i*y) + 1/2*exp(2*i*x+0*i*y)) + c2*(1/2*exp(-2*i*x-1*i*y) + 1/2*exp(2*i*x+1*i*y)) + ...
c3*(1/2*exp(-2*i*x-2*i*y) + 1/2*exp(2*i*x+2*i*y)) + c4*(1/2*exp(-2*i*x-3*i*y) + 1/2*exp(2*i*x+3*i*y)) + ...
d1*(1/2*exp(-3*i*x-0*i*y) + 1/2*exp(3*i*x+0*i*y)) + d2*(1/2*exp(-3*i*x-1*i*y) + 1/2*exp(3*i*x+1*i*y)) + ...
d3*(1/2*exp(-3*i*x-2*i*y) + 1/2*exp(3*i*x+2*i*y)) + d4*(1/2*exp(-3*i*x-3*i*y) + 1/2*exp(3*i*x+3*i*y))

---

If I run imshow(ifft2(Fourier)) I do get the original image with the initiall 4 x 4 setup.

So I assume I am on the right track, but I can't match the coefficients in the 2D DFT.

I have been able to do it for the DC component:

g = @(x, y, a1)...
dot(F(:), (a0* cos(0*x + 0*y))(:))
g(X,Y,1)
% ans = 5.72275

Which performs the Frobenius product of the image matrix multiplied by the matrix of the DC component:

$$a_0 \cos(x + y).$$

The result is 5.72275.

But I haven't been able to do the same with, say, the coefficient for

$$b_3\cos(x +3y)$$

with

q = @(x, y, b3)...
dot(F(:), (b3* cos(1*x + 3*y))(:))
q(X,Y,6)
% ans = 449.363

which yields 449.363, while there are no coefficients greater than 80 in the Fourier matrix Fourier. Evidently I am messing up treating the complex values correctly...

Answer 1

Your misconception is to assume that you can compute Fourier coefficients with the fft. You cannot (really). Fourier coefficients are found by integrating a continuous time (periodic) function after multiplying it with a harmonic kernel. The cosine function gives two real-valued Fourier coefficients equal to half its amplitude in this case.

You are considering a matrix though which means that you have sampled your cosine function at a discrete grid (which is not a problem in your case) and truncated it to a finite width (which is exactly the problem you have). Then you are computing the FFT which is a fast way of computing a DFT (Discrete Fourier Transform). The DFT (of discrete vectors) is not the same as a Fourier transform (which is defined on continuous functions).

That said, your DFT can approximate Fourier coefficients, as the Fourier integral can be approximated by a sum, which is precisely the definition of the DFT then. The approximation errors that happen in this case are sampling and truncation.

For functions that are strictly bandlimited (like your cosine), sampling does not cause any error as long as you satisfy the Nyquist criterion. Conversely, for functions that are strictly time-limited, truncation does not cause any error as long as the whole function fits in your window. Unfortunately, it can be shown that no function can be time-limited and band-limited at the same time. For this reason, the DFT can never accurately predict the result of the Fourier transform.

There is hope though, since you're not after a Fourier transform but instead you're after Fourier coefficients. These are defined on periodic continuous-time functions. For these, the DFT can give you accurate predictions of the Fourier coefficients, as long as you make sure that the samples you are considering contain a whole period of the signal. This needs to be accurate, i.e., if you periodify your extracted window it should look like the original function. For instance, if the period of your function is 1 and you want to extract four samples, they need to be at 0, 0.25, 0.5, 0.75 (i.e., not 0, 1/3, 2/3, 1, which is a very common mistake in this setting).

Your Fourier coefficients will then be real-valued if and only if your function is symmetric, which in the sampled case means x[n] = x[N-n] for n=1, 2, ..., N-1. If you sample your cosines correctly, this should work.