Phase spectrum of 2D real function

Question

I am trying to calculate the phase spectrum of a 2D function, $f(x,y)$, where $x$, $y$, and $f(x,y)$ are real, using MATLAB's fft2. Thus my operation is $f(x,y)\rightarrow F(k_x,k_y)$. My function is such that $k_x$ only has positive values and $k_y$ has both positive and negative values.

Since my function is real, I only need the one-sided spectrum. I can calculate the one-sided amplitude spectrum by doubling the spectrum where $k_x>0$. The result I get is exactly what I expect.

When it comes to the phase spectrum, I get a function that is not continuous. The phases appear to jump around with an approximate difference of $\pi$. Is this related to phase wrapping? How do I get a continuous phase plot in positive $k_x$?

More information on my image and code: Data can be found here

My signal in $x$ and $y$ looks like this:

figure; contour(y,x,eta_0)
xlabel('y')
ylabel('x')

The amplitude spectrum looks like this:

dx = mean(gradient(x));
dy = mean(gradient(y));
Nx = 64;
Ny = 32;

Fsy = 1/dy;
Fsx = 1/dx;
kx_fft = 2*pi * Fsx/2 * linspace(0, 1, size(eta_0,1)/2+1);
ky_fft = 2*pi * Fsy/2 * linspace(0, 1, size(eta_0,2)/2+1);
ky_fft = [-fliplr(ky_fft(2:end-1)) ky_fft];

amp_xy = circshift(abs(fftshift(fft2(eta_0)))/Nx/Ny,-1,2);
figure; plot(ky,kx,amp_xy)

I know that my function does not contain any negative $k_x$ values and hence I can take the single-sided amplitude spectrum by doubling the positive $k_x$ values, to get this:

amp_xy = 2*amp_xy([Nx/2+1:Nx 1],:);
figure; imagesc(ky_fft,kx_fft,amp_xy);

The phase spectrum looks like this: phase_xy = circshift(angle(fftshift(fft2(eta_0))),-1,2); This clearly displays jumps from around zero (which is where I expect my phases to be) to $\pi$ and $-\pi$. Why is that and how can I avoid it?

To illustrate my point, I can take a cut through my signal at the centreline of $y$ and check the 1D amplitude and phase spectra. The 1D phase spectrum shows most of the phases close to zero, with the lowest frequencies having a phase close to $\pi$/$-\pi$ as I expected (from theory):

Now if I was to take a cut through the 2D phase spectrum where $k_y=0$, I see something similar, yet not continuous. It is as if I need to shift some phases by $\pi$ up or down:

In essence, the question is: how do I get a smooth 2D phase spectrum in positive values of $k_x$?

Edit: If I calculate the phase angle using atan, rather than angle, I get a smoothly varying phase plot. That is expected though, since atan limits angles to $-\pi$ to $\pi$ whereas angle uses atan2, which spans the entire circle.

Answer 1

A phase measurement needs to be referenced to some point.

The standard FFT computes a circular phase relative to the first sample, or the sample at the 2D origin. For data that is not exactly integer periodic in the FFT aperture (typically true for any rectangularly windowed data), this is at a circular discontinuity in the data vector. Thus the result is an alternation of FFT phase results to represent this edge discontinuity.

An edge discontinuity can be represented by an odd function, and odd functions show up in the imaginary (or sine) components of a DFT of this discontinuity.

You can move this circular or edge discontinuity away from the FFT’s phase zero reference point by doing an fftshift before performing the FFT. This moves the phase reference point to the center of your data, which hopefully does not have an unintended discontinuity. This shift most often results in a more continuous phase spectrum across FFT result bins of windowed data. However a pre-fftshift produces the same result as post flipping the phase of all the odd numbered FFT result bins without the shift, which may be easier to do.