How to workaround this problem with angle?

(I started with MATLAB today, so I'm new). I created a script to calculate the fourier transform of a function, the amplitude spectrum and the phase spectrum. For testing purposes i tried calculating using as the input function the rectangular impulse. The fourier transform and the amplitude spectrum are right, the phase spectrum is not. When i use the angle function or the atan2 to calculate the angle i get as output alternating +pi and -pi where the output should be +pi. Watching the inputs and outputs of the angle function everything seems fine, the angle function works as intended. I think that the problem is related to the sign of the imaginary part, and maybe to the fact that they are really small number (e^-18) . I know that the code is all based on an approximation so it could be normal but i want to know if there is a workaround. (I know there is already a function(fft))

%total range of omega
delta_omega=20;
%intervals of omega that are used to approximate
omega_precision=0.01;
%initializatoin of the results vector
res=0:omega_precision:delta_omega;
%all the values of omega
val_omega=-delta_omega/2:omega_precision:delta_omega/2;
%for every omega calculate it's fourier transform (from -delta_omega/2 to delta_omega/2)
for n=1:delta_omega/omega_precision+1
    res(n)=X((n-(delta_omega/omega_precision/2))*omega_precision);
end
%plot of the real value of the fourier tranform
plot(val_omega,real(res),"red");
hold on;
%plot og the immaginary values of the fourier tranform
plot(val_omega,imag(res),"cyan");
%amplitude spectrum values
V=abs(res)/pi;
hold on;
%plot of the amplitude spectrum
plot(val_omega,V,"green");
xlabel("omega");
%phase spectrum values
Fi=-angle(res);
%plot of the phase spectrumm values
plot(val_omega,Fi,"blue")

%fourier transform of function x(t)
function X=X(omega)
    handle= @(t) x(t).*exp(-1i*omega*t);
    X=integral(handle,-1100,1100);
end
%function x(t) (rectangular impulse)
function x = x(t)
    %amplitude
    A=1;
    tao=pi;
    x=1:length(tao);
    for j=1:length(t)
        if(t(j)<=tao/2 && t(j)>=-tao/2)
            x(j)=A;
        else
            x(j)=0;
        end
    end
end

Answer 1

[Chunru]

You can try unwrap function.

%total range of omega
delta_omega=20;
%intervals of omega that are used to approximate
omega_precision=0.01;
%initializatoin of the results vector
res=0:omega_precision:delta_omega;
%all the values of omega
val_omega=-delta_omega/2:omega_precision:delta_omega/2;
%for every omega calculate it's fourier transform (from -delta_omega/2 to delta_omega/2)
for n=1:delta_omega/omega_precision+1
    res(n)=X((n-(delta_omega/omega_precision/2))*omega_precision);
end
%plot of the real value of the fourier tranform
plot(val_omega,real(res),"red");
hold on;
%plot og the immaginary values of the fourier tranform
plot(val_omega,imag(res),"cyan");
%amplitude spectrum values
V=abs(res)/pi;
hold on;
%plot of the amplitude spectrum
plot(val_omega,V,"green");
xlabel("omega");
%phase spectrum values
Fi=-angle(res);
%plot of the phase spectrumm values
% plot(val_omega,Fi,"blue")
plot(val_omega,unwrap(Fi),"blue")

%fourier transform of function x(t)
function X=X(omega)
    handle= @(t) x(t).*exp(-1i*omega*t);
    X=integral(handle,-1100,1100);
end
%function x(t) (rectangular impulse)
function x = x(t)
    %amplitude
    A=1;
    tao=pi;
    x=1:length(tao);
    for j=1:length(t)
        if(t(j)<=tao/2 && t(j)>=-tao/2)
            x(j)=A;
        else
            x(j)=0;
        end
    end
end

[SIMONE]

I have 2 problems:

1) I already tried to unse unwrap and my output was different from yours, does some know why it looke like this when i run it?

2) This would be wrong, maybe it's my fault and i did something else wrong, as the phase spectrum should look like this from what i found online:

Thanks for rsponding and trying ti help me

Answer 2

[VBBV]

You have real and imaginary componens for variable Fi, so use only real angles

%total range of omega
delta_omega=20;
%intervals of omega that are used to approximate
omega_precision=0.01;
%initializatoin of the results vector
res=0:omega_precision:delta_omega;
%all the values of omega
val_omega=-delta_omega/2:omega_precision:delta_omega/2;
%for every omega calculate it's fourier transform (from -delta_omega/2 to delta_omega/2)
for n=1:delta_omega/omega_precision+1
    res(n)=X((n-(delta_omega/omega_precision/2))*omega_precision);
end
%plot of the real value of the fourier tranform
plot(val_omega,real(res),"red");
hold on;
%plot og the immaginary values of the fourier tranform
plot(val_omega,imag(res),"cyan");
%amplitude spectrum values
V=abs(res)/pi;
hold on;
%plot of the amplitude spectrum
plot(val_omega,V,"green");
xlabel("omega");
%phase spectrum values
Fi=angle(real(res));
%plot of the phase spectrumm values
plot(val_omega,Fi,"blue")

%fourier transform of function x(t)
function X=X(omega)
    handle= @(t) x(t).*exp(-1i*omega*t);
    X=integral(handle,-1100,1100);
end
%function x(t) (rectangular impulse)
function x = x(t)
    %amplitude
    A=1;
    tao=pi;
    %x=1:length(tao);
    for j=1:length(t)
        if(t(j)<=tao/2 & t(j)>=-tao/2)
            x(j)=A;
        else
            x(j)=0;
        end
    end
end

[SIMONE]

This would work for this specific function x(t), another solution for that would be abs(angle(res)). The problem is that it is not generalized, if i try using a function x(t) that has the imaginary part in the fourier transform then the phase spectrum is wrong. An example is using x(t)={A*exp(-abs(t)/t_0) for t>=0 , 0 for t<0}. With this funtion -angle(res) gives the right phase spectrum and angle(real(res)) equals 0. Thanks for your response and your time.

function x = x(t)
%amplitude
A=1;
t_0=1;
for j=1:length(t)
    if t(j)>=0 
        x(j)=A*exp(-abs(t(j))/t_0);
    else
        x(j)=0;
    end
end
end

[VBBV]

If you consider the imaginary part of Fourier spectrum, then use the imag(res) . For a real periodic signal, the phase information of imaginary component of signal has no significance

%total range of omega
delta_omega=20;
%intervals of omega that are used to approximate
omega_precision=0.01;
%initializatoin of the results vector
res=0:omega_precision:delta_omega;
%all the values of omega
val_omega=-delta_omega/2:omega_precision:delta_omega/2;
%for every omega calculate it's fourier transform (from -delta_omega/2 to delta_omega/2)
for n=1:delta_omega/omega_precision+1
    res(n)=X((n-(delta_omega/omega_precision/2))*omega_precision);
end
%plot of the real value of the fourier tranform
plot(val_omega,real(res),"red");
hold on;
%plot og the immaginary values of the fourier tranform
plot(val_omega,imag(res),"cyan");
%amplitude spectrum values
V=abs(res)/pi;
hold on;
%plot of the amplitude spectrum
plot(val_omega,V,"green");
xlabel("omega");
%phase spectrum values
Fi= angle(imag(res));
%plot of the phase spectrumm values

plot(val_omega,Fi,"blue")

%fourier transform of function x(t)
function X=X(omega)
    handle= @(t) x(t).*exp(-1i*omega*t);
    X=integral(handle,-1100,1100);
end
%function x(t) (rectangular impulse)
function x = x(t)
%amplitude
A=1;
t_0=1;
for j=1:length(t)
    if t(j)>=0 
        x(j)=A*exp(-abs(t(j))/t_0);
    else
        x(j)=0;
    end
end
end