2D Discrete Fourier Transform implementation in MATLAB

2
votes

I am implementing the 2D Discrete Fourier Transform in Matlab using matrix multiplications.

I realize that this can be a separable operation, so I am creating a matrix for 1D DFT and multiplying it with the columns of an input image and then the rows of the image.

However, the values of the resulting 2D DFT have a large difference from the DFT that is calculated using the built-in function in MATLAB (i.e. fft2). Due to this, when performing the inverse DFT to recreate the image, the resultant image is not recreated correctly (i.e. it is not same as the original image, but it's the same if I use the fft2 function).

This is the code I wrote.

%% signal is a matrix of MxN size
function res=myDFT(signal)
    signal=double(signal);

    l=size(signal,1);   
    x=[1:l];

    %% u and x are matrices of MxM size
    [x u]=meshgrid(x,x);
    M1=l-1;

    pre_dft=exp(1i*(-2*pi)./M1)/sqrt(M1);

    pre_dft=(pre_dft.^(u.*x));

    %the below matrix will be multiplied with the rows of the signal
    post_dft=pre_dft;        

    % res is the resultant DFT of the signal matrix
    % 1D DFT matrix is first multiplied with columns and then the    
    % rows of signal matrix

    res=pre_dft*signal*post_dft;
end

I would really appreciate if someone could point out anything useful to edit to my code or point out a flaw in my theoretical understanding.

1
matlabimage-processingsignal-processingdft
You subtract the mean from the image, so you won't get the same results no matter what. – ThP
Comparing your result to the output of dftmtx from the Signal Processing Toolbox might help debug your code. – A. Donda
@ThP I have removed that line from the code – bumpk_sync
@A.Donda Yes I am trying that now – bumpk_sync

1 Answers

1
votes

OK, you've got a few errors that we need to fix... if you want this to work exactly like how fft2 works in MATLAB.

Error #1 - Wrong normalization factor

fft and ultimately fft2 have no normalization factors when computing the transform. You can get rid of that division statement when computing the first part of the DFT matrix. Also, the first part of the DFT matrix, in the exponent, you need to divide by the total length of the signal, not the length subtracted by 1.

Error #2 - Power operation for completing the last part of the DFT matrix is wrong

You calculate a meshgrid from 1 to l, but the power operations require that you start from 0. Therefore, you need to subtract u and x by 1 before you do the power calculation.

Error #3 - Application of FFT to rows and columns

As you have already noted, the 2D FFT is separable and can be done by performing the 1D FFT on the rows first, followed by the columns. The operation that you're doing unfortunately isn't correct. post_dft needs to be transposed as you want to apply the resulting intermediate operations over the columns.

As such, with all of the fixes I mentioned above, your corrected code is:

function res=myDFT(signal)
    signal=double(signal);

    l=size(signal,1);   
    x=[1:l];

    [x, u]=meshgrid(x,x);

    %// Error #1
    pre_dft=exp(1i*(-2*pi)./l); %// Change

    %// Error #2
    pre_dft=(pre_dft.^((u-1).*(x-1))); %// Change

    %// Error #3
    post_dft = pre_dft.'; %// Change

    res = pre_dft*signal*post_dft;
end

Testing the above with some random data, and comparing with fft2:

rng(123123);
in = rand(7);
out1 = myDFT(in);
out2 = fft2(in);

out1 contains the corrected custom implementation while out2 contains the result of MATLAB's fft2 algorithm. We get:

>> out1

out1 =

  Columns 1 through 4

  26.2182 + 0.0000i  -1.3805 + 1.0956i   2.2881 - 0.4435i  -0.8005 + 1.5133i
  -1.3067 + 0.3236i  -0.5703 - 1.3884i  -0.7127 + 0.1303i   3.2689 - 0.5995i
  -0.6136 - 2.0731i  -0.2776 + 1.2926i   0.1587 - 1.4504i  -0.5305 + 1.4018i
  -3.0703 + 0.3715i  -1.5094 - 0.0216i   0.2573 - 3.2934i  -1.1100 + 1.4180i
  -3.0703 - 0.3715i  -0.9940 + 0.1328i   1.7269 + 0.2915i   0.2814 + 2.2212i
  -0.6136 + 2.0731i   1.8351 + 1.0629i   1.2891 + 1.7418i   0.4402 + 1.8756i
  -1.3067 - 0.3236i  -0.4974 - 0.9678i  -2.2419 + 2.0839i  -1.6844 + 0.9781i

  Columns 5 through 7

  -0.8005 - 1.5133i   2.2881 + 0.4435i  -1.3805 - 1.0956i
  -1.6844 - 0.9781i  -2.2419 - 2.0839i  -0.4974 + 0.9678i
   0.4402 - 1.8756i   1.2891 - 1.7418i   1.8351 - 1.0629i
   0.2814 - 2.2212i   1.7269 - 0.2915i  -0.9940 - 0.1328i
  -1.1100 - 1.4180i   0.2573 + 3.2934i  -1.5094 + 0.0216i
  -0.5305 - 1.4018i   0.1587 + 1.4504i  -0.2776 - 1.2926i
   3.2689 + 0.5995i  -0.7127 - 0.1303i  -0.5703 + 1.3884i

>> out2

out2 =

  Columns 1 through 4

  26.2182 + 0.0000i  -1.3805 + 1.0956i   2.2881 - 0.4435i  -0.8005 + 1.5133i
  -1.3067 + 0.3236i  -0.5703 - 1.3884i  -0.7127 + 0.1303i   3.2689 - 0.5995i
  -0.6136 - 2.0731i  -0.2776 + 1.2926i   0.1587 - 1.4504i  -0.5305 + 1.4018i
  -3.0703 + 0.3715i  -1.5094 - 0.0216i   0.2573 - 3.2934i  -1.1100 + 1.4180i
  -3.0703 - 0.3715i  -0.9940 + 0.1328i   1.7269 + 0.2915i   0.2814 + 2.2212i
  -0.6136 + 2.0731i   1.8351 + 1.0629i   1.2891 + 1.7418i   0.4402 + 1.8756i
  -1.3067 - 0.3236i  -0.4974 - 0.9678i  -2.2419 + 2.0839i  -1.6844 + 0.9781i

  Columns 5 through 7

  -0.8005 - 1.5133i   2.2881 + 0.4435i  -1.3805 - 1.0956i
  -1.6844 - 0.9781i  -2.2419 - 2.0839i  -0.4974 + 0.9678i
   0.4402 - 1.8756i   1.2891 - 1.7418i   1.8351 - 1.0629i
   0.2814 - 2.2212i   1.7269 - 0.2915i  -0.9940 - 0.1328i
  -1.1100 - 1.4180i   0.2573 + 3.2934i  -1.5094 + 0.0216i
  -0.5305 - 1.4018i   0.1587 + 1.4504i  -0.2776 - 1.2926i
   3.2689 + 0.5995i  -0.7127 - 0.1303i  -0.5703 + 1.3884i

Looks fine to me!