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:
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.)
, 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.
