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:

enter image description here

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))(:))
% 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)$$


q = @(x, y, b3)...
dot(F(:), (b3* cos(1*x + 3*y))(:))
% 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...

  • $\begingroup$ what i would do, if you're using the 2-dim FFT (instead of some flavor of the DCT), i would double the image size in both height and width by reflecting the image over both the vertical and horizontal axes. that way, when the FFT periodically extends any horizontal or vertical scan line, that scan line has no discontinuity at the edges of the image. (i realize there are discontinuities in the image, but that's what you're looking at, dude.) in fact, this reflection operation is what you have to do to get only cosines (which are even symmetry). $\endgroup$ – robert bristow-johnson Apr 25 at 21:02

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.

  • $\begingroup$ Any leads at a practical implementation level, using Matlab? $\endgroup$ – Antoni Parellada Apr 27 at 22:15
  • $\begingroup$ Depends on what you want to show exactly. So far you only said "demonstrate how the 2D Fourier decomposition of an image works". To do that, I'd take an empty matrix M, put a single 1 somewhere, and then show real(fft2(M)). Make it 256x256 to look nice. Place it top-left: result is constant (DC->const). Place it (i,1) for i>1: harmonic along rows. Place it (1,j) harmonic along columns. Place it (i,j) 2-D harmonic. Note that the FFT is linear and every image can be decomposed like that. So the harmonics you are showing represent the basis for all transformed images... $\endgroup$ – Florian Apr 29 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.