ˇ aritmetické operace Složitejší ˇ SCS Ing. Jakub Št’astný, Ph.D.1 1
[email protected] FPGA Laboratoˇr/Laboratoˇr zpracování biologických signálu˚ ˇ Katedra teorie obvodu˚ FEL CVUT Technická 2, Praha 6, 166 27 http://amber.feld.cvut.cz/fpga
ˇ 21. kvetna 2014
Varování
Pozor! Tato verze pˇrednášky obsahuje materiál tˇretích stran. Je urˇcena pouze k internímu použití a vystavení na WWW stránkách Laboratoˇre biologických signálu, ˚ Katedry teorie ˇ obvodu˚ CVUT FEL, Praha!
ˇ funkce 1 Složitejší Druhá odmocnina Logaritmus Goniometrické funkce
Motivace
• druhá odmocnina – vedeckotechnické ˇ výpoˇcty (nemáme to
ˇ :-) ), multimediální a grafické výpoˇcty rádi, ale lidé to chtejí (velikost vektoru)
• logaritmus – pˇrevádí složité na jednodušší operace • harmonické funkce – všude kolem nás
Druhá odmocnina
Výpoˇcet druhé odmocniny – triviální metoda
√
• y = x • jedoduchý algoritmus 1 y =0 √ 2 if y × y > x , konec a ⌊ x⌋ = y − 1 3 y = y + 1, goto 2 • nevhodné – póóóómalééééé
Výpoˇcet druhé odmocniny – triviální metoda 2 • y =
√
x
• pozorování: ∆ = (x + 1)2 − x 2 = x 2 + 2x + 1 − x 2 = 2x + 1 • pˇríklad: 1 = 12 ,1 + 3 = 22 ,1 + 3 + 5 = 32 , ...
PN−1
+ 1) = N 2 √ • odhad hodnoty x •
x=0 (2x
1 2 3
y = 0, acc = 0 √ if acc > x , konec a ⌊ x ⌋ = y − 1 acc = acc + 2y + 1, y = y + 1, goto 2
• stále póóómaléééé, ale už bez násobení
• pro odmocninu z 32 bitového cˇ ísla potˇrebuji až 65536
kroku! ˚
Výpoˇcet druhé odmocniny – triviální metoda 3 • muže ˚ mi pomoci rozklad po ˇrádech √ • pˇríklad: druhá odmocnina z 32 bit. cˇ ísla, y = x • krok 1 √ x 1 xMSB = 65536 , urˇcím yMSB = xMSB 2 xMSB < 65536, takže maximáln √ eˇ 256 kroku˚ 3 yLSB = yMSB ∗ 256, protože 65536 = 256 • krok 2 – yLSB použiji jako startovací hodnotu pro iterativní
• • • •
algoritmus, v dalších max. 256 krocích najdu pˇresnou hodnotu celkem cca 512 kroku˚ max. místo 65536! lze urychlit rozkladem na více základu˚ (historie – ENIAC – implementace odmocniny) stále pomalé, ale už lepší pro celá cˇ ísla, ale desetinná cˇ ísla mohu škálovat (×22N a potom /2N )
Babylónská metoda
• poprvé zdokumentována v práci Herona Alexandrijského
(≈10-70 n.l.)
• zaˇcínáme s odhadem y0 • opakujeme postup dokud se yi “významneˇ mení" ˇ x 1 1 yi+1 = 2 (yi + y ) i • jednoduchá myšlenka – druhá odmocnina musí být nekde ˇ
mezi aktuálním odhadem a
x yi
• poˇcet platných cˇ íslic se zdvojnásobuje v každé iterace
(dukaz ˚ za domácí úkol)
Technika bisekce
• triviální postup – zaˇcínáme s odhadem intervalu napˇr.
< y0D , y0H >=< 0, x > • opakujeme postup dokud se velikost yiH − yiD “významneˇ ˇ mení" 1 2
pokud je ( 21 (yiD + yiH ))2 > x , použij jako další odhad intervalu < yi+1D , yi+1H >=< yiD , ( 12 (yiD + yiH )) > pokud je ( 21 (yiD + yiH ))2 <= x , použij jako další odhad intervalu < yi+1D , yi+1H >=< ( 12 (yiD + yiH )), yiH >
• poˇcet platných dvojkových smíst roste o jedno v každé
iteraci
• pomalý výpoˇcet...
Krátké opakování Newtonovy metody
• ˇrešíme rovnici f (x) = 0 iteraˇcním postupem
• zanebdáváme zde podmínky pro konvergenci!
• zaˇcínáme s x0 které je “dostateˇcneˇ blízko" poˇcátku • opakujeme postup dokud se xi “významneˇ mení" ˇ 1 spoˇ cteme teˇcnu v bodeˇ f (xi ) a urˇcíme pruseˇ ˚ cík s osou x 2 tento pruseˇ ˚ cík bude lepší aproximací správného ˇrešení – použijeme jako xi+1 • lze odvodit vztah xi+1 = xi −
f (xi ) f ′ (xi )
Výpoˇcet druhé odmocniny – iteraˇcní metoda 1 • y =
√
x
• ˇrešíme numerickým ˇrešením rovnice y 2 − x = 0
Newtonovou metodou
• tedy yi+1 = yi −
f (yi ) f ′ (yi )
• po dosazení dostaneme iteraˇcní vzorec yi+1 = 21 ( yx + yi )),
y0 > 0
i
• v každém kroku zhruba získáme dvojnásobek platných
míst
• metoda je nepraktická, protože obsahuje delení ˇ
• všimnete ˇ si: je to vlastneˇ babylónská metoda (ovšem
Heron samozˇrejmeˇ neznal Newtonovu metodu)
Výpoˇcet druhé odmocniny – iteraˇcní metoda 2, princip ˇ pˇri výpoˇctu, nalezneme Abychom se vyhnuli nutnosti delit pˇribližnou hodnotu √1x a pomocí algoritmu pro pˇrevrácenou √ hodnotu urˇcíme x. √ • y = x, • rˇešíme numerickým rˇešením rovnice 12 − x = 0 y Newtonovou metodou • tedy yi+1 = yi −
f (yi ) f ′ (yi )
• po dosazení dostaneme iteraˇcní qvzorec
yi+1 = 21 yi (3 − xyi2 )), 0 < y0 <
3 x
Výpoˇcet druhé odmocniny – iteraˇcní metoda 2, analýza chyby
• y =
√
x,
• pokud je yi =
√1 (1 x
+ δ) lze odvodit, že
1 √1 2 x (1
+ δ)(1 − 2δ − δ2 ) √ • necht’ ∆yi = yi − x yi+1 =
• podíly chyb yi a yi+1 lze urˇcit jako • napˇr. pro δ = 0.01 je
chybu, než yi
∆yi ∆yi+1
∆yi ∆yi+1
=
2 1 δ 3+δ
≈ 66 a yi+1 má cca 66× menší
• opet ˇ cca 2× více desetinných míst v každém kroku
Logaritmus
Proˇc logaritmus?
• pˇrevod násobení na sˇcítání log2 (A × B) = 2(log2 (A)+log2 (B)) • delení ˇ na odeˇcítání!
• vhodné pro práci s extrémneˇ velkými (x!) a malými
(Παi , |αi | < 1) cˇ ísly
• pracujeme obvykle s dvojkovým logaritmem
Obecný princip – použití tabulky pro interpolaci
0
f (xa ) = f (x0 ) +
f (x1 )−f (x0 ) x1 −x0 (xa
− x0 )
Pˇribližný výpoˇcet segmentací na úseky poˇcítáme logaritmus log2 (y) cˇ ísla y z intervalu < 0; 2N > • najdeme i takové, aby 2i ≤ y < 2i+1 (potom jisteˇ i ≤ log2 (y) < i + 1), nepotˇrebuji ani tabulku • zpˇresnení ˇ lineární aproximací: i+1
i
2 (2 ) (y − 2i ) = log2 (y) ≈ log2 (2i ) + log2 (22i+1)−log −2i −i i ˇ jen 2 – výhodné! = i − 1 + y2 , zde se delí • pˇríklad pro 8 bit. cˇ íslo
• y = 105 • 64 ≤ y < 128, i = 6 • log2 (y ) ≈ 5 + 105/64 = 6.64, pˇresná hodnota je 6.71
• pro násobení potˇrebujeme urˇcit dva logaritmy, souˇcet a
posléze odlogaritmovat výsledek (postup opaˇcný k uvedenému)
• max. chyba výsledku násobení 11.6%, logaritmus lze dále
zpˇresnit
Goniometrické funkce
Výpoˇcet hodnoty goniometrické funkce • algoritmus “s náhodným pˇrístupem" • cíl: pro zadaný úhel urˇcit hodnotu sin, cos, . . . • ruzné ˚ techniky vhodné pro HW i SW realizaci • aplikace – obvyklé výpoˇcty • generování spojitého prub ˇ ˚ ehu • cíl: generovat spojitý prub ˇ harmonické funkce, rychle, ˚ eh pˇresneˇ • algoritmus “s náhodným pˇrístupem" lze použít, ale není to cˇ asto výhodné • naopak, speciální algoritmy pro generování nejsou použitelné pro výpoˇcet “s náhodným pˇrístupem" • aplikace – generování prub ˇ u˚ (meˇ ˇ rení, buzení, digitální ˚ eh komunikace)
Tabulková metoda • v reálném cˇ ase neprovádíme
ideálneˇ žádné výpoˇcty • pˇresnost nepˇrímo úmerná ˇ velikosti tabulky
• lze interpolovat (potom cˇ asoveˇ
ˇ nároˇcnejší)
• velikost tabulky obvykle 2N
• uložíme hodnoty pro jeden
kvadrant – ROM
• triviální hardwarová
implementace
y
−,+
+,+
180 − α
α x 0
α − 180 −,−
360 − α +,−
Expanze ˇrad – Tayloruv ˚ rozvoj • vhodné spíše pro softwarovou implementaci • Taylorova ˇrada f (x) =
P∞
• Tayloruv ˚ rozvoj sinu sin(x) • cos(x) = sin(x +
π 2)
f (n) (x0 ) n n! (x − x0 ) 5 7 3 = x − x3! + x5! − x7!
n=0
+ ...
• pˇrímoˇcará implementace napˇríklad v C
• pˇresnost volím tak, že ˇrídím poˇcet sˇcítaných cˇ lenu˚
• na intervalu x ∈< −π/2; π/2 > konverguje rychleji než na
zbytku cˇ íselné osy, takže je vhodné úhel transformovat do tohoto intervalu ˇ • další techniky – aproximace Cebyševovými polynomy, Newton, Lagrange, spline interpolace, atd. – použitelné i na obecné funkce
Interpolace a aproximace – i pro jiné funkce
• Tayloruv ˚ rozvoj
CORDIC – co je to?
• COordinate Rotation Digital Computer
• technika výpoˇctu hodnot goniometrických funkcí s malými
nároky na hardware
• použitelné i pro širší škálu problému: ˚ cyklometrické funkce,
pˇrevod pravoúhlých souˇradnic na polární, hyperbolické a hyperbolometrické funkce, logaritmus, exponenciála, absolutní hodnota a fáze komplexního cˇ ísla, a další.
CORDIC – základní myšlenka • hledám x = r cos β, y = r sin β
y
• znám x1 = r cos α, y1 = r sin α • β =α+ϕ
2
• x2 = r cos β = r cos(ϕ + α) =
ϕ 1
r cos ϕ cos α − r sin ϕ sin α, y2 = r sin β = r sin(ϕ + α) = r sin ϕ cos α + r cos ϕ sin α
β α
0
x 1
• x2 = x1 cos αi+1 − y1 sin ϕ =
cos ϕ(x1 − y1 tan ϕ)
• y2 = x1 cos αi+1 + y1 sin ϕ =
cos ϕ(x1 + y1 tan ϕ)
CORDIC – základní myšlenka
y
• chci zjistit cos α a sin α,
α ∈< 0; 90 >
i+1 αi+1
• intuitivneˇ – rotace vektoru po
krocích až do dosažení cílové pozice
i
• xi+1 = cos αi+1 (xi − yi tan αi+1 ),
yi+1 = cos αi+1 (xi + yi tan αi+1 )
α1 + α2 + . . . + αi 0
x r
• použijeme tabelované hodnoty
αi+1 spolu s jejich tangenty
CORDIC – tabulka hodnot xi+1 = cos αi+1 (xi − yi tan αi+1 ), yi+1 = cos αi+1 (xi + yi tan αi+1 ) 1 – redukce množství násobení Volím αi aby tan αi = 2i−1 Poˇcet hodnot v tabulce daný žádanou pˇresností vyjádˇrení α α ≈ α′ = α0 ± α1 ± α2 . . . ± αN−1 i 0 1 2 3 4 5 6 .. .
tan αi 1.0000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156 .. .
αi 45.00 26.57 14.04 07.13 03.58 01.79 00.90 .. .
cos αi 0.7071 0.8944 0.9701 0.9923 0.9981 0.9995 0.9999 .. .
N −1
2−(N−1)
...
...
CORDIC – tabulka hodnot, pˇríklad α ≈ α′ = α0 ± α1 ± α2 . . . ± αN−1 α = 38′ i i 0 1 2 3 4 5 6
tan αi
αi
cos αi
±
1.0000 0.5000 0.2500 0.1250 0.0625 0.0313 0.0156
45.00 26.57 14.04 07.13 03.58 01.79 00.90
0.7071 0.8944 0.9701 0.9923 0.9981 0.9995 0.9999
00.00<38,+ 45.00>38,18.43<38,+ 32.47<38,+ 39.60>38,36.02<38,+ 37.81>38,+
α′ 0 45.00 18.43 32.47 39.60 36.02 37.81 38.71
CORDIC – iteraˇcní výpoˇcet • α ≈ α′ = α0 ± α1 ± α2 . . . ± αN−1
• rekurentní vztah xi+1 = cos αi+1 (xi ∓ yi tan αi+1 ), x0 = 1,
y0 = 0 yi+1 = cos αi+1 (xi ± yi tan αi+1 )
• násobení tan αi+1 je nyní jen bitový posuv
• násobení cos αi+1 budeme behem ˇ iterace ignorovat,
použijeme rekurentní vztah xi+1 = (xi ∓ yi tan αi+1 ), yi+1 = (xi ± yi tan αi+1 )
• výsledek pak bude tˇreba na konci vynásobit faktorem
K = cos α0 × cos α1 × . . . cos αN−1 , ten je ale konstantní dokud je konstantní poˇcet kroku˚ – lze ho pˇredpoˇcítat dopˇredu
• funguje pro α ∈< 0; 90 >
CORDIC – rozšíˇrení na celou kružnici y
−,+
180 − α
α
+,+ x
0
α − 180 −,−
360 − α +,−
Goniometrické funkce – generování
Použití IIR filtru
• rezonátor s póly na jednotkové kružnici
• y[n] = −a1 y[k − 1] − a2 y[k − 2] + b0 x[k] • b0 = A sin(ω), a1 = −2 cos(ω), a2 = 1 • ω=
2πf fs
• výhoda – extrémneˇ jednoduché
• nevýhoda – postupná akumulace numerických chyb
(kvantovací efekty) “niˇcí" kvalitu generované sinusovky
DDS – vlastnosti Direct Digital Synthesis
Výhody • vysoké frekvenˇcní rozlišení (až µHz) • malé fázové zkreslení
• jednoduchá hardwarová implementace
• žádné ustalování systému – rychlé zmeny ˇ frekvence • v tabulce lze mít libovolnou periodickou funkci
Nevýhody • pro velkou šíˇrku pásma potˇrebuji rychlý hardware
• pro velkou šíˇrku pásma bude vysoká spotˇreba energie
• muže ˚ být problém s odezvou (rychlostí) DA pˇrevodníku
Pˇrímá cˇ íslicová syntéza – základní princip
Figure copied from http://www.cs.washington.edu/homes/diorio/Talks/InvitedTalks/MTT97/index.html
Faze harmonicke
Fazovy prirustek
DDS – akumulátor
• generuje fázovou informaci; frekvence pˇretékání urˇcuje fout • generovaná frekvence je fout =
FP×fclk 2N
• poˇcet bitu˚ akumulátoru dává rozlišení ve frekvenci • frekvenˇcní krok ∆f =
fclk 2N
• pro fclk = 100MHz a N = 32 pˇr. ∆f = 23mHz
DDS – akumulátor
• zmenou ˇ ˇrídím zmenu ˇ frekvence ˚ fázového pˇrírustku • okamžite, ˇ bez pˇrechodového deje ˇ • pˇri zmen ˇ eˇ frekvence spojitá funkce (fáze se nemení) ˇ • jednoduché generování napˇríklad FSK signálu (alternace ˚ u) ˚ dvou fázových pˇrírustk
DDS – konvertor fáze–amplituda
• pˇresnost konverze jen taková, aby zkreslení výstupu bylo
urˇceno hlavneˇ DA pˇrevodníkem a následným zpracováním ˇ ˇ než umožnuje analogový (nemá smysl být pˇresnejší, výstupní díl!) • mužeme ˚ použít libovolný algoritmus s “náhodným pˇrístupem" • mimoˇrádneˇ vhodná je tabulková implementace (LUT, Look Up Table)
DDS – LUT konvertor fáze–amplituda • výstup akumulátoru – fázová informace má N bitu˚
• teoreticky potˇrebuji 2N položek v pameti, ˇ budu indexovat
fází a vyˇctu hodnotu sin
• staˇcí mi jen
2N 4
– symetrie harmonické po kvadrantech, mírná komplikace s obracením znamének (ale žádný problém)
• nekdy ˇ muže ˚ mít i méneˇ položek (velké rozlišení cˇ ítaˇce
kvuli ˚ malému frekvenˇcnímu kroku, malé rozlišení generované sinusovky kvuli ˚ horšímu analogovému výstupnímu dílu a DA pˇrevodníku)
• zmenšení poˇctu položek lze dosáhnout napˇríklad
interpolací
DDS – LUT konvertor fáze–amplituda • další možnost redukce velikosti tabulky goniometrickými
identitami
• sin(ϕ) = sin(α ∗ 2L + β) =
sin(α ∗ 2L ) cos(β) + cos(α ∗ 2L ) sin(β)
• rozdelím ˇ fázi na dveˇ slova, F = FMSB ∗ 2L + FLSB • vytvoˇrím dveˇ tabulky hodnot sin o
2N−L 4
a
2L 4
prvcích
• podle identity složím dohromady výslednou hodnotu
• menší tabulky (pro N = 16 místo 16384 položek v jedné
tabulce pro L = 8 staˇcí 64 a 64 položek!), snadný výpoˇcet ˇ bez podmínených operací
• dílˇcí tabulky zˇrejmeˇ musí být s vetší ˇ šíˇrkou slova
