Why the calculation error enhances as the phase angle increases when using atan2?

%% 
clear all;clc;
alpha = 0:pi/10:100*pi;
phi = atan2(sin(alpha),cos(alpha));
phi_up = unwrap(phi);
%%
figure
plot(phi_up,'bo--')
grid on;
figure
plot(alpha,'r')
grid on;
figure
plot(phi_up - alpha,'r')
grid on;

Answer 1

[dpb]

alpha1 = 0:pi/10:100*pi;
alpha2=linspace(0,100*pi,numel(alpha1));
plot(abs(alpha1-alpha2),'r-')

It's not where you think it is -- the problem is in how you generated your alpha vector; the colon operator rounded pi/10 and used that approximation repetitively as n*delta. linspace() ototh, calculated the closest approximation to the locations between the start and end points with the same number of elements in the two vectors.

So, it's actually the alpha values that diverge from best approximation, not atan2 being less accurate; it reproduces the input it is given as best approximation; if that value is not quite what is expected, it's the input that is to blame, not atan2

[zhouqing wang]

Thank you very much for your reply!

However when using the linspace() , the result is still the same that the calculation error enhances as the phase angle increase.

%%
clear all;clc;
alpha1 = 0:pi/10:100*pi;
alpha2=linspace(0,100*pi,numel(alpha1));
% plot(abs(alpha1-alpha2),'r-')

    phi = atan2(sin(alpha2),cos(alpha2));

    phi_up = unwrap(phi);

%%
figure
plot(phi_up,'bo--')
grid on;

figure
plot(alpha2,'r')
grid on;

figure
plot(phi_up - alpha2,'r')
grid on;

In the other way, my key question is How to make the error convergence when using atan2() and unwrap()?

[dpb]

The problem of argument reduction from multiples of hundreds of times the fundamental period of the trigonometric functions has been discussed on Answers before I remember, but I couldn't seem to find the particular threads that talk about it specifically.

The problem is the trig functions are computed over a range of something like [0 pi/4] from specialized routines using combinations of lookup and various approximations; to compute them out side of that range requires "range reduction" to reduce the input argument into the range actually used by the algorithm. This introduces an error besides just the inherent issues of floating point precision and rounding.

To see what happens, compare the full-blown series of the sin() to that computed over the first period--let the code return the full vector and then build a second version that repeats the first 0-2pi section as many times as needed. The following figure shows that you can't see the difference in actual value on a plot, but if you subtract the two, then the differences show up,

alpha1 = 0:pi/10:100*pi;
alpha2=linspace(0,100*pi,numel(alpha1));
phi = atan2(sin(alpha2),cos(alpha2));
s1=sin(alpha2); c1=cos(alpha2);
s2=[repmat(s1(1:20),1,50), s1(11)];
subplot(2,1,1)
plot(alpha2.',[s1;s2].'), legend('S1','S2')
xlim([0 100*pi])
title('sin over 0-2pi and 0-100pi')
subplot(2,1,2)
plot(alpha2.',[s1-s2].')
xlim([0 100*pi])
title('Difference between sin 0-2pi and 0-100pi')

I'm sorry, but I don't understand the follow-up question to know how to try to address it. What are you trying to converge on? The magnitude of the difference is still in the range of multiples of the resolution of a double precision and so would rarely, if ever, be actually observed.

Can you illustrate where you're having a problem, specifically?

[zhouqing wang]

Thank you for your reply !

To present my question more clearly, I want to make the difference between "intial setting phase angle values "and "the phase angle results of (atan2+unwrap) processing " as small as possible. However when I run the code as below, the difference between them become larger as the phase angle values increase, like a ladder as shown in the figure below. Why it like this? How to make the difference small and convergence?

%% 
clear all;clc;
alpha = 0:pi/10:1000*pi;

phi = atan2(sin(alpha),cos(alpha));
phi_up = unwrap(phi);

%%
figure
plot(phi_up,'b')
grid on;

figure
plot(alpha,'r')
grid on;

figure
plot(phi_up - alpha,'r')
grid on;

Answer 2

[Bruno Luong]

It doesn't seem weird to me.

Both of the (atan2+unwrap) and (sin+cos) do somekind of undo modulo 2*pi differently.

unwrap makes a cumulative correction by adjusting 2*pi when there is "wrong" jump multiple time, on other hand sin/cos perhaps do modulo (2*pi) by remainder after division.

I don't know the detail but sin and cos probably come from processor FPU implementation and the wrapping by 2*pi is probably carried out in smaller intervals using trigonometric identities with higher internal precision and a multiple of such implementation details.

So it is normal that the error increases. Note that 2*pi cannot be exatctly represented by binary-bases floating point number.

It is hard to know which is more acurate (my guess is sin/cos). I don't think one can do anything to match methods, beside replicate exactly what the FPU does, which is a daunting task.

[Bruno Luong]

Actually the match is not bad. Only the least significant bit has a jitter as showed in this plot. So nothing to worry about (see eps function)

%% 
nround = 1e6;
alpha = 0:pi/10:nround*pi;
phi = atan2(sin(alpha),cos(alpha));
phi_up = unwrap(phi);
%%
err = alpha - phi_up;

close all
figure
plot(err,'.r')
hold on
plot(eps(alpha),'g')
grid on;
legend('error = (alpha - phi\_up)', 'eps(alpha)', 'location', 'northwest')

[dpb]

Yeah, that's what was pointing out to OP that it's on order of rounding and with floating point precision is same number of digits so as the magnitude of the argument goes up by orders of magnitude, the absolute accuracy simply has to go down; rounding is a larger numeric value but the same relative magnitude.

[Bruno Luong]

@dpb I did not read careful your post, I started to read then I see you mentioned linspace vs colon. I have no idea why you compare with linspace that doesn't involve at all in OP question.

Same phenomenal occurs regardless colon, linspace, logspace or whatever command to generate alpha with large value

nround = 1e6;
alpha = logspace(0,log10(nround*pi),nround*20+1);
phi = atan2(sin(alpha),cos(alpha));
phi_up = unwrap(phi);
%%
err = alpha - phi_up;

figure(1)
plot(err,'.')
hold on
semilogx(eps(alpha),'r')
grid on;
legend('error (alpha - phi\_up)', 'eps(alpha)', 'location', 'northwest')

Now that I see this curve I think the cause is not internal/external modulo, not alpha but because atan2(sin(alpha),cos(alpha)) does not return exactly alpha even for alpha in (-pi, pi) - that is expected. If I do the same test in this range because that does not require any unrapping for comparison of angle. Below is the error plot. The error is again either 0 or 1 LSB.

figure(2)
nround = 1e6;
alpha = linspace(-pi,pi,1e6);
alpha([1 end]) = [];
phi = atan2(sin(alpha),cos(alpha));
%%
err = alpha - phi;

figure
plot(abs(err),'.')
hold on
plot(eps(alpha),'r')
grid on;
legend('error (alpha - phi)', 'eps(alpha)', 'location', 'best')

For large range, unwrap just expands the LSB error to the range of alpha, thus the bit error is eps(alpha) and increases with alpha. That theory alone explains evrything.

[dpb]

Yes, two different issues -- colon and linspace end up with different approximations to alpha but the same effect is present with both regarding the unwrapping.