# fft2 results don't match the analytical Fourier transform results

3 views (last 30 days)

Show older comments

Hello everyone,

I am trying to implement fft2 and compare the results with the analytical Fourier trasnform. For that I have chosen a constant rectangular field with constant value (the constant c = 1) which its analytical Fourier trasnform already exist. However the results of my fft2 and analytival Fourier trasnform do not match each other. Can you help me to find out my mistake?

clear

clc

% Define the dimensions of the square field

Lx = 2; % length of the rectangle

Lz = 1; % width of the rectangle

x = 0:0.005:1; % x vector

z = 0:0.005:1; % z vector

xgrd = length(x);

zgrd = length(z);

kxgrd = length(x);

kzgrd = length(z);

dkx = 1/(x(2)-x(1)); % Sampling frequency in x

dkz = 1/(z(2)-z(1)); % smapling frequency in z

kxval = (-length(x)/2:length(x)/2-1)*(dkx/length(x)); % Frequency vector

kzval = (-length(z)/2:length(z)/2-1)*(dkz/length(z));

% Define the size of the rectangular field and the constant value

c = 1; % constant value

rect = zeros(xgrd,zgrd); % initialize the rectangular field to all zeros

rect(xgrd/4+1:3*xgrd/4,zgrd/4+1:3*zgrd/4) = c; % set the central region to the constant value

% Evaluate the analytical Fourier transform on a grid of u and v values

for indkx = 1:kxgrd

for indkz = 1:kzgrd

F_analytical(indkx,indkz) = Lx*Lz*sinc(pi*Lx*kxval(indkx))*sinc(pi*Lz*kzval(indkz)); % must match the fft2 results

end

end

func = rect; % cosine signal

% Compute the Fourier transform

func_wo_shift = fft2(func);

func_shift = fftshift(fft2(func));

##### 0 Comments

### Accepted Answer

Paul
on 12 Mar 2023

Edited: Paul
on 12 Mar 2023

Hi Seyedalireza,

The code below get closer, maybe all the way there, by addressing a few issues

% Define the dimensions of the square field

Lx = 2; % length of the rectangle

Lz = 1; % width of the rectangle

I multiplied by Lx and Lz to get the correct physical dimensions

x = Lx*(0:0.005:1); % x vector

z = Lz*(0:0.005:1); % z vector

Define dx and dz for use below

dx = (x(2)-x(1));

dz = (z(2)-z(1));

xgrd = length(x);

zgrd = length(z);

kxgrd = length(x);

kzgrd = length(z);

Defined in terms of dx and dz

dkx = 1/dx; % Sampling frequency in x

dkz = 1/dz; % smapling frequency in z

The next two lines are incorrect, I believe.

kxval = (-length(x)/2:length(x)/2-1)*(dkx/length(x)); % Frequency vector

kzval = (-length(z)/2:length(z)/2-1)*(dkz/length(z));

The middle of the frequency vector should be zero, which is important in this problem because that's where all the action occurs. But 0 doesn't show up.

kxval(99:103)

With numel(x) = numel(z) = 201, the frequency vectors should be

kxval = (-100:100)/201*dkx;

kzval = (-100:100)/201*dkz;

kxval(99:103)

% Define the size of the rectangular field and the constant value

c = 1; % constant value

rect = zeros(xgrd,zgrd); % initialize the rectangular field to all zeros

rect(xgrd/4+1:3*xgrd/4,zgrd/4+1:3*zgrd/4) = c; % set the central region to the constant value

That warning is troublesome, but it looks like rect comes out the way you want.

figure

pcolor(x,z,rect.'),shading interp

xlabel('x'),ylabel('Z')

I think this is the correct expression for F_analytical, assuming we're only interested in its amplitude and not its phase. Note that sinc is defined as sinc(x) = sin(pi*x)/(pi*x);

% Evaluate the analytical Fourier transform on a grid of u and v values

for indkx = 1:kxgrd

for indkz = 1:kzgrd

% F_analytical(indkx,indkz) = Lx*Lz*sinc(pi*Lx*kxval(indkx))*sinc(pi*Lz*kzval(indkz)); % must match the fft2 results

F_analytical(indkx,indkz) = Lx/2*Lz/2*sinc(Lx/2*kxval(indkx))*sinc(Lz/2*kzval(indkz)); % must match the fft2 results

end

end

func = rect; % cosine signal

% Compute the Fourier transform

func_wo_shift = fft2(func);

func_shift = fftshift(fft2(func));

figure

pcolor(kxval,kzval,abs(F_analytical.')),colorbar

xlabel('Kx'); ylabel('Kz')

axis([-10 10 -10 10])

I'm going to assume that F_analytical is supposed to be the Continous Time Fourier Transform. If so, the result from fft2 has to be multiplied by dx*dz for comparison

figure

pcolor(kxval,kzval,abs(dx*dz*func_shift.')),colorbar

xlabel('Kx');ylabel('Kz')

axis([-10 10 -10 10])

### More Answers (1)

Benjamin Thompson
on 10 Mar 2023

Your problem is 201 by 201 samples in size? Making the size a power of 2. So use 256.

##### 3 Comments

Benjamin Thompson
on 13 Mar 2023

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!