Search code examples
matlab2dfftdft

Bug in 2D Fourier Transform implementation

I am trying to implement a 2D DFT using a combination of 1D DFT in Matlab. When comparing my results against Matlab's inbuilt function (fft2) I realized that I have the following issues:

  1. The sign of the imaginary part is being inverted. i.e + changed to - and vice versa.
  2. The rows are being ordered in descending order except for the first one.

This image shows the comparison between the two results. The red numbers on the side are there to indicate the reordering issue.

My code is as follows:

x = imread('img.png');
x = double(x);
x = x(1:12,1:5)

% FFT
Xw = complex(zeros(size(x)));
for row = 1:size(x,1)
    Xw(row,:) = fft(x(row,:));
end

for col = 1:size(x,2)
    Xw(:,col) = fft(Xw(:,col)');
end

Could someone kindly indicate where my problem is please? Thanks

Solution

  • The ' operator is for the complex transpose, which means that the matrix is transposed and the conjugate of the values is taken. This means that the signs of the complex values are reversed. You don't notice this with real numbers because technically, the imaginary component is zero.

    You want to use .' instead to preserve the sign of the imaginary component as performing the intermediate FFT per row will result in complex values. Using just ' will change the imaginary components of the intermediate results, which will give you the wrong results:

    for col = 1:size(x,2)
    Xw(:,col) = fft(Xw(:,col).');
end

    BTW, as a minor note, there is no need to transpose the intermediate columns of your result at all. If you provide a single vector, fft will operate along the first non-singleton dimension, and so using a transpose is superfluous. It would actually help if you didn't transpose your result at all:

    for col = 1:size(x,2)
    Xw(:,col) = fft(Xw(:,col));
end

    Here's an example. I've generated a random 10 x 10 matrix:

    rng(123);
x = rand(10);

    Using your code with correcting the transpose (and without the indexing at the beginning), we get:

    >> Xw

Xw =

  Columns 1 through 5

  50.1429 + 0.0000i  -0.4266 - 0.2624i   0.8803 + 0.9311i  -0.0526 + 1.7067i   0.7187 + 0.7161i
   0.5066 - 2.4421i   2.7551 - 1.7421i  -1.9994 + 0.6052i   0.2891 + 3.4182i   0.5300 + 2.4417i
   0.1956 + 0.1790i   4.1461 + 1.9648i  -1.3781 - 1.0303i  -0.6872 - 1.0103i  -1.2184 - 0.5783i
  -0.3645 - 1.6193i  -1.8470 - 1.3445i   4.1555 + 0.7432i   2.3707 + 3.8265i  -1.9526 + 1.9464i
  -3.1136 - 0.3704i   2.4132 - 1.0795i  -0.2255 + 1.3062i   0.8436 - 0.5157i  -0.3493 - 0.9994i
  -1.5962 + 0.0000i  -0.3780 - 1.4055i   1.6242 - 0.4842i   0.4457 + 0.4718i  -0.1794 - 2.0014i
  -3.1136 + 0.3704i   0.0134 + 0.1267i   1.0630 + 1.4563i  -0.8864 - 0.3174i  -0.5720 + 1.3041i
  -0.3645 + 1.6193i   0.7028 + 0.2797i   0.1064 + 2.0705i   2.1644 + 0.1685i   0.3095 + 0.7426i
   0.1956 - 0.1790i   2.3511 + 2.1440i   0.7301 - 0.8264i  -1.1974 - 0.3794i  -2.4981 + 1.2363i
   0.5066 + 2.4421i  -3.5897 + 0.7444i   1.2191 - 3.6386i  -2.9659 - 1.6626i  -2.0339 + 0.0880i

  Columns 6 through 10

   2.0373 + 0.0000i   0.7187 - 0.7161i  -0.0526 - 1.7067i   0.8803 - 0.9311i  -0.4266 + 0.2624i
  -1.8782 - 0.9047i  -2.0339 - 0.0880i  -2.9659 + 1.6626i   1.2191 + 3.6386i  -3.5897 - 0.7444i
   2.3752 + 1.8811i  -2.4981 - 1.2363i  -1.1974 + 0.3794i   0.7301 + 0.8264i   2.3511 - 2.1440i
   4.5213 + 0.9237i   0.3095 - 0.7426i   2.1644 - 0.1685i   0.1064 - 2.0705i   0.7028 - 0.2797i
  -1.2259 + 2.1690i  -0.5720 - 1.3041i  -0.8864 + 0.3174i   1.0630 - 1.4563i   0.0134 - 0.1267i
   6.2411 + 0.0000i  -0.1794 + 2.0014i   0.4457 - 0.4718i   1.6242 + 0.4842i  -0.3780 + 1.4055i
  -1.2259 - 2.1690i  -0.3493 + 0.9994i   0.8436 + 0.5157i  -0.2255 - 1.3062i   2.4132 + 1.0795i
   4.5213 - 0.9237i  -1.9526 - 1.9464i   2.3707 - 3.8265i   4.1555 - 0.7432i  -1.8470 + 1.3445i
   2.3752 - 1.8811i  -1.2184 + 0.5783i  -0.6872 + 1.0103i  -1.3781 + 1.0303i   4.1461 - 1.9648i
  -1.8782 + 0.9047i   0.5300 - 2.4417i   0.2891 - 3.4182i  -1.9994 - 0.6052i   2.7551 + 1.7421i

    I've also verified that this works without transposing. You'll get the same thing as if you were to perform .'. Now, using fft2, we get:

    >> Xw2 = fft2(x)

Xw2 =

  Columns 1 through 5

  50.1429 + 0.0000i  -0.4266 - 0.2624i   0.8803 + 0.9311i  -0.0526 + 1.7067i   0.7187 + 0.7161i
   0.5066 - 2.4421i   2.7551 - 1.7421i  -1.9994 + 0.6052i   0.2891 + 3.4182i   0.5300 + 2.4417i
   0.1956 + 0.1790i   4.1461 + 1.9648i  -1.3781 - 1.0303i  -0.6872 - 1.0103i  -1.2184 - 0.5783i
  -0.3645 - 1.6193i  -1.8470 - 1.3445i   4.1555 + 0.7432i   2.3707 + 3.8265i  -1.9526 + 1.9464i
  -3.1136 - 0.3704i   2.4132 - 1.0795i  -0.2255 + 1.3062i   0.8436 - 0.5157i  -0.3493 - 0.9994i
  -1.5962 + 0.0000i  -0.3780 - 1.4055i   1.6242 - 0.4842i   0.4457 + 0.4718i  -0.1794 - 2.0014i
  -3.1136 + 0.3704i   0.0134 + 0.1267i   1.0630 + 1.4563i  -0.8864 - 0.3174i  -0.5720 + 1.3041i
  -0.3645 + 1.6193i   0.7028 + 0.2797i   0.1064 + 2.0705i   2.1644 + 0.1685i   0.3095 + 0.7426i
   0.1956 - 0.1790i   2.3511 + 2.1440i   0.7301 - 0.8264i  -1.1974 - 0.3794i  -2.4981 + 1.2363i
   0.5066 + 2.4421i  -3.5897 + 0.7444i   1.2191 - 3.6386i  -2.9659 - 1.6626i  -2.0339 + 0.0880i

  Columns 6 through 10

   2.0373 + 0.0000i   0.7187 - 0.7161i  -0.0526 - 1.7067i   0.8803 - 0.9311i  -0.4266 + 0.2624i
  -1.8782 - 0.9047i  -2.0339 - 0.0880i  -2.9659 + 1.6626i   1.2191 + 3.6386i  -3.5897 - 0.7444i
   2.3752 + 1.8811i  -2.4981 - 1.2363i  -1.1974 + 0.3794i   0.7301 + 0.8264i   2.3511 - 2.1440i
   4.5213 + 0.9237i   0.3095 - 0.7426i   2.1644 - 0.1685i   0.1064 - 2.0705i   0.7028 - 0.2797i
  -1.2259 + 2.1690i  -0.5720 - 1.3041i  -0.8864 + 0.3174i   1.0630 - 1.4563i   0.0134 - 0.1267i
   6.2411 + 0.0000i  -0.1794 + 2.0014i   0.4457 - 0.4718i   1.6242 + 0.4842i  -0.3780 + 1.4055i
  -1.2259 - 2.1690i  -0.3493 + 0.9994i   0.8436 + 0.5157i  -0.2255 - 1.3062i   2.4132 + 1.0795i
   4.5213 - 0.9237i  -1.9526 - 1.9464i   2.3707 - 3.8265i   4.1555 - 0.7432i  -1.8470 + 1.3445i
   2.3752 - 1.8811i  -1.2184 + 0.5783i  -0.6872 + 1.0103i  -1.3781 + 1.0303i   4.1461 - 1.9648i
  -1.8782 + 0.9047i   0.5300 - 2.4417i   0.2891 - 3.4182i  -1.9994 - 0.6052i   2.7551 + 1.7421i

    We can show these two are equal by ensuring that all of the element-wise differences between the two matrices is less than a small threshold... say 1e-10:

    >> all(abs(Xw(:)-Xw2(:)) <= 1e-10) 

ans =

     1

    Take away message

    Never use ' for the transpose.... ever. Luis Mendo is a strong advocate of this very fact.