Evitare l'aritmetica in virgola mobile

Question

Ho scritto un piccolo sintetizzatore software per iPhone.
Per ottimizzare ulteriormente le prestazioni ho misurato la mia applicazione con Shark e ho scoperto che sto perdendo molto tempo nelle conversioni float / SInt16.
Quindi ho riscritto alcune parti per aggirare le conversioni pre-calcolando le tabelle di ricerca che restituiscono " pronto per l'uso " SInt16 campioni. Finora funziona bene.
Attualmente sto cercando di riscrivere alcuni filtri e l'implementazione della mia busta ADSR per usare solo l'aritmetica intera ma potrei usare alcuni consigli su come eseguire moltiplicazioni / divisioni senza float.
Sto prendendo di mira formato canonico iPhone :

• LPCM
• campioni interi a 16 bit

Quali sono i buoni approcci per applicare un'ampiezza al mio campione finale senza usare un float?

Modifica
L'unica cosa che ho capito finora è che posso dividere per poteri di 2 spostando a destra il mio campione attuale.

inBuffer[frame] = wavetable[i % cycleLengthInSamples] >> 4;

Ma non riesco a pensare a un modo elegante per creare una busta ADSR liscia con quello.

Edit2: Grazie per tutte le tue fantastiche risposte!
Il mio approccio attuale:

• porta tutti i miei valori di busta ADSR nella gamma SInt16 positiva
• moltiplica per il valore corrente della wavetable (memorizza gli intermedi come SInt32)
• sposta il risultato di 16 a destra

sembra funzionare :)

Answer 1

Il punto fisso è buono, poiché in questo caso stai usando 16 bit. Il modo più semplice è moltiplicare per una potenza di 10 a seconda della precisione di cui hai bisogno. Se puoi utilizzare ints a 32 bit come intermedio, dovresti essere in grado di ottenere una precisione decente. Alla fine puoi tornare a un int a 16 bit, arrotondando o troncando come preferisci.

Modifica: Volete spostare a sinistra, per aumentare i valori. Memorizza il risultato del turno in un tipo con maggiore precisione (32 o 64 bit a seconda delle necessità). il semplice spostamento non funzionerà se si utilizzano tipi firmati

Fai attenzione se stai moltiplicando o dividendo due numeri a virgola fissa. Il moltiplicarsi finisce per essere (a * n) * (b n) e finirai con un b n ^ 2 invece di un b n. La divisione è (a n) / (b n) che è (a / b) anziché ((a n) / b). Ecco perché ho suggerito di usare potenze di 10, rende facile trovare i tuoi errori se non hai familiarità con il punto fisso.

Quando hai finito i tuoi calcoli, torni indietro a destra per tornare a un int 16 bit. Se vuoi divertirti, puoi anche fare un arrotondamento prima di cambiare.

Ti suggerisco di fare qualche lettura se sei veramente interessato a implementare un punto fisso efficiente. http://www.digitalsignallabs.com/fp.pdf

Answer 2

Le risposte a questa domanda SO sono piuttosto completo in termini di implementazione. Ecco un po 'più di una spiegazione di quella che ho visto laggiù:

Un approccio è forzare tutti i tuoi numeri in un intervallo, diciamo [-1.0,1.0). Quindi si mappano quei numeri nell'intervallo [-2 ^ 15, (2 ^ 15) -1]. Ad esempio,

Half = round(0.5*32768); //16384
Third = round((1.0/3.0)*32768); //10923

Quando moltiplichi questi due numeri ottieni

Temp = Half*Third; //178962432
Result = Temp/32768; //5461 = round(1.0/6.0)*32768

La divisione per 32768 nell'ultima riga è il punto Patros si moltiplica necessitando di un ulteriore passaggio in scala. Questo ha più senso se scrivi esplicitamente il ridimensionamento 2 ^ N:

x1 = x1Float*(2^15);
x2 = x2Float*(2^15);
Temp = x1Float*x2Float*(2^15)*(2^15);
Result = Temp/(2^15); //get back to 2^N scaling

Quindi questa è l'aritmetica. Per l'implementazione, nota che la moltiplicazione di due numeri interi a 16 bit richiede un risultato a 32 bit, quindi Temp dovrebbe essere a 32 bit. Inoltre, 32768 non è rappresentabile in una variabile a 16 bit, quindi tieni presente che il compilatore renderà immediati a 32 bit. E come hai già notato, puoi passare a moltiplicare / dividere per poteri di 2 in modo da poter scrivere

N = 15;
SInt16 x1 = round(x1Float * (1 << N));
SInt16 x2 = round(x2Float * (1 << N));
SInt32 Temp = x1*x2;
Result = (SInt16)(Temp >> N);
FloatResult = ((double)Result)/(1 << N);

Ma supponi che [-1,1) non sia il giusto intervallo? Se preferisci limitare i tuoi numeri a, diciamo, [-4.0,4.0), puoi usare N = 13. Quindi hai 1 bit di segno, due bit prima del punto binario e 13 dopo. Questi sono chiamati rispettivamente 1,15 e 3,13 tipi frazionari a punto fisso. Scambia precisione nella frazione per headroom.

L'aggiunta e la sottrazione di tipi frazionari funziona bene finché si cerca la saturazione. Per dividere, come diceva Patros, il ridimensionamento si annulla effettivamente. Quindi devi farlo

Quotient = (x1/x2) << N;

o, per preservare la precisione

Quotient = (SInt16)(((SInt32)x1 << N)/x2); //x1 << N needs wide storage

Moltiplicare e dividere per numeri interi funziona normalmente. Ad esempio, per dividere per 6 puoi semplicemente scrivere

Quotient = x1/6; //equivalent to x1Float*(2^15)/6, stays scaled

E in caso di divisione per un potere di 2,

Quotient = x1 >> 3; //divides by 8, can't do x1 << -3 as Patros pointed out

L'aggiunta e la sottrazione di numeri interi, tuttavia, non funziona in modo ingenuo. Devi prima vedere se l'intero si adatta al tuo tipo x.y, crea il tipo frazionario equivalente e procedi.

Spero che questo aiuti con l'idea, guarda il codice in quell'altra domanda per implementazioni pulite.

Answer 3

Dai un'occhiata a questa pagina che descrive gli algoritmi di moltiplicazione rapida.

http://www.newton.dep.anl.gov /askasci/math99/math99199.htm

Answer 4

In generale, supponiamo che utilizzerai una rappresentazione a punto fisso 16.16 firmata. In modo che un numero intero a 32 bit avrà una parte intera a 16 bit con segno e una parte frazionaria a 16 bit. Quindi non so quale lingua viene utilizzata nello sviluppo di iPhone (forse Objective-C?), Ma questo esempio è in C:

#include <stdint.h>

typedef fixed16q16_t int32_t ;
#define FIXED16Q16_SCALE 1 << 16 ;

fixed16q16_t mult16q16( fixed16q16_t a, fixed16q16_t b )
{
    return (a * b) / FIXED16Q16_SCALE ;
}

fixed16q16_t div16q16( fixed16q16_t a, fixed16q16_t b )
{
    return (a * FIXED16Q16_SCALE) / b ;
}

Si noti che quanto sopra è un'implementazione semplicistica e non fornisce alcuna protezione dall'overflow aritmetico. Ad esempio in div16q16 () I multiplo prima della divisione per mantenere la precisione, ma a seconda degli operandi l'operazione potrebbe traboccare. È possibile utilizzare un intermedio a 64 bit per ovviare a questo. Inoltre, la divisione viene sempre arrotondata perché utilizza la divisione intera. Ciò offre le migliori prestazioni, ma può influire sulla precisione dei calcoli iterativi. Le correzioni sono semplici ma aggiungono al sovraccarico.

Nota che quando si moltiplica o si divide per una potenza costante di due, la maggior parte dei compilatori individuerà la banale ottimizzazione e utilizzerà un turno. Comunque C non definisce il comportamento per uno spostamento a destra di un intero con segno negativo, quindi l'ho lasciato al compilatore per risolverlo per sicurezza e portabilità. YMV in qualunque lingua tu stia utilizzando.

In un linguaggio OO, fixed16q16_t sarebbe naturalmente un candidato per una classe con sovraccarico dell'operatore in modo da poterlo usare come un normale tipo aritmetico.

Potresti trovare utile la conversione tra tipi:

double fixed16q16_to_double( fixed16q16_t fix )
{
    return (double)fix / FIXED16Q16_SCALE ;
}

int fixed16q16_to_int( fixed16q16_t fix )
{
    // Note this rounds to nearest rather than truncates
    return ((fix + FIXED16Q16_SCALE/2)) / FIXED16Q16_SCALE ;
}

fixed16q16_t int_to_fixed16q16( int i )
{
    return i * FIXED16Q16_SCALE ;
}

fixed16q16_t double_to_fixed16q16( double d )
{
    return (int)(d * FIXED16Q16_SCALE) ;
}

Queste sono le basi, è possibile diventare più sofisticati e aggiungere trig e altre funzioni matematiche.

Aggiunte e sottrazioni fisse funzionano con gli operatori + e - integrati e le loro varianti.