2D phase of Fourier transform is not behaving as expected

Hello. I am trying to compute the 2D FFT of a very simple image, a 2000x2000 matrix of zeros with a small 100x100 square of 1s in the center. I can obtain the complex amplitude correctly (a 2D sinc function), but I am having trouble when computing the phase. First, I have seen that I obtain a phase that is wrapped between -pi and pi. I have therefore tried unwrapping the phase first along rows and then unwrap that along the columns as well. Still, I obtain a phase which looks like a 45 degrees gradient. My understanding is that the phase of a centered square should be zero across the entire image. My suspicion is that MATLAB considers the zero frequency component of the image to be on the top left corner, and a centered rectangle would mean a shift in both vertical and horizontal direction. I have tried applying fftshift to get rid of this problem, but it persists. I think I have tried almost all combinations, using fftshift, ifftshift, unwrapping in different orders, etc... No matter what I do, I do not get the expected result. What can I try?

Here is my code for the complex phase:

canvas = zeros(2000,2000); % Create empty matrix of 2000x2000 zeros
canvas(950:1050,950:1050)=1; % Place 100x100 rectangle in the center

% Calculate 2D FFT

canvas_FFT2=(fft2(canvas)); % Take the bidimensional FFT of the canvas, shift it such that it's in the center of the image with fftshift
canvas_FFT2_ampl=abs(fftshift(canvas_FFT2)); % Extract the complex amplitude
canvas_FFT2_phase=(angle(fftshift(canvas_FFT2))); % Extracts the complex phase
canvas_FFT2_phase1=unwrap(canvas_FFT2_phase,[],1); % Unwrap it along the rows
canvas_FFT2_phase=unwrap(canvas_FFT2_phase1,[],2); % Unwrap the result along the columns to complete the 2D unwrapping

% Display results

figure(1)
subplot(1,2,1) % Plot the real space on the left in figure 1
imagesc(canvas)
colormap parula
colorbar
title("Real space")
xlabel("Pixel x")
ylabel("Pixel y")
axis equal
xlim([1,2000])
ylim([1,2000])

subplot(1,2,2) % Plot the reciprocal space on the right (complex phase) in figure 1
imagesc((canvas_FFT2_phase));
colormap parula
colorbar
title("Reciprocal space (complex phase)")
xlabel("Pixel x")
ylabel("Pixel y")
axis equal
xlim([1,2000])
ylim([1,2000])

Thank you very much!

Answer 1

[Daniel]

There are two issues compounding here to make for the phase shift you're seeing.

Issue 1: Careful alignment with the middle of the FFT. The symmetry has to be about the first element, so your construction from elements 950 to 1050 on a length-2000 axis is not quite correct. Ignoring the first element, the next 948 elements are zero on the low side, and then the last 950 elements are zero on the high side. Those numbers need to match. So part 1 of the solution is to set

canvas(951:1051,951:1051)=1

in line 2. This is the only real issue, but after you apply that change you will still see a strange nonzero unwrapped phase.

This is because negative numbers have a phase of depending on noise in the imaginary component, rather than . Since fft2 operates by running fft along each axis in turn, take a look at the figures below:

pulse = zeros(20,1);
pulse(6:16) = 1;

PULSE = fft(pulse);

subplot(2,2,1)
plot(pulse)
axis padded
title('Time domain')
subplot(2,2,2)
plot(abs(PULSE))
title('PSD-ish')
subplot(2,2,3)
plot(real(PULSE))
hold on
plot(imag(PULSE))
title('Real-imag DFT')
subplot(2,2,4)
plot(angle(PULSE))
title('Wrapped phase of DFT')

Again, when the real part is negative and the imaginary part is , the phase will be rather than 0. (The sign is determined by computational noise in the imaginary part.)

Now let's look at that in two dimensions.

canvas = zeros(40,40);
canvas(11:31,11:31) = 1;

CANVAS = fft(canvas);
CANVAS2 = fft2(canvas);

figure
imagesc(canvas)
title("Signal")

figure
imagesc(real(CANVAS))
colorbar
title("DFT along each column, real component")

figure
imagesc(real(CANVAS2))
colorbar
title("2-D DFT, real component")

All those negative real components lead to non-fixed phases, which lead you to a somewhat strange-looking unwrap pattern.

figure
imagesc(angle(CANVAS2))
title("Wrapped phases")
colorbar

figure
imagesc(unwrap(angle(CANVAS2),[],1))
title("Unwrapping along one dimension")
colorbar

figure
imagesc(unwrap(unwrap(angle(CANVAS2),[],1),[],2))
title("Unwrapping along both dimensions")
colorbar

You should see that same sort of streaky pattern along your 2000x2000 image once you've shifted the indices to get rid of the deterministic phase gradient you were seeing before.

[Andrea]

Thank you very much for this answer, it was very elucidating! I am wondering if there is any way to fix this artifact (i.e. looping through the imaginary part and change when the phase is -pi or pi to 0?).

[Daniel]

@Andrea, there are a number of ways to "fix" this (though it's not really a problem, it's more of a mathematical inconvenience :) ).

1. You could do as you suggest, and manually force anything to 0.
2. A slightly more refined approach would be to look at the real part of the fft2 output. If the real part is negative, then add π to the angle if the angle is negative, and subtract π if the angle is positive. That would give you a bunch of numbers near 0.

But either approach would still look funny when subjected to imagesc, since it will happily normalize to cover the full range of colors across the range of your input, even if the input range is tiny.

Which leads to

3. You could just show that all the imaginary components are near zero directly, and forget all about phases and unwrapping.

(By the way, there's no need to loop, since MATLAB is perfectly happy to operate on vectors with logical indexing. Here's an example of that for the loop checking. Of course in your example, since the values should all be very near 0 or , the "fixed" values should all be extremely near 0.)

rng default
a = randn(10,2)*[1;j];
angles = angle(a);
figure
stem(angles)
angles(angles<-pi/2) = angles(angles<-pi/2)+pi;
angles(angles>pi/2) = angles(angles>pi/2)-pi;
hold on
stem(angles)
yline(-pi/2)
yline(pi/2)
legend("Before casting to [-\pi,\pi] range","After casting to [\pi,\pi] range","","")

Answer 2

[ZQ]

Use odd*odd matrix instead of even*even.(i.e. 2001*2001)

[Andrea]

I have tried now to change the matrix size but the result is unchanged...