jak obliczyć sinusoidę o dokładności pomiaru czasu

Question

przypadków użycia jest generowanie fali sinusoidalnej syntezy cyfrowego tak, trzeba obliczyć wszystkie wartości sin (dt) gdzie:

t jest liczbą całkowitą numer, reprezentujący numer próbki. To jest zmienne. Zakres wynosi od 0 do 158 760 000 na 1 godzinę dźwięku o jakości CD.

d jest podwójny, reprezentujący delta kąta. To jest stałe. Zakres wynosi: większy niż 0, mniejszy niż pi.

Celem jest osiągnięcie wysokiej dokładności z tradycyjnymi typami danych int i podwójnymi. Wydajność nie jest ważna.

naiwna implementacja jest:

double next() 
{ 
    t++; 
    return sin(((double) t) * (d)); 
} 

Ale problem jest, gdy t zwiększa dokładność zostanie zmniejszona, ponieważ duże numery dostarczane do funkcji "grzechu".

Ulepszona wersja jest następująca:

double next() 
{ 
    d_sum += d; 
    if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2); 

    return sin(d_sum); 
} 

Tutaj I upewnij się, aby zapewnić numery w zakresie od 0 do 2 * pi do funkcji „grzechu”.

Ale teraz jest problem, gdy d jest mały, istnieje wiele małych dodatków, które zmniejsza dokładność za każdym razem.

Pytanie tutaj brzmi: jak poprawić dokładność.

---

Dodatek 1

"dokładność zostanie zmniejszona, ponieważ duże cyfry przewidziane do "" funkcji" grzech

#include <stdio.h> 
#include <math.h> 

#define TEST  (300000006.7846112) 
#define TEST_MOD (0.0463259891528704262050786960234519968548937998410258872449766) 
#define SIN_TEST (0.0463094209176730795999323058165987662490610492247070175523420) 

int main() 
{ 
    double a = sin(TEST); 
    double b = sin(TEST_MOD); 

    printf("a=%0.20f \n" , a); 
    printf("diff=%0.20f \n" , a - SIN_TEST); 
    printf("b=%0.20f \n" , b); 
    printf("diff=%0.20f \n" , b - SIN_TEST); 
    return 0; 
} 

wyjściowa:

a=0.04630944601888796475 
diff=0.00000002510121488442 
b=0.04630942091767308033 
diff=0.00000000000000000000 

Answer 1

Można spróbować podejście, które jest używane, to niektóre implementacje szybkiej transformacji Fouriera n. Wartości funkcji trygonometrycznych są obliczane na podstawie poprzednich wartości i wartości delta.

Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d) 

Mamy tu do przechowywania i aktualizacja wartości cosinus zbyt i przechowywania stałe (dla danego delta) współczynników cos (d) i sin (d).

Teraz o precyzji: cosinus (d) dla małych d jest bardzo zbliżony do 1, więc istnieje ryzyko utraty precyzji (jest tylko kilka cyfr znaczących w liczbach takich jak 0,99999987).Aby rozwiązać ten problem, możemy przechowywać stale czynniki jak

dc = Cos(d) - 1 = - 2 * Sin(d/2)^2 
ds = Sin(d) 

używając innej formuły, aby zaktualizować wartość prądu
(tutaj sa = Sin(A) do wartości bieżącej, ca = Cos(A) dla aktualnej wartości)

ts = sa //remember last values 
tc = ca 
sa = sa * dc + ca * ds 
ca = ca * dc - ts * ds 
sa = sa + ts 
ca = ca + tc 

PS: Niektóre implementacje FFT okresowo (co K kroków) odnawiają wartości sa i poprzez trig. funkcje w celu uniknięcia kumulacji błędów.

Przykładowy wynik. Obliczenia w podwójnych.

d=0.000125 
800000000 iterations 
finish angle 100000 radians 

          cos    sin 
described method  -0.99936080743598 0.03574879796994 
Cos,Sin(100000)   -0.99936080743821 0.03574879797202 
windows Calc   -0.9993608074382124518911354141448 
          0.03574879797201650931647050069581   

Answer 2

sin (x) = sin (x + 2N ∙ π), więc problem można sprowadzić do dokładnego znalezienie małej ilości, która jest równa w dużej liczbie x modulo 2π.

Na przykład -1,61059759 ≅ 256 mod 2π, można obliczyć sin(-1.61059759) z większą precyzją niż sin(256)

Warto więc wybrać pewną liczbę całkowitą do pracy, 256. Po pierwsze znaleźć małe numery, które są równe uprawnienia 256, modulo 2π:

// to be calculated once for a given frequency 
// approximate hard-coded numbers for d = 1 below: 
double modB = -1.61059759; // = 256 mod (2π/d) 
double modC = 2.37724612; // = 256² mod (2π/d) 
double modD = -0.89396887; // = 256³ mod (2π/d) 

a następnie podzielić indeks jako numer w bazie 256:

// split into a base 256 representation 
int a = i   & 0xff; 
int b = (i >> 8) & 0xff; 
int c = (i >> 16) & 0xff; 
int d = (i >> 24) & 0xff; 

Obecnie można znaleźć wiele mniejszą liczbę x która jest równa i modulo 2π/d

// use our smaller constants instead of the powers of 256 
double x = a + modB * b + modC * c + modD * d; 
double the_answer = sin(d * x); 

dla różnych wartości d będziesz musiał obliczyć różne wartości modB, modC i modD, które są równe tym mocom 256, ale modulo (2π/d). Przy tych kilku obliczeniach można użyć biblioteki o wysokiej precyzji.

Answer 3

Scale się ten okres do 2^64, a nie mnożenie za pomocą arytmetyki liczb całkowitych:

// constants: 
double uint64Max = pow(2.0, 64.0); 
double sinFactor = 2 * M_PI/(uint64Max); 

// scale the period of the waveform up to 2^64 
uint64_t multiplier = (uint64_t) floor(0.5 + uint64Max * d/(2.0 * M_PI)); 

// multiplication with index (implicitly modulo 2^64) 
uint64_t x = i * multiplier; 

// scale 2^64 down to 2π 
double value = sin((double)x * sinFactor); 

Dopóki okres nie jest miliardy próbek, precyzja multiplier będzie wystarczająco dobry.

Answer 4

Poniższy kod zachowuje dane wejściowe funkcji sin() w niewielkim zakresie, jednocześnie nieco zmniejszając liczbę małych dodatków lub odejmowań z powodu potencjalnie bardzo małego przyrostu fazy.

double next() { 
    t0 += 1.0; 
    d_sum = t0 * d; 
    if (d_sum > 2.0 * M_PI) { 
     t0 -= ((2.0 * M_PI)/d); 
    } 
    return (sin(d_sum)); 
} 

Answer 5

hiper Dla dokładności, PO ma 2 problemy:

1. namnażania d przez n i utrzymanie większej precyzji niż double.Odpowiedzi udzielono w pierwszej części poniżej.

2. Wykonywanie mod okresu. Prostym rozwiązaniem jest użycie stopni, a następnie mod 360, łatwe do zrobienia dokładnie. W tym celu 2*π z dużymi kątami jest trudne, ponieważ wymaga to wartość 2*π z około 27 bitów, w dokładności niż (double) 2.0 * M_PI

---

użyciu 2 double s reprezentuje d.

Załóżmy 32-bitowe int i binary64double. Tak więc double ma 53-bitową dokładność.

co stanowi około 2 ^27,2. Ponieważ double może obsługiwać 53-bitowe liczby całkowite bez znaku w sposób ciągły i dokładnie, 53-28 -> 25, dowolne double z tylko 25 znaczącymi bitami może być pomnożone przez n i nadal być dokładne.

Segment d na 2 double s dmsb,dlsb, 25 najbardziej znaczących cyfr i co najmniej 28.

int exp; 
double dmsb = frexp(d, &exp); // exact result 
dmsb = floor(dmsb * POW2_25); // exact result 
dmsb /= POW2_25;    // exact result 
dmsb *= pow(2, exp);   // exact result 
double dlsb = d - dmsb;  // exact result 

Następnie każdy mnożenia (lub kolejne dodawanie) dmsb*n będzie dokładne. (to jest ważna część.) dlsb*n będzie tylko błąd w najmniejszej liczbie bitów.

double next() 
{ 
    d_sum_msb += dmsb; // exact 
    d_sum_lsb += dlsb; 
    double angle = fmod(d_sum_msb, M_PI*2); // exact 
    angle += fmod(d_sum_lsb, M_PI*2); 
    return sin(angle); 
} 

Uwaga: Przewiduje się fmod(x,y) Wyniki należy podać dokładny dokładny x,y.

---

#include <stdio.h> 
#include <math.h> 

#define AS_n 158760000 
double AS_d = 300000006.7846112/AS_n; 
double AS_d_sum_msb = 0.0; 
double AS_d_sum_lsb = 0.0; 
double AS_dmsb = 0.0; 
double AS_dlsb = 0.0; 

double next() { 
    AS_d_sum_msb += AS_dmsb; // exact 
    AS_d_sum_lsb += AS_dlsb; 
    double angle = fmod(AS_d_sum_msb, M_PI * 2); // exact 
    angle += fmod(AS_d_sum_lsb, M_PI * 2); 
    return sin(angle); 
} 

#define POW2_25 (1U << 25) 

int main(void) { 
    int exp; 
    AS_dmsb = frexp(AS_d, &exp);   // exact result 
    AS_dmsb = floor(AS_dmsb * POW2_25); // exact result 
    AS_dmsb /= POW2_25;     // exact result 
    AS_dmsb *= pow(2, exp);    // exact result 
    AS_dlsb = AS_d - AS_dmsb;   // exact result 

    double y; 
    for (long i = 0; i < AS_n; i++) 
    y = next(); 
    printf("%.20f\n", y); 
} 

wyjścia

0.04630942695385031893 

---

Użyj stopni

Poleca używać stopni jak 360 stopni jest dokładnych okresie i M_PI*2 radianów jest aplikacja roximation. C nie może reprezentować π dokładnie.

Jeśli PO nadal chce korzystać radianach do dalszego wglądu na wykonywanie mod Õ zobacz Good to the Last Bit