0
$\begingroup$

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?

enter image description here

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);
shading interp
axis tight
title ('Rect 3D');
rect=ifftshift(rect);
figure, surf(X,Y,rect);
shading interp
axis tight
title ('Rect 3D shifted');

R = fft2((rect));
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R)); 
shading interp
axis tight
title('Fourier Transform of Rectangular function');
%plot real part
figure;
surf(X,Y,real(R)); 
shading interp
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);
shading flat
axis tight
title ('Phase of the rect');
$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – Vittorio Todisco Jun 8 '17 at 19:15
  • $\begingroup$ 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 $\endgroup$ – Vittorio Todisco Jun 8 '17 at 19:19
  • $\begingroup$ 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. $\endgroup$ – Peter K. Jun 8 '17 at 19:23
  • $\begingroup$ 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. $\endgroup$ – Vittorio Todisco Jun 8 '17 at 19:47
  • $\begingroup$ s24.postimg.org/irsdwgz51/linearshifting.jpg $\endgroup$ – Vittorio Todisco Jun 8 '17 at 19:48
0
$\begingroup$

@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:

enter image description here

enter image description here

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));
shading interp
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)); 
shading interp
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);
shading flat
axis tight
title ('phase of rect pulse');
$\endgroup$
0
$\begingroup$

The checkerboard pattern is due to not using an fftshift. Since your pulse is not centered at FFT index (0,0), the beginning of the FFT array, which is different from spatial coordinate (0,0), all your phases are twisted (alternately negated for a half aperture shift from center).

$\endgroup$
0
$\begingroup$

Unfortunately I don't use MatLab so can't help with that code. I did quickly put something together with my own code that seems close to what you were aiming for. A square pulse of width 10 centered at the origin within a 512 field. The attached is my magnitude and phase spectra. enter image description here

enter image description here

$\endgroup$

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.