Polynomials and FFT

因爲計算機只能表示離散的信號，所以不少連續的函數與信號在計算機中不得不採用離散信號來近似。由Tylor公式等能夠知道通常平滑連續信號均可以表示爲一個N項的Tylor多項式，該多項式能夠在任意精度上擬合平滑信號。所以Polynomials在計算機計算中顯的由爲重要。

1.Polynomials

1.1 Polynomials的通用形式：

\(A(x) = \sum_{j=0}^{n-1}a_jx^j\)

對於等號右邊具備n項的Polynomial，能夠用一個係數向量來簡單表示，即

\((a_0,a_1,a_2,...,a_{n-1})\)

能夠用此向量來計算關於多項式的各類操做，如加減乘除等。假設以上多項式的最高項爲\(x^{n-1}\),且其係數不爲0，能夠說這個多項式的度(degree)爲n-1,degree-bound是另一個數，設爲N，是多項式

degree的一個上界，即

\(Degree(n) is always <= c * Degree-bound(n)\).

1.2 係數向量上的操做

若是求上面這樣一個多項式在某個點處的值，能夠在線性時間完成：

\(y_0 = a_0+x_0*(a_1+x_0*(a_2+x_0*(a_3+(...))));\)

若是有兩個多項式係數向量:

\(A = (a_0,a_1,a_2,...,a_{n-1});\)

\(B = (b_0,b_1,b_2,...,b_{n-1});\)

在此基礎上求多項式,\(C=(c_0,c_1,c_2,...,c_{n-1});\)，C若是爲上面兩個多項式的和，那麼計算是方便的，僅需計算\(c_j = a_j+b_j\),僅需線性的時間，即O(n)，且

\(degree(c) = max(degree(a),degree(b))\);

若是c是這兩個多項式的積，計算就沒這麼簡單了，兩個多項式相乘，須要O(n^2)的時間複雜度。用公式表示這樣的過程：\(C=A \otimes B\)

平方級的時間複雜度通常是咱們但願優化的，由於它能處理數據的上限僅爲百萬的數據量。

1.3 Polynomials的另外一種表示:Point-value

若是對於通用形式給出的多項式，咱們有其在n個點處的值，即:

\(\{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\};\)

是否是也能夠惟一的表示一個Polynomial;答案是能夠的，由於有以上n對點值，咱們能夠寫出一個線性方程組，並轉換爲矩陣形式：

左邊是一個範的蒙矩陣，這個線性方程組若是有惟一解，即\((a_0,a_1,...,a_{n-1});\)列向量，就獲得了通用形式的係數向量；事實上這個假設是明顯成立的，由於範德蒙矩陣

的行列式不爲0，方程組有惟一的解：\(a = V(x_0,x_1,...,x_{n-1})^{-1}y\);

也就是說，咱們能夠用兩種方式來惟一的表示一個通用形式的多項式，對於第二種Point-value的表示方法，考察其上的多項式操做：

仍然以A,B,C多項式爲例，這裏的A,B是基於Point-value表示的：

\(A = \{(x_0,y_0),(x_1,y_1),...,(x_{n-1},y_{n-1})\};\)

\(B = \{(x_0,z_0),(x_1,z_1),...,(x_{n-1},z_{n-1})\};\)

\(C = \{(x_0,t_0),(x_1,t_1),...,(x_{n-1},t_{n-1})\};\)

若是\(C=A+B\)，僅需計算每項的value值，即\(t_j=y_j+z_j\),時間複雜度O(n),而若是\(C=A \otimes B\) ,仍然僅需計算value值，此處爲\(t_j = y_j * z_j\),在O(n)的時間內完成。

看起來貌似使多項式乘法的複雜度從平方級降到了線性，但其實遠非沒有那麼簡單，首先，有上述計算\(C=A \otimes B\)的方式獲得的C是不可用的，由於按此方式計算出的C的Point-value表示

僅有n對點值，與A,B是相同的，但僅僅計算\((x+1)*(x+1) = x^2+2x+1\)就能夠看出問題，左邊兩個乘數多項式的degree均爲1，所以能夠用有一個長度爲2的係數向量\(a_0,a_1\)表示，

也就能夠用2個點值對錶示。

等號右邊卷積結果爲2，係數向量長度爲3,\((a_0,a_1,a_2)\),所以必須至少用3個點值對。

解決這個問題的方法很簡單，由於Degree_bound(n)*Degree_bound(n) = Degree_bound(2n),也就是說，兩個degree_bound爲n的多項式的卷積結果的degree不會大於2n，係數向量的長

度也不會超過2n，所以徹底可使用2n個點值對，乘數也用2n個點值對，來表示長度爲2n的係數向量，原係數向量長度僅爲n，超出的部分使用0填充。

2.加速多項式乘法：

如上述所述，若是以point-value表示一個多項式，是能夠線性的時間內計算兩個point-value表示的多項式的卷積結果(仍然以point-value)表示，但實際的狀況並不是如此美好，通常狀況下不會

給定一個Polynomial的point-value表示，一般能夠獲得的是其係數向量的表示形式，並且卷積結果通常也要求以係數向量表示。畢竟係數向量是惟一的，且是一個Polynomial的直接體現。

但在係數向量表示下，是沒法加速卷積的。

能夠採用這樣的策略，獲得給的的係數向量後，轉換爲相應的point-value表示，在point-value下計算卷積結果的point-value,而後將結果再轉換爲係數向量；

這個策略是可行的，係數向量到point-value的轉換成爲evaluate,由point-value到係數向量的轉換稱爲interpolation，即插值。前者爲連續到離散的過程，後者爲離散到連續的過程。

注：插值在計算機領域的應用是很是普遍的，例如在計算機圖形學三維重建中，有離散點集重建連續平滑的三維表面，就須要大量的插值。經常使用的插值算法有拉格朗日插值算法等;

可行方案使人鼓舞的同時仍然有一個重要的問題亟待思考：evaluate和interpolation過程的時間複雜度? 若是evaluate或interpolation的複雜度高於O(n^2),那麼整個方案將是徒勞的;

2.1 evaluate and interpolation

若是給定n的係數向量，要獲得n的point-value,只須要選擇n個離散點point，計算Polynomial在該點處的值，計算每一個點時間複雜度爲O(n),整個evaluate就是O(n^2);而interpolation是evaluate

的反過程，時間複雜度也可爲O(n^2);

如此，轉換瓶頸限制了整個加速策略的時間複雜度，所以目標就是但願可以加快轉換速度；

2.2 選擇特殊的點集

可以加速evaluate和interpolation的辦法是選擇特殊的點，獲得特殊的Point-value表示，在這些點上的轉換可以在時間複雜度O(nlgn)時間內完成;整個過程的大體步驟以下:

圖2 加速乘法過程

 

如上圖所示，選擇在特殊點\((w_{2n}^0,w_{2n}^1,...,w_{2n}^{2n-1})\)進行evaluation和interpolation,\(w_n^k\)是這樣的點：

\( w_n^k = e^{2{\pi}ki/n} = cos(2{\pi}k/n)+i*sin(2{\pi}k/n)\)

\((w_n^k)^n = 1\)

w爲1的n次復根(complex nth root of unity),具備週期性和一些對稱的性質,如圖:

圖3：8個1的8次復根

列出其一些特有的性質，能夠自行證實：

\(w_{d*n}^{d*k} = w_n^k\)

\((w_n^k)^2 = w_n^{2k} = w_{n/2}^k\)

\(w_n^{k+n/2} = -w_n^k\)

其餘一些性質也能夠從圖中看出;

因而，在evalutaion過程當中,只需計算\(y_k = A(w_n^k)\),如過將\( w_n^k = e^{2{\pi}ki/n}\)代入,將發現這就是DFT離散傅立葉變換公式，至於二者之間的淵源，能夠追本溯源的考察一番，此處不作贅述；

表面上看evaluation仍須要n次計算point，每次O(n),但能夠改進計算的過程，利用特殊點對稱週期等的特性，將O(n^2)降到O(nlgn);

從新明確一下目標：在O(nlgn)時間內完成計算:

\(y_k = A(w_n^K) = \sum_{j=0}^{n-1} a_j*w_n^{kj},  k = 0,1,2,...,n-1\);

若是用分治的方法計算結果，能夠把計算過程作這樣的安排：

\(y_k = A(x) = A^{[0]}(x^2)+x*A^{[1]}(x^2);\)

\(A^{[0]}(x) = a_0+a_2x+a_4x^2+...+a_{n-2}x^{n/2-1};\)

\(A^{[1]}(x) = a_1+a_3x+a_5x^2+...+a_{n-1}x^{n/2-1};\)

\(y_k = A(w_n^k) = A^{[0]}(w_n^{2k})+w_n^k*A^{[1]}(w_n^{2k});\)

 以n = 8爲例:

\(y = A(x) = a_0+a_1x+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6+a_7x^7=(a_0+a_2x^2+a_4x^4+a_6x^6)+x(a_1+a_3x^2+a_5x^4+a_7x^6)\)

              \( = (((a_0+a_4x^4)+x^2(a_2+a_6x^4))+x((a_1+a_5x^4)+x^2(a_3+a_7x^4))) ;\)

即將全多項式求值過程拆成全部偶數項和奇數項的和,計算偶數項或奇數項時能夠按此方法繼續拆，知道多項式只有一項，設計算整個多項式值的時間複雜度爲\(T(n)\),則有:

\(T(n) = T(n/2) + n/2;\)

\(T(1) = 1;\)

所以，\(T(n) = O(nlgn); \),若是在僅計算\(x_0\)處多項式的值就需O(nlgn),那麼計算完n個採樣點，須要O(n^2lgn),比不採用分治花費更多的時間，彷佛並非咱們所須要的，確實在以前求多項式在某一點值也並無採用此法，但若是把某一點變爲特殊的點，如\(w_n^k\),狀況就大不相同了,在此情形下，不用逐個計算\(y_k\),採用蝶形算法，能夠在一個迭代內計算全部\(y_k\),時間複雜度O(nlgn);

選擇特殊點後:

\(y_k = a_0 + a_1w_8^k + a_2w_8^{2k} + a_3w_8^{3k} + a_4w_8^{4k} + a_5w_8^{5k} + a_6w_8^{6k} + a_7w_8^{7k} \)

\(      = (((a_0 + a_4w_8^{4k}) + w_8^{2k}(a_2 + a_6w_8^{4k})) + w_8^k(( a_1 + a_5w_8^{4k}) + w_8^{2k}(a_3 + a_5w_8^{4k}))) \)

\(      = (((a_0 + a_4w_2^k) + w_4^k(a_2 + a_6w_2^k)) + w_8^k((a_1 + a_5w_2^k) + w_4^k(a_3 + a_5w_2^k))) \)

上式給出了\(y_k,k=0,1,2...7\)的通項形式，不一樣的\(y_k\)均可由上式獲得，由於\(w_n^k\)是以n爲週期的,而且半週期正負易號，因此上述公式中不一樣的\(y_k\)計算存在不少相同的成分,這些成分是不須要重複計算的。

繼續回到分治上來，上述公式是由分治策略推導而來的,因此\(y_k\)應該看做是8個長度爲1的項的歸併，公式中也是這樣體現的:首先\((a_0,a_4),(a_2,a_6),(a_1,a_5),(a_3,a_7)\)一一歸併，歸併係數\(w_2^k\),

接着\(((a_0,a_4),(a_2,a_6)),((a_1,a_5),(a_3,a_7))\)兩兩歸併，歸併係數\(w_4^k\),最後四四歸併,歸併係數\(w_8^k\),對於不一樣的k,\(y_k\)的不一樣只取決以不一樣

的歸併係數,因爲歸併係數具備週期性,半週期易號的特性,因此不一樣的k,歸併過程有有不少相同的歸併步驟，因此能夠無需逐個計算\(y_k\);

圖4  歸併係數半週期易號

上圖是歸併的圖形表示,其中用到的歸併係數半週期易號的特性，例如對於全部偶數k,歸併有\(a_0 + a_4w_2^0 = a_0 + a_4\),對於奇數,則有\(a_0+a_4w_2^1=a_0 - a_4\);

使用自底向上的歸併,即略去切分的過程,假設已經知道歸併的方向,那麼向上的歸併過程能夠用下圖描述:

圖5 分治過程

 

圖6 自底向上的蝶形歸併過程

 圖6展現瞭如何由自底向上的歸併在\(O(nlgn)\)時間內完成evaluation;

對於interpolation過程，其實和evaluation是相似的,

\(a = V(x_0,x_1,...,x_{n-1})^{-1} * y;\) interpolation

\(y = V(x_0,x_1,...,x_{n-1}) * a;\) evaluation

 對於特殊點\(w_n^k = e^{2{\pi}ki/n}\),其逆矩陣只需將\(w_n^k\)改成\(w_n^k = e^{-2{\pi}ki/n}\)並數乘\(1/n\),計算過程是徹底相同的;

其實上述evaluation就是DFT離散傅立葉轉換過程,使用蝶形歸併計算就是著名的FFT快速傅立葉算法,interpolation即逆DFT;係數向量a能夠看做時域(或空域)信號序列,y看做頻域信號序列;

注:上述算法過程假設了序列長度是的2的指數,對於長度非2的指數的序列,可使用零填充成2的指數,就整理本文時,還未對非填充序列的FFT進行了解;

注:其C++實如今另外一篇博文二維離散傅立葉變換圖像頻域處理有完整描述;關於其應用FFT的二維擴展即並行改進將在下一篇文章中詳述;