MATLAB reconstruction of FFT 1D with an image

1
votes

I'm trying to reconstruct the fft function in MATLAB. Referring to the manual the fft function does the following: If X is a matrix, then fft(X) treats the columns of X as vectors and returns the Fourier transform of each column.

So I construct after reading the image and converting it to double a temporary tab containing the current column to calculate. I then assemble the results in a new array (named FinalMatrix). (By the way the image is 256x256 pixels)

For the moment I get the following code but the output doesn't look near the output of the fft function in MATLAB.

 A = imread('start.jpg','jpg');
 Agray = rgb2gray(A);
 x=double(Agray);

 FourierMachine1 = fft(x);
 FourierMachine2 =fft2(x);
 figure;imshow(real(FourierMachine1));
 figure;imshow(real(FourierMachine2));

 [largeur,hauteur]= size(Agray);
 FinalMatrix= zeros(largeur,hauteur);

 X = zeros(largeur,1);
 for k = 1: largeur
    PremiereColonne= x(1:hauteur,k);
    for n = 1:hauteur
    X(k) = X(k) + PremiereColonne(n)*exp(-1i*pi/2*n*k); 
    end
    FinalMatrix(1:hauteur,k)=X;
 end
 S=fftshift(FinalMatrix);

 figure;imshow(real(S));

What am I exactly missing? Thanks for any advice.

enter image description here

1
matlabimage-processingsignal-processingfft
Could you post your sample image, please? – kkuilla
how do I post it here? I don't see a buton to edit the question. It's the famous girl with the head. (I cropped a 256x256 part out, the head.) – moutaindwarf
added it, stupid me. – moutaindwarf

1 Answers

2
votes

The main problem is that you're not calculating the 1D DFT for each column properly. You also aren't iterating over each column correctly. Therefore, you need to have an additional loop to iterate over each column, then you'll need a pair of nested for loops. One loop will control the access of where to write in the output vector and another loop iterates over each sample of an individual column.

Also, the exponential term for the DFT is incorrect. If you recall, the equation of the DFT is:


(source: alwayslearn.com)

That isn't the exponential that I see in your code. In addition, because MATLAB starts indexing at 1, you must subtract the k and n in the exponential by 1 so that we start from 0, not at 1, as per the definition of the DFT.

With these changes, we get:

for k = 1: hauteur %// Change
    X = zeros(largeur,1);
    PremiereColonne = x(1:hauteur,k);
    for y = 1 : largeur %// Change
        for n = 1 : largeur %// Change
            X(y) = X(y) + PremiereColonne(n)*exp(-1i*2*pi*(n-1)*(y-1)/largeur); %// Change
        end
    end
    FinalMatrix(1:hauteur,k)=X;
 end

Here's proof to show you that it works. Let's generate a random 10 x 10 matrix with a random seed set to 123:

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

x =

    0.6965    0.3432    0.6344    0.0921    0.6240    0.1206    0.6693    0.0957    0.3188    0.7050
    0.2861    0.7290    0.8494    0.4337    0.1156    0.8263    0.5859    0.8853    0.6920    0.9954
    0.2269    0.4386    0.7245    0.4309    0.3173    0.6031    0.6249    0.6272    0.5544    0.3559
    0.5513    0.0597    0.6110    0.4937    0.4148    0.5451    0.6747    0.7234    0.3890    0.7625
    0.7195    0.3980    0.7224    0.4258    0.8663    0.3428    0.8423    0.0161    0.9251    0.5932
    0.4231    0.7380    0.3230    0.3123    0.2505    0.3041    0.0832    0.5944    0.8417    0.6917
    0.9808    0.1825    0.3618    0.4264    0.4830    0.4170    0.7637    0.5568    0.3574    0.1511
    0.6848    0.1755    0.2283    0.8934    0.9856    0.6813    0.2437    0.1590    0.0436    0.3989
    0.4809    0.5316    0.2937    0.9442    0.5195    0.8755    0.1942    0.1531    0.3048    0.2409
    0.3921    0.5318    0.6310    0.5018    0.6129    0.5104    0.5725    0.6955    0.3982    0.3435

Using this, I get this for your resulting matrix:

>> FinalMatrix

FinalMatrix =

  Columns 1 through 5

   5.4420 + 0.0000i   4.1278 + 0.0000i   5.3795 + 0.0000i   4.9542 + 0.0000i   5.1894 + 0.0000i
  -0.7167 + 0.5845i   0.3827 - 0.0441i   0.6872 - 1.1141i  -0.1564 + 0.9087i  -0.3029 + 0.8021i
   0.2819 - 0.0768i   0.6751 + 0.0040i   0.2472 + 0.1070i  -1.2778 + 0.1311i  -0.2934 + 0.6208i
   1.0166 + 0.1215i  -1.1997 - 0.5153i   0.0443 - 0.0726i  -0.2362 - 0.4714i   1.0213 - 0.3459i
  -0.2040 - 0.2060i  -0.0361 + 0.0325i  -0.5435 + 0.1292i  -0.1884 - 0.0683i  -0.1153 + 0.8681i
   0.7670 + 0.0000i  -0.3402 - 0.0000i   0.0941 - 0.0000i  -0.3156 - 0.0000i   0.4307 - 0.0000i
  -0.2040 + 0.2060i  -0.0361 - 0.0325i  -0.5435 - 0.1292i  -0.1884 + 0.0683i  -0.1153 - 0.8681i
   1.0166 - 0.1215i  -1.1997 + 0.5153i   0.0443 + 0.0726i  -0.2362 + 0.4714i   1.0213 + 0.3459i
   0.2819 + 0.0768i   0.6751 - 0.0040i   0.2472 - 0.1070i  -1.2778 - 0.1311i  -0.2934 - 0.6208i
  -0.7167 - 0.5845i   0.3827 + 0.0441i   0.6872 + 1.1141i  -0.1564 - 0.9087i  -0.3029 - 0.8021i

  Columns 6 through 10

   5.2262 + 0.0000i   5.2544 + 0.0000i   4.5066 + 0.0000i   4.8248 + 0.0000i   5.2380 + 0.0000i
   0.3612 + 0.2466i   0.1933 - 0.8737i   0.2852 - 0.7816i  -0.5467 - 1.0722i   0.3197 - 1.0983i
  -1.1157 - 0.2910i   0.2011 + 0.0622i   0.0105 - 0.6416i   0.8486 + 0.3168i   0.6180 - 0.0535i
  -0.5658 - 0.4700i   0.8047 + 0.4189i  -0.7276 + 0.9442i  -0.8086 - 0.4696i   0.2864 - 0.7590i
  -0.4355 - 0.3588i  -0.9470 + 0.0380i  -0.5385 - 0.5152i  -0.3600 + 0.0700i   0.2547 - 0.3598i
  -0.5083 - 0.0000i   0.9345 + 0.0000i  -1.6087 - 0.0000i   0.0961 - 0.0000i  -1.1459 - 0.0000i
  -0.4355 + 0.3588i  -0.9470 - 0.0380i  -0.5385 + 0.5152i  -0.3600 - 0.0700i   0.2547 + 0.3598i
  -0.5658 + 0.4700i   0.8047 - 0.4189i  -0.7276 - 0.9442i  -0.8086 + 0.4696i   0.2864 + 0.7590i
  -1.1157 + 0.2910i   0.2011 - 0.0622i   0.0105 + 0.6416i   0.8486 - 0.3168i   0.6180 + 0.0535i
   0.3612 - 0.2466i   0.1933 + 0.8737i   0.2852 + 0.7816i  -0.5467 + 1.0722i   0.3197 + 1.0983i

Verifying with MATLAB's fft function, where each column is processed individually, I get:

>> out = fft(x)

out =

  Columns 1 through 5

   5.4420 + 0.0000i   4.1278 + 0.0000i   5.3795 + 0.0000i   4.9542 + 0.0000i   5.1894 + 0.0000i
  -0.7167 + 0.5845i   0.3827 - 0.0441i   0.6872 - 1.1141i  -0.1564 + 0.9087i  -0.3029 + 0.8021i
   0.2819 - 0.0768i   0.6751 + 0.0040i   0.2472 + 0.1070i  -1.2778 + 0.1311i  -0.2934 + 0.6208i
   1.0166 + 0.1215i  -1.1997 - 0.5153i   0.0443 - 0.0726i  -0.2362 - 0.4714i   1.0213 - 0.3459i
  -0.2040 - 0.2060i  -0.0361 + 0.0325i  -0.5435 + 0.1292i  -0.1884 - 0.0683i  -0.1153 + 0.8681i
   0.7670 + 0.0000i  -0.3402 + 0.0000i   0.0941 + 0.0000i  -0.3156 + 0.0000i   0.4307 + 0.0000i
  -0.2040 + 0.2060i  -0.0361 - 0.0325i  -0.5435 - 0.1292i  -0.1884 + 0.0683i  -0.1153 - 0.8681i
   1.0166 - 0.1215i  -1.1997 + 0.5153i   0.0443 + 0.0726i  -0.2362 + 0.4714i   1.0213 + 0.3459i
   0.2819 + 0.0768i   0.6751 - 0.0040i   0.2472 - 0.1070i  -1.2778 - 0.1311i  -0.2934 - 0.6208i
  -0.7167 - 0.5845i   0.3827 + 0.0441i   0.6872 + 1.1141i  -0.1564 - 0.9087i  -0.3029 - 0.8021i

  Columns 6 through 10

   5.2262 + 0.0000i   5.2544 + 0.0000i   4.5066 + 0.0000i   4.8248 + 0.0000i   5.2380 + 0.0000i
   0.3612 + 0.2466i   0.1933 - 0.8737i   0.2852 - 0.7816i  -0.5467 - 1.0722i   0.3197 - 1.0983i
  -1.1157 - 0.2910i   0.2011 + 0.0622i   0.0105 - 0.6416i   0.8486 + 0.3168i   0.6180 - 0.0535i
  -0.5658 - 0.4700i   0.8047 + 0.4189i  -0.7276 + 0.9442i  -0.8086 - 0.4696i   0.2864 - 0.7590i
  -0.4355 - 0.3588i  -0.9470 + 0.0380i  -0.5385 - 0.5152i  -0.3600 + 0.0700i   0.2547 - 0.3598i
  -0.5083 + 0.0000i   0.9345 + 0.0000i  -1.6087 + 0.0000i   0.0961 + 0.0000i  -1.1459 + 0.0000i
  -0.4355 + 0.3588i  -0.9470 - 0.0380i  -0.5385 + 0.5152i  -0.3600 - 0.0700i   0.2547 + 0.3598i
  -0.5658 + 0.4700i   0.8047 - 0.4189i  -0.7276 - 0.9442i  -0.8086 + 0.4696i   0.2864 + 0.7590i
  -1.1157 + 0.2910i   0.2011 - 0.0622i   0.0105 + 0.6416i   0.8486 - 0.3168i   0.6180 + 0.0535i
   0.3612 - 0.2466i   0.1933 + 0.8737i   0.2852 + 0.7816i  -0.5467 + 1.0722i   0.3197 + 1.0983i

We can see that these do match. If you absolutely want to be sure, let's check every element between the two arrays and show that their difference in magnitude is less than some threshold... say... 1e-10:

>> all(abs(FinalMatrix(:) - out(:)) < 1e-10)

ans =

     1