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])
