Livres / ressources pour implémenter diverses fonctions mathématiques en arithmétique à virgule fixe à des fins DSP

8

Je recherche des livres ou des ressources qui couvrent en détail les éléments suivants:

  • mettre en œuvre des fonctions mathématiques (par exemple, logarithme, exponentielle, sinus, cosinus, inverse) en arithmétique à virgule fixe à des fins DSP.

  • des techniques telles que l'utilisation de tables de recherche, de séries Taylor, etc.

Je suis assez familier avec la programmation C et je suis plus intéressé par les algorithmes sur la façon de mettre en œuvre efficacement diverses fonctions mathématiques.

RuD
la source
1
Ce n'est qu'une astuce, mais c'est très utile. Il s'agit de calculer la fonction atan2, c'est-à-dire de calculer l'argument d'un nombre complexe à partir de ses parties réelles et imaginaires.
Matt L.
1
je peux vous donner quelques séries optimisées d'ordre fini que j'ai développées il y a plus de dix ans. à part le client initial, je n'ai pas obtenu d'argent pour cela car je pense que je peux aussi bien le rendre du domaine public. à l'origine, il a été développé pour virgule flottante dans un contexte où le client ne pouvait pas inclure le stdlib dans la construction. et l'évaluation n'est pas du tout parfaite. c'est-à-dire qu'il y a une erreur, mais très petite et optimisée pour la fonction particulière évaluée. lemme trouver ce fichier, je vais retirer les coefficients et poster les séries.
robert bristow-johnson
@ robertbristow-johnson j'ai hâte de profiter de la série et de voir comment ça se passe! Merci!
RuD
Ruchir, j'ai posté les séries ci-dessous. il y a des choses de bon sens que vous devez faire pour étendre la portée. comme
2x=2k×2xk
si kxk+1. chose similaire pourlog2()et les sinusoïdes périodiques. pour exp et log, cela nécessitera un décalage arithmétique des bits de ce qui sort (pour exp) ou de ce qui entre (pour log). je pense que vous pouvez comprendre cela. et, bien sûr, pour exp et log de différentes bases (commee), vous mettez à l'échelle avec la constante appropriée ce qui se passe dans 2x et ce qui sort log2(x).
robert bristow-johnson

Réponses:

9

la forme polynomiale générale est:

f(u)=n=0N an un=a0+(a1+(a2+(a3+...(aN2+(aN1+aNu)u)u ...)u)u)u

cette dernière forme utilise la méthode de Horner , ce qui est fortement recommandé, surtout si vous le faites en virgule flottante simple précision.

puis pour quelques fonctions spécifiques:

racine carrée:

f(x1)x1x2N=4une0=1.0une1=0,49959804148061une2=-0.12047308243453une3=0,04585425015501une4=-0,01076564682800

si 2X4, utilisez ce qui précède pour évaluer X2 et multiplier ce résultat avec 2 obtenir x. comme aveclog2(x), appliquez la puissance de 2 mise à l'échelle pour mettre l'argument à l'échelle dans la plage nécessaire.

logarithme en base 2:

xf(x1)log2(x)1x2N=5a0=1.44254494359510a1=0.7181452567504a2=0.45754919692582a3=0.27790534462866a4=0.121797910687826a5=0.02584144982967

exponentielle base 2:

f(x)2x0x1N=4a0=1.0a1=0.69303212081966a2=0.24137976293709a3=0.05203236900844a4=0.01355574723481

sinus:

XF(X2)péché(π2X)-1X1N=4une0=1.57079632679490une1=-0,64596406188166une2=0,07969158490912une3=-0,00467687997706une4=0,00015303015470

cosinus (utiliser sinus):

cos(πX)=1-2péché2(π2X)

tangente:

bronzer(X)=péché(X)cos(X)

tangente inverse:

XF(X2)arctan(X)-1X1N=4une0=1.0une1=0,33288950512027une2=-0,08467922817644une3=0,03252232640125une4=-0,00749305860992

arctan(X)=π2-arctan(1X)1X

arctan(X)=-π2-arctan(1X)X-1

sinus inverse:

arcsin(X)=arctan(X1-X2)

cosinus inverse:

arccos(X)=π2-arcsin(X)=π2-arctan(X1-X2)
robert bristow-johnson
la source
Cela semble assez utile! Merci beaucoup de partager cela. Mais la première entrée de référence man soxest toujours la meilleure;)
jojek
ne sais pas sox. qu'en dit le manuel?
robert bristow-johnson
2
Tout simplement [1] R. Bristow-Johnson, Cookbook formulae for audio EQ biquad filter coefficients, http://musicdsp.org/files/Audio-EQ-Cookbook.txt:)
jojek
BTW, la série minimise l' erreur pondérée maximale . l'erreur est pondérée d'une manière qui a du sens pour la fonction. erreur maximale uniforme pourJournal2(). erreur maximale proportionnelle pourX et 2X. quelque chose de similaire à proportionnel pourpéché() ayant à voir avec le calcul des coefficients pour les filtres résonants de sorte que l'erreur maximale de Journal(F0)est minimisé. je ne me souviens pas pour quel critère j'ai utiliséarctan(). et pour une raison quelconque, je ne trouve pas mon fichier qui me dit quelles sont les erreurs maximales, étant donné la plage deX. quelqu'un avec MATLAB peut le découvrir.
robert bristow-johnson
1
vous devriez tracer l'erreur avec votre python, si vous le pouvez. np.max(np.abs(sqrt_1px(xp)-np.sqrt(1+xp)))pourrait aussi être à la place np.max(np.abs((sqrt_1px(xp)-np.sqrt(1+xp))/np.sqrt(1+xp))) et même pour 2**x la pondération d'erreur pour le péché est différent et je vais devoir trouver comment je l'ai fait. j'ai d'anciens scripts MATLAB qui fonctionnaient auparavant dans Octave, mais maintenant je ne peux même plus faire tracer Octave sur mon ancien ordinateur portable Mac G4.
robert bristow-johnson
2

Bien qu'il ne soit pas spécifique au point fixe, je recommande fortement le livre "Math Toolkit for Real-Time Programming" de Jack Crenshaw. Il est livré avec un CD avec le code source.

David
la source
2

TI possède des bibliothèques IQMath pour tous ses microcontrôleurs à virgule fixe. Je les ai trouvés être une mine d'or de fonctions mathématiques et DSP à virgule fixe qui ne se limitent pas nécessairement aux puces TI.

MSP430 C28X

Otto Hunt
la source
Je suis plus intéressé par les algorithmes que par la simple implémentation des fonctions
RuD
1

L'approximation de Chebyshev peut aider à calculer des coefficients polynomiaux qui sont proches de l'optimum pour approximer une fonction sur une plage finie. Vous exécutez la routine d'approximation sur un PC pour obtenir un ensemble particulier de coefficients polynomiaux, que vous pouvez ensuite appliquer sur la plate-forme que vous souhaitez (par exemple intégrée / DSP). Les petits caractères sont plus ou moins les suivants:

  • Cela ne fonctionne que pour les fonctions d'une variable; si vous avez une fonctionz=F(X,y) alors l'approximation de Chebyshev ne va pas vous aider.
  • La fonction que vous approximez doit être de type "polynomial". Les coins, les virages serrés et de nombreux mouvements nécessitent des polynômes d'ordre supérieur pour atteindre un niveau de précision donné.
  • Il est important de maintenir l'ordre polynomial bas - au-dessus du 5e environ, vous pouvez commencer à voir des erreurs numériques.
  • L'approximation de Chebyshev utilise les polynômes de Chebyshev évalués sur le domaine [-1,1], donc si la plage d'entrée de votre fonction est significativement différente, vous devrez peut-être mettre à l'échelle / compenser l'entrée de manière appropriée avant de déterminer les coefficients et avant de les appliquer. Par exemple, si je me soucie d'une fonction sur la plage d'entréeX[0,20] alors je pourrais définir u=(X×0,1)-1 pour que u varie de -1 à +1 et je peux évaluer un polynôme dans upour calculer le résultat requis. Sans cette mise à l'échelle, vous pouvez rencontrer des erreurs numériques plus facilement - ou déclaré d'une autre manière, pour la même précision, vous devrez peut-être des longueurs de bits plus élevées pour calculer les valeurs intermédiaires, ce qui est généralement indésirable.
Jason S
la source
Jason, je n'ai pas utilisé de polynômes de Tchebyshev, mais j'ai utilisé une erreur MinMax pondérée (parfois appelée "Lnorm "on error), qui est, je pensais, la même que l'approximation de Tchebyshev. la méthode que j'ai utilisée était l'algorithme d'échange Remez.
robert bristow-johnson
Remez est supérieur (quoique plus complexe) à Chebyshev. Chebyshev se rapproche juste de la condition minimax, mais généralement c'est assez bon.
Jason S