How to workaround this problem with angle?

6 次查看(过去 30 天)
(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

回答(2 个)

Chunru
Chunru 2024-3-5
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
  1 个评论
SIMONE
SIMONE 2024-3-5
编辑:SIMONE 2024-3-5
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

请先登录,再进行评论。


VBBV
VBBV 2024-3-5
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
  2 个评论
SIMONE
SIMONE 2024-3-5
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
VBBV 2024-3-5
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

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

产品


版本

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by