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