Unexpected phases in fresnel diffraction using fft2

Hello everyone,

I'm relative new to Matlab and encountered trouble trying to realize a Fresnel-diffraction using the Fast-Fourier Algorithm. Please apologize my somehow lacking English skills. Since it's not my mother tongue, some of my explanations may sound a little bit weird.

I want to propagate a N x N matrix, whereas N is a power of 2. The matrix I use as a test case describes a square which is symmetric to the middel of the matrix, I generate the matrix with the following code snippet:

      A = zeros(256,256);
      A(119:138, 119:138) = 100;
      fresnelpropagation(A, 100, 0.682, 0.39, 0.39); %function call

I get a small 20px x 20px square in the middel of a 256 x 256 matrix. The function I call looks like this:

function Aout=fresnelpropagation2( Ain ,z , lambda, dx , dy )
[nx, ny] = size(Ain);

k = 2*pi / lambda;

%dx and dy in distance z
dx2 = z*lambda/(nx  * dx);
dy2 = z*lambda/(ny  * dy);

%"coordinates" of matrix Ain
x1 = (-nx/2*dx2 + dx2/2 : dx2 : nx/2*dx2 - dx2/2);
y1 = (-ny/2*dy2 + dy2/2 : dy2 : ny/2*dy2 - dy2/2);
[X1, Y1] = meshgrid(x1 *dx,y1 *dy);

%Phase of matrix Ain and FFT of Ain
AFtPart = Ain.*exp(1i*k*z + (1i*k/(2*z)).*(X1.^2 + Y1.^2)); 
AFt = fftshift(fft2(ifftshift(AFtPart)));

%"coordinates" of matrix Aout
x2 = (-nx/2*dx2 + dx2/2 : dx2 : nx/2*dx2 - dx2/2);
y2 = (-ny/2*dy2 + dy2/2 : dy2 : ny/2*dy2 - dy2/2);
[X2, Y2] = meshgrid(x2, y2);

%compare Sf(x,y) in formula
AFoldCore = (1i*k /(2*pi*z)).*exp((1i*k/(2*z)).*(X2.^2 + Y2.^2));
Aout = AFoldCore.*AFt;
end

The formula this algorithm is based on:

Image 1: Formula for diffraction

My result for the amplitude of Aout is as expected, the problem lies within the phases, they are not symmetrical (obvious asymmetrical part marked in the picture):

Image 2: Amplitude and Phase of Aout

So I tried to find out where this asymetrical part came into existence. AFoldCore proved to be symmetric, AFtPart appears to be symmetric as well, but the phases of AFt itself are not symmetric:

Image 3: Amplitude and Phase of AFt

It appears as if a Sine with a rather low frequency was mulitplied into the phase. It ´"starts" in the top left corner and goes down to the bottom right corner.

After some try and error I decided to use the same algorithm for a Matrix with an odd size:

>> A = zeros(257,257);
>> A(119:139, 119:139) = 100;
>> fresnelpropagation(A, 100, 0.682, 0.39, 0.39);

The phase of Aout is symmetric, furthermore the phase of AFt is symmetric as well. Because I need to use the algorithm a lot of times in my program I intend to work with matrixes, with a size of N x N, whereas N is a power of 2 for perfomance reasons. Due to that simply "upsacling" to N+1 x N+1 is not a possible solution.

I think I know where the error is to be found: The ifftshift is not able to flip along a Pixel with the coordinates [0,0] if N is even. Due to that somehow the sine enters the phase.

My question is how am I able to surpress this sine in my phase without upscaling my Matrix? Or is there something wrong with my code?

Thank you for your time and attention.

Best regards,

Adrian

Answer 1

[David Goodmanson]

Hello Adrian, Could you live with a power-of-2 fft and an odd size square? Changing 119:138 to 120:138 puts the center of the square exactly at (129,129) and improves things, although the result is still not quite symmetric. Then changing

x1 = (-nx/2*dx2 + dx2/2 : dx2 : nx/2*dx2 - dx2/2)
to
x1 = (-nx/2*dx2         : dx2 : nx/2*dx2 - dx2  );

and a similar change in three other places puts the zero of coordinates (after fftshift etc.) where the fft thinks it is. The result is about as symmetric as I think you are likely to get.

[Adrian Schwarzenberger]

Thank you for your answer, I changed my code according to your suggestions and the result improved quite a bit. When using this code in the context of my bigger setup the asymmetrical part is still visable but the phases aren't completely random anymore and closer to what I expected them to be like.

Answer 2

[科秀]

There is a mistake in the diffraction formula, omitting the phase term exp(ikz) inside the FT(...).Frenel diffraction

[David Goodmanson]

I agree, you don't need exp(ikz) twice