MATLAB phase of 2D rectangular pulse's Fourier transform

I'm trying to plot the graph of the phase of the Fourier transform of a 2D rectangular pulse. I've been able to evaluate the FFT but I'm not sure if the phase is correct because there are some tilts I can't explain.

• What are those diagonal lines that appear in the phase?
• How can I eliminate them?

Here there is the code:

clc
close all
clear
npoints=512;
perc=1;
dt=6*1E-7/(npoints*perc);  % Tempo di campionamento
df=1/(npoints*dt); % Frequenza di campionamento

t(1)=0;
f(1)=0;
for k=2:npoints/2
t(k)=(k-1)*dt;
t(npoints-k+2)=-t(k);
f(k)=(k-1)*df;
f(npoints-k+2)=-f(k);
end
t(npoints/2+1)=t(npoints/2+2)-dt;
f(npoints/2+1)=f(npoints/2+2)-df;

ts=ifftshift(t);
fs=fftshift(f);

figure
[X,Y]=meshgrid(ts,ts);

D = npoints/2;        % to indicate origin at the center of the function
a = 5;          % change it to enlarge or reduce the pulse
y = repmat(1:npoints,npoints,1);
x = y';
rect = zeros(npoints);
rect(D-a:D+a-1,D-a:D+a-1) = ones(2*a);
rect=(rect);
surf(X,Y,rect);
axis tight
title ('Rect 3D');
rect=ifftshift(rect);
figure, surf(X,Y,rect);
axis tight
title ('Rect 3D shifted');

R = fft2((rect));
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R));
axis tight
title('Fourier Transform of Rectangular function');
%plot real part
figure;
surf(X,Y,real(R));
axis tight
title('Real part');

Rm=abs(R);
imm=imag(R);
re=real(R);
re(abs(re) < 1e-8) = 0;
imm(abs(imm) < 1e-8) = 0;
R=re+imm;
phase=angle(R);

%plot phase
[X,Y]=meshgrid(fs,fs);
figure
surf(X,Y,phase);
axis tight
title ('Phase of the rect');


The issue is that the phase value found is between $-\pi$ and $\pi$ (or $0$ and $2\pi$) but it needs to be "unwrapped" to be continuous.

In 1D, the unwrap function will help. Your mileage may vary when applying it to a 2D Fourier transform.

Also, there is a bug in your code that forces the R variable to be real valued:

R=re+imm;
phase=angle(R);


is not going to get you the right phase at all.

• The phase only assumes two values, 0 and $\pi$ which is what I expect from a sinc function. s17.postimg.org/jmiz1v39r/0pi.jpg But I don't get those diagonal lines... Jun 8 '17 at 19:15
• But rising the parameter "a" to almost npoints/2 and zooming on the top of the shifted gaussian time domain pulse I found out that the rect is missing some points. Could it be the cause? s17.postimg.org/da3ts107j/zoom_on_top.jpg Jun 8 '17 at 19:19
• You're not doing anything to account for fft2's 0 to $N-1$ / 0 to $M-1$ axes, so I can't see how you'd expect the FFT result to be purely real. Ah. R=re+imm; phase=angle(R); means you're forcing the R FFT to be real-valued. Don't do that.
– Peter K.
Jun 8 '17 at 19:23
• You're right i forgot to multiply for the imaginary unit... I fixed that but know it seems that the phase is multiplied for a linear shift. Jun 8 '17 at 19:47
• s24.postimg.org/irsdwgz51/linearshifting.jpg Jun 8 '17 at 19:48

@peter I solved the problem by changing the method to create the rectangular pulse, I just used the rectpulse matlab function. This is the phase of the pulse:

And this is the working code:

clear
clc
close all

npoints=256;
perc=1;
dt=6*1E-7/(npoints*perc);  % time of sampling
df=1/(npoints*dt); % Sampling frequency

t(1)=0;
f(1)=0;
for k=2:npoints/2
t(k)=(k-1)*dt;
t(npoints-k+2)=-t(k);
f(k)=(k-1)*df;
f(npoints-k+2)=-f(k);
end
t(npoints/2+1)=t(npoints/2+2)-dt;
f(npoints/2+1)=f(npoints/2+2)-df;

ts=ifftshift(t);
fs=fftshift(f);

a = 5;          % change it to enlarge or reduce the pulse. Values from 1 to 50
T=a*1E-8;
rect=rectpuls(t,T);     %build 1D rect
rect = rect'*rect;      %build 2D rect

figure;                 %plot 2D rect
[X,Y]=meshgrid(ts,ts);
surf(X,Y,ifftshift(rect));
axis tight
title ('3D rect pulse');

R = fft2(rect);         %Fourier transform of the pulse
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R));
axis tight
title('Fourier Transform of Rectangular function');

% Modulo
Rm=abs(R);
imm=imag(R);
re=real(R);
re(abs(re) < 1e-12) = 0;
imm(abs(imm) < 1e-12) = 0;
R=re+1j*imm;
faserect=angle(R);

[X,Y]=meshgrid(fs,fs);
figure
surf(X,Y,faserect);