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.
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.
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.
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.
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!
Collected from the Internet
Please contact [email protected] to delete if infringement.
Comments