Pi phase shift b/w adjacent pixels while using fft2

Question

I have a 2D Gaussian beam at its waist (where the phase is zero throughout the transverse plane). When I use fft2 to find the 2D spatial Fourier transform and plot the phase, I observe that there is a pi phase shift between any 2 adjacent data points. However, I don't observe this when I use for loops to calculate the Fourier transform instead of using fft2.

Is this due to phase wrapping? How do I overcome this?

Thanks.

Edit: I'm posting the code for the fft of a circular aperture, since the same results are observed, and because it is much simpler.

Nx = 200; Ny = Nx;

%creating coordinate grids
x = -Nx/2:Nx/2 - 1; y = -Ny/2:Ny/2 - 1;
[X,Y] = meshgrid(x,y);

r = 15; %radius of aperture

Eip = ((X.^2 + Y.^2 ) <= r^2); %aperture

figure;pcolor(abs(Eip));axis square; shading flat; colorbar;
figure;pcolor(angle(Eip));axis square; shading flat; colorbar;

Cip = fftshift(fft2(Eip)); %FFT

figure;pcolor(abs(Cip));axis square; shading flat; colorbar;
figure;pcolor(angle(Cip));axis square; shading flat; colorbar; 

Answer 1

Short answer

Cip turns out to be real. You're just seeing sign changes between contiguous points.

Long answer

Eip is obviously real. Besides, it exhibits the following symmetry along the x and y axes. Let's take the x axis first. For any fixed y, consider Eip as a finite signal defined at x values 0,1,...,Nx-1. If that finite signal is extended into a periodic sequence, then this periodic sequence turns out to be even. This is true for your specific definition of Eip. Together with the fact that the values are real, this implies¹ that the DFT along x is also real and has the same type of symmetry. Now, regarding the y axis, the same holds. The final result is that the 2D-DFT Cip is real and has the mentioned x and y symmetries.

(As a side note, Cip as obtained in your code is not real but complex. However, the imaginary part is of the order of eps, and thus negligible. You can check:

>> max(real(Cip(:)))
ans =
   709
>> max(imag(Cip(:)))
ans =
  4.6262e-014

So the imaginary part is only an artifact of finite numerical precision. We can do

Cip = real(Cip);

to remove that spurious imaginary part.)

Now that we know that Cip is real (up to numerical precision), it's natural that all its values have phase either 0 or pi (or equivalently -pi). That is, what you are seeing is just changes of sign between consecutive values in the DFT. To illustrate, see the following graph, which shows the variation of Cip along a line parallel to the y axis (corresponding to x equal to 120; just as an example):

stem(1:Nx, Cip(120,:))

^{1 See Discrete-time signal processing, Oppenheim et al., 2nd ed., pp. 568--570.}

Answer 2

To avoid this, you need an additional ifftshift. This is due to the algorithm Matlab uses and the results you expect. To avoid the pi phase jumps use this instead:

fftshift(fft2(ifftshift(Eip)));