Come calcolare un'onda sinusoidale con precisione nel tempo

Question

Caso d'uso è quello di generare un'onda sinusoidale per la sintesi digitale, in modo, abbiamo bisogno di calcolare tutti i valori del peccato (dt) dove:

t è un numero intero numero, che rappresenta il numero del campione. Questo è variabile. La gamma va da 0 a 158.760.000 per un suono di un'ora di qualità CD.

d è doppio, rappresenta il delta dell'angolo. Questo è costante. E l'intervallo è: maggiore di 0, minore di pi.

obiettivo è di ottenere un'elevata precisione con tradizionali int e doppie tipi di dati. Le prestazioni non sono importanti.

Naive implementazione è:

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

Ma, il problema è quando t aumenta, la precisione si riduce a causa grandi numeri forniti alla funzione "peccato".

Una versione migliorata è la seguente:

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

    return sin(d_sum); 
} 

Ecco, faccio in modo di fornire i numeri nella gamma da 0 a 2 * pi alla funzione di "peccato".

Ma, ora, il problema è quando d è piccolo, ci sono molte piccole aggiunte che diminuiscono la precisione ogni volta.

La domanda qui è come migliorare la precisione.

---

Appendice 1

"precisione viene ridotto perché i grandi numeri forniti a "" funzione" peccato:

#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; 
} 

uscita:

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

Answer 1

È possibile provare un approccio che viene utilizzato è alcune implementazioni di trasformata veloce di Fourier n. I valori della funzione trigonometrica sono calcolati in base a valori e delta precedenti.

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

Qui abbiamo memorizzare e aggiornare il valore del coseno troppo e memorizzare costante (per proposta delta) fattori Cos (d) e Sin (d).

Ora circa la precisione: coseno (d) per piccola d è molto vicino a 1, quindi c'è il rischio di perdita di precisione (ci sono solo poche cifre significative in numeri come 0.99999987).Per superare questo problema, è possibile memorizzare costanti come

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

utilizzando altre formule per aggiornare il valore corrente
(qui sa = Sin(A) per il valore corrente, ca = Cos(A) per il valore corrente)

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

P.S. Alcune implementazioni FFT periodicamente (ogni passi K) rinnovano i valori sa e ca tramite trig. funzioni per evitare l'accumulo di errori.

Esempio di risultato. Calcoli in doppio.

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 ∙ π), in modo che il problema può essere riassunta in trovare accuratamente un piccolo numero che è uguale ad un gran numero x modulo 2π.

Per esempio, -1,61059759 ≅ 256 mod 2π, e si può calcolare sin(-1.61059759) con più precisione rispetto sin(256)

Quindi cerchiamo di scegliere un numero intero con cui lavorare, 256. In primo luogo trovare piccoli numeri, che sono pari a poteri di 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) 

e poi dividere l'indice come un numero in base di 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; 

ora è possibile trovare un numero molto inferiore x che è pari a 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); 

per diversi valori di d dovrete calcolare i valori diversi modB, modC e modD, che sono uguali a quelle di 256, ma modulo (2π/d). Potresti usare una libreria di alta precisione per questi due calcoli.

Answer 3

Scala fino al periodo di 2^64, e fare la moltiplicazione utilizzando aritmetica intera:

// 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); 

Finché il periodo non è miliardi di campioni, la precisione di multiplier sarà abbastanza buono.

Answer 4

Il seguente codice mantiene l'input della funzione sin() all'interno di un intervallo ridotto, riducendo al tempo stesso il numero di piccole aggiunte o sottrazioni dovute a un incremento di fase potenzialmente molto piccolo.

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

Per la precisione iper, OP ha 2 problemi:

1. moltiplicando d da n e mantenere più precisione rispetto double.Questo è stato risposto nella prima parte qui sotto.

2. Esecuzione di un mod del periodo. La soluzione semplice è usare i gradi e poi mod 360, abbastanza facile da fare esattamente. Per fare 2*π di grandi angoli è difficile in quanto necessita di un valore di 2*π con circa 27 più bit di precisione rispetto (double) 2.0 * M_PI

---

Usa 2 double s per rappresentare d.

Ipotizziamo 32 bit int e binary64double. Quindi double ha una precisione di 53 bit.

0 <= n <= 158,760,000 che è di circa 2 ^27,2. Poiché double è in grado di gestire numeri interi senza segno a 53 bit in modo continuo ed esatto, 53-28 -> 25, qualsiasi double con solo 25 bit significativi può essere moltiplicato per n ed essere ancora esatto.

segmento d in 2 double s dmsb,dlsb, le 25 cifre più significative e il 28- minimo.

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 

Poi ogni moltiplicazione (o successiva aggiunta) di dmsb*n saranno esatti. (questa è la parte importante.) dlsb*n errore solo nei suoi minimi alcuni bit.

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); 
} 

Nota: fmod(x,y) risultati sono attesi per l'esattezza dare esatto 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); 
} 

uscita

0.04630942695385031893 

---

Utilizzare i gradi

consiglia di utilizzare gradi come 360 gradi è il esatte radianti epoca e M_PI*2 è un app roximation. C non può rappresentare π esattamente.

Se OP ha ancora voglia di usare radianti, per una visione più su come eseguire il mod di π, vedere Good to the Last Bit