用 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 數列, 實在是有趣。

--