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

3 visualizzazioni (ultimi 30 giorni)
Seyedalireza Abootorabi il 10 Mar 2023
Commentato: Paul il 13 Mar 2023
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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

Paul il 12 Mar 2023
Modificato: Paul il 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)
ans = 1×5
-1.2438 -0.7463 -0.2488 0.2488 0.7463
With numel(x) = numel(z) = 201, the frequency vectors should be
kxval = (-100:100)/201*dkx;
kzval = (-100:100)/201*dkz;
kxval(99:103)
ans = 1×5
-0.9950 -0.4975 0 0.4975 0.9950
% 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
Warning: Integer operands are required for colon operator when used as index.
Warning: Integer operands are required for colon operator when used as index.
That warning is troublesome, but it looks like rect comes out the way you want.
figure
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])
##### 1 CommentoMostra -1 commenti meno recentiNascondi -1 commenti meno recenti
Seyedalireza Abootorabi il 13 Mar 2023
thank you for the very clear answer.

Accedi per commentare.

### Più risposte (1)

Benjamin Thompson il 10 Mar 2023
Your problem is 201 by 201 samples in size? Making the size a power of 2. So use 256.
##### 3 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Benjamin Thompson il 13 Mar 2023
Yes, and the difference is probably due to fft2 padding your data out to 256 because it is required to operate on data sizes that are a power of 2, while your analytical method has not such restriction. The different data sizes causes the results of fft2 to be changed.
Paul il 13 Mar 2023
Based on the doc page fft2, it doesn't pad the input unless specified via the optional second and third input arguments.

Accedi per commentare.

### Categorie

Scopri di più su Transforms in Help Center e File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by