定點數表示法

筆者今天要來介紹定點數表示法, 也就是 fixed-point representation. 讀者可能會很納悶, 現在有浮點數的機器非常普遍, 而且速度都很快, 比如說 Intel Core i7 處理器,開啟 AVX SIMD 的情況下, 可以一次處理四個 64 位元浮點數的乘法。但是在很多地方, 比如說訊號處理、控制系統、或者很多手機內部的處理器, 都欠缺浮點運算能力,這時候就必須使用定點數。...全文 PDF

用 Linear Filter(IIR) 來實作 Fibonacci 序列

PyCon.TW 2012 在日前圓滿落幕, 會中的 keynote speaker Travis Oliphant 講
了一個使用 linear filter 實作 Fibonacci 的例子, 在演講中只有一張投影片,
在這裡我嘗試做分解動作解釋一下。投影片內容在 http://www.slideshare.net/pycontw/largescale-arrayoriented-computing-with-python

我們知道 Fibonacci 序列是 x(n) = x(n-1) + x(n-2), 而 x(0) = 0, x(1)
= 1. Fibonacci 序列的產生也常常是教科書介紹 recursive 的好範例, 不過就
如 Travis Oliphant 所示, 用 recursive 實作的複雜度是 exponential 成長的,
而我們也可以用 SciPy 裡面的 linear filter 來實作這件事, 方法如下：

from scipy.signal import lfilter
from numpy import zeros
b = array([1.0])
a = array([1., -1, -1])
zi = array([0, 1.0])
def fib3a(N):
   y, zf = lfilter(b, a, zeros(N, dtype=float), zi=zi)
   return y

我來慢動作分解一下上面的實作方法。首先 lfilter 在 b = array([1.0]) 以及
a=array([1.,-1,-1]) 的情況下會產生b(0)x(n) = a(0)y(n) + a(1)y(n-1) +
a(2)y(n-2) 的式子, 代入 a,b, 我們可以得到 x(n) = y(n) - y(n-1) -
y(n-2),而 lfilter 的第三個參數就是 x 的序列, 是 zeros(N, dtype=float),
也就是說輸入全為零,上述的式子就會變成 0 = y(n) - y(n-1) - y(n-2), 把 y(n) 移到等式的左方, 就可以得到 y(n) = y(n-1) + y(n-2), 這就是 Fibonacci 數列的表示法啦! 最後的 y(n) 就是 Fibonacci 的第 n 個數列。另外,
zi=array([0,1.0]) 就是當 lfilter 開始執行的時候, delay element 裡面的元
素, 也就是 y(0) 與 y(1) 的值。

以濾波器的角度而言, 這個濾波器在輸入全為零的狀況下, 自己產生
0,1,1,2,3,5,… 這種發散的數列, 等於自己在震盪, 是我們在實際上不會去使
用的濾波器, 拿來實作 Fibonacci 數列, 實在是有趣。

--

LaTeX PSTricks 訊號處理範例: z-Transform

上一篇的 moving average filter y(n) = (15/16)*y(n-1) + (1/16)*x(n) 的
z-transform 有一個 pole 在 z=(15/16, 0) 的點上, 要畫這個 z-Transform 在
z-plane 上的表示圖, 可以用下列的 PSTricks 來表示。

\usepackage{pstricks}
\usepackage{pst-sigsys}

圖形如下表示:

\begin{center}
  \begin{pspicture}[showgrid](-2,-2)(2,2)
    \pscircle[linecolor=gray](0,0){1} % unit circle
    \pspole(0.9375,0){z1}
    \nput{0}{z1}{$(\frac{15}{16},0)$}
  \end{pspicture}
\end{center}

--

LaTeX PSTricks 訊號處理範例: Moving Average Filter

今天要介紹的是如何用 PSTricks 來繪製訊號處理的 functional diagram, 假設
我們要描述的訊號處理方程式是: y(n) = (15/16)*y(n-1) + (1/16)*x(n). 這是
一個動平均濾波器, y(n) 是目前的輸出值, x(n) 是目前的輸入值, y(n-1) 則是上一次的輸出值,
那他的 functional diagram 會是如下圖所示:

上述圖形使用了 PSTricks 以及 PSTricks 的 pst-sigsys 套件，所以要在 preamble 的部份如下宣告:

\usepackage{pstricks}
\usepackage{pst-sigsys}

描述這個濾波器的 LaTeX source code 如下:

\begin{center}
  % To describe y(n) = (15/16)*y(n-1) + (1/16)*x(n)
  \begin{pspicture}(-2,-1)(10,2)
    \pssignal(0,1){x}{$x(n)$}
    \pscircleop[oplength=0.25,operation=times](2,1){op1}
    \pssignal(2,0){coefx}{$\frac{1}{16}$}
    \pscircleop[oplength=0.25](4,1){op2}
    \pssignal(8,1){y}{$y(n)$}
    \dotnode(6,1){ydot}
    \psblock(6,0){delay}{$z^{-1}$}
    \pscircleop[oplength=0.25,operation=times](4,0){op3}
    \pssignal(3,0){coefy}{$\frac{15}{16}$}
    \nclist{->}{ncline}{x,op1,op2,y}
    \nclist{->}{ncline}{coefx,op1}
    \nclist{->}{ncline}{ydot,delay,op3,op2}
    \nclist{->}{ncline}{coefy,op3}
  \end{pspicture}
\end{center}

--

gcc for Andes compile 出來的 assembly

為了要做 24-bit signed integer to 32-bit signed integer 的 signed extension 測試，我寫了以下的一小段 code 來測試, 在這裡順便稍微提一下 optimization, 我們先看 c code:

void testSign()
{
        unsigned int test=0x00800000;
        int test_s;
        test_s = ((int)test << 8 ) >> 8;
}

除去 prologue 或 epilogue 的部份，我們來看主體:

00000614 <testSign>:
......
 61c:   46 00 08 00     sethi $r0,#2048
 620:   14 0e 7f fd     swi $r0,[$fp+#-12]
 624:   04 0e 7f fd     lwi $r0,[$fp+#-12]
 628:   40 00 20 08     slli $r0,$r0,#0x8
 62c:   90 08           srai45 $r0,#0x8
 62e:   14 0e 7f fe     swi $r0,[$fp+#-8]

解釋如下：

sethi $r0,#2048

$r0 = 0x00800000

swi $r0,[$fp+#-12]

把 $r0 存回 test 的儲存空間

lwi $r0,[$fp+#-12]

把變數 test 值抓到 $r0

slli $r0,$r0,#0x8

把 $r0 做 logical left shift by 8 bits, $r0 = 0x80000000

srai45 $r0,#0x8

把 $r0 做 arithmetic right shift by 8 bits, $r0 = 0xff800000

swi $r0,[$fp+#-8]

存回 test_s 的儲存空間。

而如果我們把這個 c code 改成如下：

int testSign(unsigned int value)
{
        return ((int)test << 8) >> 8;
}

產生出來的 assembly 就會變成:

000005f8 <testSign>:
 5f8:   40 00 20 08     slli $r0,$r0,#0x8
 5fc:   90 08           srai45 $r0,#0x8
 5fe:   dd 9e           ret5 $lp

這個版本完全沒有任何 local variables, 傳入的參數也只有一個, 當傳入參數
只有一個的時候, 會使用 $r0 這個暫存器傳入, 回傳值也是在 $r0 回傳,
所以這樣子最省空間了,也是我們 assembly programmer 會寫出來的樣子。注意
如果沒有做 (int) 的 explicit cast, 編譯器並不會使用 srai45 這個
arithmetic right shift 指令製造 signed extension, 即便是最後的回傳值是
有號數。

--
My Emacs Files At GitHub