2017-10-17 137 views
4

我在寫一段代碼時遇到了問題,這段代碼是在集成過程中圍繞一個角度進行的,並且是我正在處理的一個小模擬的一部分。所以基本上這個想法是通過確保它始終具有理智的價值來防止角度變大。我已經嘗試了三種不同的方法,我期望得到相同的結果。他們大部分時間都是這樣。但是前兩個在角度值環繞的地方產生了僞影。當我從角度值生成波形時,由於這些精度錯誤,我會得到不理想的結果。atan2f vs fmodf vs只是簡單的減法

所以第一個方法是這樣的(極限角度-8PI + 8PI範圍):

self->state.angle = atan2f(sinf(angle/8), cosf(angle/8)) * 8; 

這產生僞影,看起來像這樣:

enter image description here

二的方法

self->state.angle = fmodf(angle, (float)(2.f * M_PI * 8)) 

創建相同的結果: enter image description here

但是,如果我只是做這樣的:

float limit = (8 * 2 * M_PI); 
if(angle > limit) angle -= limit;   
if(angle < 0) angle += limit;    
self->state.angle = a; 

然後,它工作正常,沒有任何文物: enter image description here

那我在這裏失蹤?爲什麼其他兩種方法會產生精度錯誤?我希望他們都能產生相同的結果(我知道角度的範圍是不同的,但是當角度進一步傳遞到sin函數時,我會期望結果是相同的)。

編輯:小試

// g++ -o test test.cc -lm && ./test 
#include <stdio.h> 
#include <stdlib.h> 
#include <math.h> 
#include <stdint.h> 

int main(int argc, char **argv){ 
    float a1 = 0; 
    float a2 = 0; 
    float a3 = 0; 
    float dt = 1.f/7500.f; 

    for(float t = -4.f * M_PI; t < (4.f * M_PI); t+=dt){ 
     a1 += dt; 
     a2 += dt; 
     a3 += dt; 

     float b1 = a1; 
     if(b1 > 2.f * M_PI) b1 -= 2.f * M_PI; 
     if(b1 < 0.f) b1 += 2.f * M_PI; 
     float b2 = atan2f(sinf(a2), cosf(a2)); 
     float b3 = fmodf(a3, 2 * M_PI); 

     float x1 = sinf(b1); 
     float x2 = sinf(b2); 
     float x3 = sinf(b3); 

     if((x1 * x2 * x3) > 1e-9){ 
      printf("%f: x[%f %f %f],\tx1-x2:%f x1-x3:%f x2-x3:%f]\n", t, x1, x2, x3, (x1 - x2) * 1e9, (x1 - x3) * 1e9, (x2 - x3) * 1e9); 
     } 
    } 

    return 0; 
} 

輸出:

-9.421306: x[0.001565 0.001565 0.001565],  x1-x2:0.000000 x1-x3:0.000000 x2-x3:0.000000] 
-9.421172: x[0.001431 0.001431 0.001431],  x1-x2:0.000000 x1-x3:0.000000 x2-x3:0.000000] 
-9.421039: x[0.001298 0.001298 0.001298],  x1-x2:0.000000 x1-x3:0.000000 x2-x3:0.000000] 
-9.420905: x[0.001165 0.001165 0.001165],  x1-x2:0.000000 x1-x3:0.000000 x2-x3:0.000000] 
-9.420772: x[0.001032 0.001032 0.001032],  x1-x2:0.000000 x1-x3:0.000000 x2-x3:0.000000] 
-6.275573: x[0.001037 0.001037 0.001037],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.275439: x[0.001171 0.001171 0.001171],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.275306: x[0.001304 0.001304 0.001304],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.275172: x[0.001438 0.001438 0.001438],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.275039: x[0.001571 0.001571 0.001571],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.274905: x[0.001705 0.001705 0.001705],  x1-x2:0.000000 x1-x3:174.855813 x2-x3:174.855813] 
-6.274772: x[0.001838 0.001838 0.001838],  x1-x2:0.116415 x1-x3:174.855813 x2-x3:174.739398] 
+1

'if(angle> limit)angle - = limit; '不是'while(angle> limit)angle - = limit; '所以基本上如果'angle = 800000 * M_PI'你的最後一個方法不會把你的值放在範圍內。一個[mcve]可能很有用,輸入值是預期的,而使用'fmod'的「意想不到」(忘記atan2f' ATM) –

+0

屏幕截圖看起來並不是在一臺機器上完成的,當你使用* double *而不是* float *。 –

+0

該代碼只能在32位浮點支持的嵌入式目標上運行。 – Martin

回答

3

沒有更多的信息很難作出解釋,但我想試試。

使用fmod和「普通減法」(或增加),很喜歡現在做的區別是,如果該值出路的範圍已經(如800000 * M_PI例如),然後加/減方法不更改數值(它幾乎沒有效果),並且非常大(以絕對值)角度觸及您的計算函數,沒有問題,因爲沒有看到僞影。

使用fmod(或atan2)可確保該值處於您定義的範圍內,這不是同一件事。

需要注意的是這樣做的:

float limit = (8 * 2 * M_PI); 
while(angle > limit) angle -= limit;   
while(angle < 0) angle += limit;    
self->state.angle = a; 

將相當於(大約)至fmod(但會比更糟糕fmod爲大值,因爲它引入了由於反覆的加法或減法的浮點積累誤差)。

所以如果在您的計算中輸入非常大的值會產生正確的結果,那麼您可能想知道標準化角度而不是將其留給數學庫是否明智。

編輯:這個答案的第一部分,假設這種超級出界外的情況下會發生,以及進一步的問題編輯表明,這種情況並非如此,所以......

另一個區別fmod和2之間的測試是,有沒有保證時調用fmod

舉例來說,如果實現是像value - int(value/modulus)*modulus;,浮點不準確性可能。減去較小值,以原始值的值在範圍相同的,如果已經。

使用atan2f結合sin ...如果已經在範圍內,也會更改結果。

(即使該值略微超出了範圍,添加/底塗喜歡你做不涉及分/截斷/乘)只需添加

既然你可以在範圍內調整值或subbing一次,使用fmodfatan2f在你的情況下是矯枉過正,你可以堅持你的簡單的子/添加(添加一個else會節省一個測試:如果你只是調整了一個太低的值,不需要測試,看看是否值太大)

+0

當然可以。我只是想弄清楚爲什麼不同的結果給出了上面的例子。目的是保持角度不會出現在無窮遠處,因此無需一次性將其捕捉到範圍內,因爲在每次迭代期間它永遠不會超過0.5PI。 – Martin

+0

添加了一個小例子。很明顯,1e-6值的順序存在很小的差異。這可能是因爲這些小的變化會被放大,並在幾次矩陣乘法和積分之後產生實質性的錯誤。我仍然試圖準確地確定使用atan2/fmodf時可能導致問題的原因。 – Martin

+0

好吧,我已經增強了我的答案。仍然無法找出根本原因,但試圖解釋和證明你的設計選擇。 –

3

floatdouble數學。

當然,第三種方法效果最好。它使用double數學。


看看b1, b3b3肯定是由於fmodf()調用而導致精度爲float

請注意M_PI通常是double,所以b1 -= 2.f * M_PI;很可能與double精度數學完成,並提供更準確的答案。 f2.f不強制產品2.f * M_PIfloat - 產品是double,因此是-=

b1 -= 2.f * M_PI; 
// same as 
b1 = (float)((double)b1 - (2.f * M_PI)); 

此外:與優化和FLT_EVAL_METHOD > 0,C被允許在比類型的精度更高執行FP代碼。 b1可能會在double處計算,即使代碼出現float。憑藉更高的精度,而事實上,M_PI(有理數)是不完全π(無理數),導致更準確b1fmodf(a3, 2 * M_PI);

float b1 = a1; 
if(b1 > 2.f * M_PI) b1 -= 2.f * M_PI; // double math 
if(b1 < 0.f) b1 += 2.f * M_PI;   // double math 
float b3 = fmodf(a3, 2 * M_PI); 

爲了確保float結果,用volatile float b1 = a1;做到公平比較和使用float常數如#define M_PIf ((float) M_PI)

此外。通過公平的比較,最好使用if(b1 < -2.f * M_PIf) b1 += 2.f * M_PIf;

推薦OP print FLT_EVAL_METHOD以幫助進一步討論。

#include <float.h> 
printf("%d\n", FLT_EVAL_METHOD); 

OP有2個解決方案:

  1. 使用更廣泛的數學像double對於敏感的弧度減少。

    float b3 = fmod(a3, 2 * M_PI); // not fmodf 
    
  2. 不要使用弧度,但角度測量像度或BAM並執行準確的範圍減少。在trig調用之前,角度需要度數來進行弧度轉換。

    float b3 = fmodf(a3, 360.0f); // use fmodf, a3, b3 are in degrees 
    

注:float b2 = atan2f(sinf(a2), cosf(a2));方法是不是一個合理的競爭者。

+0

關於你的最後一點,我曾建議在註釋中考慮使用'sinpi()'和'cospi()',這意味着將角度視爲π的分數,並且在非常大範圍, – njuffa

+0

@njuffa謝謝。 'sinpi(),cospi()'很有用,但不幸的是它們還不是標準C庫的一部分。 – chux

+0

這就是爲什麼我發佈合理的C實現,以紀念今年的π日[這裏](https://stackoverflow.com/questions/42792939/implementation-of-sinpi-and-cospi-using-standard-c-數學庫)。 – njuffa