十五、多項(xiàng)式乘法與快速傅里葉變換
經(jīng)典算法研究系列,已經(jīng)寫到第十五章了,本章,咱們來介紹多項(xiàng)式的乘法以及快速傅里葉變換算法。本博客之前也已詳細(xì)介紹過離散傅里葉變換(請(qǐng)參考:十、從頭到尾徹底理解傅里葉變換算法、上,及十、從頭到尾徹底理解傅里葉變換算法、下),這次咱們從多項(xiàng)式乘法開始,然后介紹FFT算法的原理與實(shí)現(xiàn)。同時(shí),本文雖涉及到不少數(shù)學(xué)公式和定理(當(dāng)然,我會(huì)盡量舍去一些與本文咱們要介紹的中心內(nèi)容無關(guān)的定理或證明,只為保證能讓讀者易于接受或理解),但盡量保證通俗易懂,以讓讀者能看個(gè)明白。
有朋友建議,算法專一種,就ok,沒必要各個(gè)都學(xué)習(xí)。但個(gè)人實(shí)在抑制不住自己的興趣,就是想寫,當(dāng)沒法做到強(qiáng)迫自己不寫時(shí),這個(gè)經(jīng)典算法研究系列便一直這么寫下來了。
我們知道,有兩種表示多項(xiàng)式的方法,即系數(shù)表示法和點(diǎn)值表示法。什么是系數(shù)表示法?所謂的系數(shù)表示法,舉個(gè)例子如下圖所示,A(x)=6x^3 + 7x^2 - 10x + 9,B(x)=-2x^3+4x-5,則C(x)=A(x)*B(x)就是普通的多項(xiàng)式相乘的算法,系數(shù)與系數(shù)相乘,這就是所謂的系數(shù)表示法。
那么,何謂點(diǎn)值表示法?。一個(gè)次數(shù)為n次的多項(xiàng)式A(x)的點(diǎn)值表示就是n個(gè)點(diǎn)值所形成的集合:{(x0,y0), (x1,y1),..., (xn-1,yn-1)}。其中所有xk各不相同,且當(dāng)k=0,1,……,n-1時(shí),有y(k)=A(xk)。
一個(gè)多項(xiàng)式可以由很多不同的點(diǎn)值表示,這是由于任意n個(gè)相異點(diǎn)x0,x1,...,xn-1組成的集合,都可以看做是這種點(diǎn)值法的表示方法。對(duì)于一個(gè)用系數(shù)形式表示的多項(xiàng)式來說,在原則上計(jì)算其點(diǎn)值表示是簡單易行的,我們只需要選取n個(gè)相異點(diǎn)x0,x1,...,xn-1,然后對(duì)k=0,1,...,n-1,求出A(xk),然后根據(jù)霍納法則,求出這n個(gè)點(diǎn)的值所需要的時(shí)間為O(n^2)。
然,稍后,你將看到,如果巧妙的選取xk的話,適當(dāng)?shù)睦命c(diǎn)值表示可以使多項(xiàng)式的乘法可以在線性時(shí)間內(nèi)完成。所以,如果我們能恰到好處的利用系數(shù)表示法與點(diǎn)值表示法的相互轉(zhuǎn)化,那么我們可以加速多項(xiàng)式乘法(下面,將詳細(xì)闡述這個(gè)過程),達(dá)到O(n*logn)的最佳時(shí)間復(fù)雜度。
前面說了,當(dāng)用系數(shù)表示法,即用一般的算法表示多項(xiàng)式時(shí),兩個(gè) n次多項(xiàng)式的乘法需要用 O(n^2)的時(shí)間才能完成。但采用點(diǎn)值表示法時(shí),多項(xiàng)式乘法的運(yùn)行時(shí)間僅為O(n)。接下來,咱們要做的是,如何利用系數(shù)表示法與點(diǎn)值表示法的相互轉(zhuǎn)化來實(shí)現(xiàn)多項(xiàng)式系數(shù)表示法的O(n*logn)的快速乘法。
在此之前,我們得明白兩個(gè)概念,求值與插值。通俗的講,所謂求值(系數(shù)->點(diǎn)值),是指系數(shù)形式的多項(xiàng)式轉(zhuǎn)換為點(diǎn)值形式的表示方式。而插值(點(diǎn)值->系數(shù))正好是求值的逆過程,即反過來,它是已知點(diǎn)值表示法,而要你轉(zhuǎn)換其為多項(xiàng)式的系數(shù)表示法(用n個(gè)點(diǎn)值對(duì)也可以唯一確定一個(gè)不超過n-1次多項(xiàng)式,這個(gè)過程稱之為插值)。而這個(gè)系數(shù)表示法與點(diǎn)值表示法的相互轉(zhuǎn)化,即無論是從系數(shù)表示法轉(zhuǎn)化到點(diǎn)值表示法所謂的求值,還是從點(diǎn)值表示法轉(zhuǎn)化到系數(shù)表示法所謂的插值,求值和插值兩種過程的時(shí)間復(fù)雜度都是O(n*logn)。
前面,我們已經(jīng)說了,在多項(xiàng)式乘法中,如果用系數(shù)表示法表示多項(xiàng)式,那么多項(xiàng)式乘法的復(fù)雜度將為O(n^2),而用點(diǎn)值表示法表示的多項(xiàng)式,那么多項(xiàng)式的乘法的復(fù)雜度將是線性復(fù)雜度為O(n),即: 適當(dāng)?shù)睦命c(diǎn)值表示可以使多項(xiàng)式的乘法可以在線性時(shí)間內(nèi)完成。此時(shí)讀者可以發(fā)揮你的想象,假設(shè)我們以下面這樣一種過程來計(jì)算多項(xiàng)式的乘法(不過在此之前,你得先把兩個(gè)要相乘的多項(xiàng)式A(x)和B(x)擴(kuò)充為次數(shù)或者說系數(shù)為2n次的多項(xiàng)式),我們將會(huì)得到我們想要的結(jié)果:
那么,綜上上述三項(xiàng)操作,系數(shù)表示法表示的多項(xiàng)式乘法總的時(shí)間復(fù)雜度已被我們降到了O(n*logn+n+n*logn)=O(N*logN)。
如下圖所示,從左至右,看過去,這個(gè)過程即是完成多項(xiàng)式乘法的計(jì)算過程。不過,完成這個(gè)過程有兩種方法,一種就是前面第一節(jié)中所說的普通方法,即直接對(duì)系數(shù)表示法表示的多項(xiàng)式進(jìn)行乘法運(yùn)算,復(fù)雜度為O(n^2),它體現(xiàn)在下圖中得Ordinary multiplication過程。還有一種就是本節(jié)上文處所述的三個(gè)步驟:1、將多項(xiàng)式由系數(shù)表示法轉(zhuǎn)化為點(diǎn)值表示法(點(diǎn)值過程);2、 利用點(diǎn)值表示法完成多項(xiàng)式乘法;3、將點(diǎn)值表示法再轉(zhuǎn)化為系數(shù)表示法(插值過程)。三個(gè)步驟下來,總的時(shí)間復(fù)雜度為O(N*logN)。它體現(xiàn)在下圖中的Evaluation + Pointwise multiplication + Interpolation 三個(gè)合過程。
由上圖中,你也已經(jīng)看到了。我們巧妙的利用了關(guān)于點(diǎn)值形式的多項(xiàng)式的線性時(shí)間乘法算法,來加快了系數(shù)形式表示的多項(xiàng)式乘法的運(yùn)算速度。在此過程中,我們的工作最關(guān)鍵的就在于以O(shè)(n*logn)的時(shí)間快速把一個(gè)多項(xiàng)式從系數(shù)形式轉(zhuǎn)換為點(diǎn)值形式(求值,我們即稱之為FFT),然后再以O(shè)(n*logn)的時(shí)間從點(diǎn)值形式轉(zhuǎn)換為系數(shù)形式(插值,我們即稱之為FFT逆)。最終達(dá)到了我們以最終的系數(shù)形式表示的多項(xiàng)式的快速乘法為O(N*logN)的時(shí)間復(fù)雜度。好不令人心生快意。
對(duì)上圖,有一點(diǎn)必須說明,項(xiàng)w2n是2n次單位復(fù)根。且不容忽視的是,在上述兩種表示法即系數(shù)表示法和點(diǎn)值表示法相互轉(zhuǎn)換的過程中, 都使用了單位復(fù)根,才得以在O(n*logn)的時(shí)間內(nèi)完成求值與插值運(yùn)算(至于何謂單位復(fù)根,請(qǐng)參考相關(guān)資料。因?yàn)闉榱吮WC本文的通俗易懂性,無意引出一大堆公式或定理)。
注:本節(jié)有相當(dāng)部分文字來自算法導(dǎo)論中文版第二版以及維基百科。
在具體介紹FFT之前,我們得熟悉知道折半定理是怎么一回事兒:如果n>0且n為偶數(shù),那么,n個(gè)n次單位復(fù)根的平方等于n/2個(gè)n/2個(gè)單位復(fù)根。我們已經(jīng)知道,通過使用一種稱為快速傅里葉變換(FFT)的方法,可以在O(n*logn)的時(shí)間內(nèi)計(jì)算出DFTn(a),而若采用直接的方法復(fù)雜度為O(n^2)。FFT就是利用了單位復(fù)根的特殊性質(zhì)。
你將看到,F(xiàn)FT方法運(yùn)用的策略為分治策略,所謂分治,即分而治之。兩個(gè)多項(xiàng)式A(x)與B(x)相乘的過程中,F(xiàn)FT用A(x)中偶數(shù)下標(biāo)的系數(shù)與奇數(shù)下標(biāo)的系數(shù),分別定義了兩個(gè)新的次數(shù)界為n/2的多項(xiàng)式A[0](x)和A[1](x):
A[0](x) =a0 +a2x +a4x2 +··· +an-2xn/2-1,
A[1](x) =a1 +a3x +a5x2 +··· +an-1xn/2-1.
注意,A[0]包含了A中所有偶數(shù)下標(biāo)的系數(shù)(下標(biāo)的相應(yīng)二進(jìn)制數(shù)的最后一位為0),而A[1]包含A中所有奇數(shù)下標(biāo)的系數(shù)(下標(biāo)的相應(yīng)二進(jìn)制數(shù)的最后一位為1)。有下式:
這樣,求A(x)在
1)求次數(shù)界為n/2的多項(xiàng)式A[0]與A[1]在點(diǎn)
2)將上述結(jié)果進(jìn)行組合。
下面,我們用N次單位根WN來表示
WN的性質(zhì):
為了簡單起見,我們下面設(shè)待變換序列長度n = 2r。 根據(jù)上面單位根的對(duì)稱性,求級(jí)數(shù)
Fodd(k) 和 Feven(k)是兩個(gè)分別關(guān)于序列
這樣,一個(gè)N點(diǎn)變換就分解成了兩個(gè)N/2點(diǎn)變換,照這樣可繼續(xù)分解下去。這就是庫利-圖基快速傅里葉變換算法的基本原理。此時(shí),我們已經(jīng)不難分析出此時(shí)算法的時(shí)間復(fù)雜度將為O(NlogN)。
ok,沒想到,本文之中還是出現(xiàn)了這么多的公式(沒辦法,搞這個(gè)FFT就是個(gè)純數(shù)學(xué)活兒。之前買的一本小波與傅里葉分析基礎(chǔ)也是如此,就是一本數(shù)學(xué)公式和定理的聚集書。不過,能看懂更好,實(shí)在無法弄懂也只權(quán)當(dāng)做個(gè)稍稍了解)。
好了,下面,咱們來簡單理解下FFT中的蝶形算法,本文將告結(jié)束。如下圖所示:
有人這樣解釋蝶形算法:對(duì)于N(2的x次方)點(diǎn)的離散信號(hào),把它按索引位置分成兩個(gè)序列,分別為0,2,4,....,2K(記為A)和1,3,5,....,2K-1(記為B),由傅立葉變換可以推出N點(diǎn)的傅立葉變換前半部分的結(jié)果為A+B*旋轉(zhuǎn)因子,后半部分的結(jié)果為A-B*旋轉(zhuǎn)因子。于是求N點(diǎn)的傅立葉變換就變成分別求兩個(gè)N/2點(diǎn)序列的傅立葉變換,對(duì)每一個(gè)N/2點(diǎn)的序列,遞歸前面的步驟,直到只有兩點(diǎn)的序列,就只變成簡單的加減關(guān)系了。把這些點(diǎn)的加減關(guān)系用線連接,看上去就是個(gè)蝶形。ok,更多可參考算法導(dǎo)論第30章。
舉一個(gè)例子,我們知道,4X4的矩陣運(yùn)算如果按常規(guī)算法的話仍要進(jìn)行64次乘法運(yùn)算和48次加法,這將耗費(fèi)較多的時(shí)間,于是在H.264中,有一種改進(jìn)的算法(蝶形算法)可以減少運(yùn)算的次數(shù)。這種矩陣運(yùn)算算法構(gòu)造非常巧妙,利用構(gòu)造的矩陣的整數(shù)性質(zhì)和對(duì)稱性,可完全將乘法運(yùn)算轉(zhuǎn)化為加法運(yùn)算。如下圖所示:
下面的代碼來自一本數(shù)字圖像處理的書上的源代碼:
待我日后繼續(xù)斟酌,補(bǔ)充。完。
updated:關(guān)于快速傅立葉變換(FFT)的C++實(shí)現(xiàn)與Matlab實(shí)驗(yàn),這里有一篇不錯(cuò)的文章,讀者可以看看:http://blog.csdn.net/rappy/article/details/1700829。2012.03.03。
聯(lián)系客服