3
$\begingroup$

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

Signal in x and y

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)

Amplitude Spectrum

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

Amplitude spectrum in positive k_x

The phase spectrum looks like this: phase_xy = circshift(angle(fftshift(fft2(eta_0))),-1,2); Phase spectrum in positive k_x 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): 1D Phase

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:

Cut from 2D phase spectrum

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.

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ You are right. I should have thought of it: for the 1D case I shift the maximum of the signal at the beginning of my array to adjust the phases, so it makes sense that for the 2D signal, I should also shift the maximum at the (0,0) location. $\endgroup$ – Mitsarien Jan 20 '18 at 13:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.