相信不少人和我一樣,第一次看到傅里葉變換是在算法書上實現(xiàn)快速高精度乘法的章節(jié),可是又看也看不懂,百度之后更加云里霧里.

今天,我要試圖用簡單但不一定正確的理解,探討快速傅里葉變換(FFT)和高精度乘法之間的關(guān)系.

傅里葉級數(shù):

在討論FFT之間,我們要說清楚一下各種傅里葉變換.

傅里葉變換最早要追朔到傅里葉級數(shù),傅里葉級數(shù)其實和冪級數(shù)是同一個玩意.

冪級數(shù)是說,我使用x^0,x^1,x^2,x^3....加上不同的系數(shù),可以表達世界上任何一個函數(shù)(夸張手法).

而傅里葉級數(shù)則說,我使用sin0x,cos0x,sinx,cosx,sin2x,cos2x,sin3x,cos3x.....加上不同的系數(shù),可以表達世界上任何一個周期函數(shù)(甚至非周期).

于是乎,傅里葉級數(shù)的系數(shù),和某個周期函數(shù)式,形成了對應(yīng)關(guān)系.

我知道傅里葉級數(shù)的所有系數(shù),就能求出對應(yīng)的函數(shù)式.

我知道某個函數(shù)式,我可以求出傅里葉級數(shù)的所有系數(shù).

這是高等數(shù)學上學到的.

離散傅里葉變換:

什么是離散傅里葉變換(DFT)?我用一個非常簡單的例子來說明.
問題1: 已知y=f(x)=2+3x+x^2,求x=0,1,2時y的值. 問題2: 已知一個二次函數(shù)f(x)經(jīng)過點(0,2),(1,6)(2,12),求f(x)的函數(shù)表達式.
相信這兩個題目,99%的人都會做.但是,我們不是要做這兩個題目,而是要尋找它們的奧妙. 在問題1,我們知道f(x)的三個系數(shù),求不同的x時y的值. 而問題2則反過來,我們知道三個不同的x時y的值,反過來求系數(shù). 顯然f(x)=2+3x+x^2和二次函數(shù)f(x)經(jīng)過(0,2),(1,6)(2,12),是在表達同一個函數(shù),只是用了不同的形式! 沒錯,離散傅里葉變換就是這樣的,函數(shù)式和傅里葉系數(shù)是周期函數(shù)的兩種不同表達形式! 而這里,多項式系數(shù)和點坐標是多項式函數(shù)的兩種不同表達形式. 我們知道多項式系數(shù),可以求點坐標.我們知道點坐標,可以反過來求多項式系數(shù).我們稱這種變換為離散傅里葉變換.
有一點要注意的是,要想通過點坐標求出所有系數(shù)是有要求的. 如果是一個二次函數(shù),它有3個系數(shù),那么我們需要三個點坐標的信息,并且這3個點不能重合,即x坐標不能一樣. 推廣的話,n次函數(shù),有n+1個系數(shù),需要n+1個點的信息. 這在線性代數(shù)里能得到有關(guān)的敘述和嚴格的證明.

時域和頻域:

傅里葉變換有兩個很重要的術(shù)語,時域和頻域. 其實,在DFT里,簡單地說,時域就是點坐標,頻域就是系數(shù). 這么說好像不容易理解,畢竟時域和頻域是在信號與系統(tǒng)那方面命名的,用在這里有點奇怪. 不過,我們就這樣記住好了.時域是點坐標,頻域是系數(shù). 要想理解時域和頻域的命名含義,可以看有關(guān)傅里葉變換的書籍.

快速傅里葉變換:

快速傅里葉變換FFT又是什么呢?我們先看離散傅里葉變換DFT. 既然說FFT是個算法,談到算法我們當然要談時間復(fù)雜度.要實現(xiàn)DFT,時間復(fù)雜度是很高的. 比如從頻域到時域(知道系數(shù),給出不同的x坐標,求y坐標),我們有O(n)個坐標要算,每個坐標是O(n)個項,復(fù)雜度就是O(n2)了. 反過來,從時域到頻域(知道坐標,求系數(shù)),我們要解一個n元一次方程組.用線性代數(shù)的高斯消元法,復(fù)雜度是O(n3). 那么,FFT就是來降低這個復(fù)雜度的,它能把DFT的這兩個轉(zhuǎn)換的復(fù)雜度降到O(nlogn). 哇,真的有這么厲害?! 是的,只是這里有一個很大的限制. 什么限制?下面再說.

多項式乘法與高精度乘法:

本文的第二個主角,高精度乘法終于出場了. 實現(xiàn)高精度乘法,其實和實現(xiàn)多項式乘法是幾乎一樣的. 這話怎說?其實,我們用某個進制去表達一個數(shù)字,恰恰和多項式的表達一模一樣. 比如15386,它其實是1*10^4+5*10^3+3*10^2+8*10^1+6*10^0. 這和函數(shù)f(x)=x^4+5*x^3+3*x^2+8*x^1+6是一個套路. 那么,兩個數(shù)字的乘法,其實就是兩個多項式的乘法. 例如12*43,我們可以理解為(x+2)(4x+3). (x+2)(4x+3)=4x^2+11x+6,所以12*43=4*10^2+11*10+6.
稍有不同的是,數(shù)字乘法要進位,所以(4,11,6)要變成(5,1,6),也就是12*43=516.
好了,我們已經(jīng)明白了多項式乘法和高精度乘法是蛇鼠一窩了. 接下來,我們看看多項式乘法的算法復(fù)雜度.分析很簡單,兩個for循環(huán),O(n2). 而FFT告訴你,我能O(nlogn)幫你搞定!

FFT實現(xiàn)多項式乘法:

FFT如何實現(xiàn)多項式乘法,我們先探討多項式乘法的本質(zhì).

多項式乘法,其實就是給出你兩個多項式的頻域,求它們的積的頻域!

我們知道兩個多項式的表達式,其實就是知道頻域(系數(shù)),而我們要求的是它們的積的表達式,也就是求頻域.
我們把上面的結(jié)論收集起來,說明如何用FFT實現(xiàn)O(nlogn)多項式乘法. 1.FFT能O(nlogn)實現(xiàn)頻域到時域 2.FFT能O(nlong)實現(xiàn)時域到頻域 3.我們知道兩個多項式的頻域,要求它們的積的頻域.
只有這三個結(jié)論是不夠的,我們還需要一個很白癡但你可能想不到的結(jié)論: 4.我們知道兩個多項式的時域,能O(n)時間求出它們的積的時域.
這個結(jié)論看似很難,但講清楚之后你會發(fā)現(xiàn)很白癡. 知道兩個多項式的時域,就是知道了f(x)經(jīng)過(x1,f1),(x2,f2)...(xn,fn),g(x)經(jīng)過(x1,g1),(x2,g2)....(xn,gn) 設(shè)k(x)=f(x)*g(x),那么k(x)在x1,x2...xn上的y坐標是多少?一想就知道是f1*g1,f2*g2....fn*gn!!!因為k(xi)=f(xi)*g(xi)啊!! 也就是說k(x)會經(jīng)過(x1,f1*g1),(x2,f2*g2)....(xn,fn*gn).這些點,就是k(x)的時域啊! 這里的復(fù)雜度不用說都知道是O(n),原來"時域相乘"是這么簡單的事情,比"頻域相乘"簡單多了.
現(xiàn)在我們可以用FFT實現(xiàn)O(nlogn)多項式乘法了! 1.用FFT在O(nlogn)時間,把f(x)和g(x)的頻域轉(zhuǎn)變?yōu)闀r域 2.用O(n)時間,把f(x)和g(x)的時域變成k(x)的時域 3.用FFT在O(nlogn)時間,把k(x)的時域轉(zhuǎn)變?yōu)轭l域

FFT算法:

前面我們說了一大堆,目的是說明為什么FFT能實現(xiàn)nlogn多項式乘法.
現(xiàn)在,我們剩下的問題,也就是最核心的問題是,如何實現(xiàn)FFT算法?

FFT算法包括頻域到時域,以及時域到頻域.
我們在這里著重討論頻域到時域,也就是上面那個圖的第一個箭頭.
現(xiàn)在我們知道的是兩個函數(shù)的頻域(n個多項式系數(shù)),我們要求它們的時域(n個點坐標).
還有一點是,這個n一定要是2的冪.如果不是,你可以強制改成是.
比如n=6,你可以加上0*n^7+0*n^6,強行讓它有8個項.

首先有一個問題,如何選擇時域的一組x坐標呢?
我們知道,時域的n個點并不是確定的,它們只需要保證x坐標互不相同就能構(gòu)成一個時域.
那么,我們是不是可以隨便找n個x坐標就可以了呢?像x=1,2,3...n這樣,行嗎?
答案是不行!前文說道FFT轉(zhuǎn)換有一個限制,這個限制就在這里:我們不能任意選擇x坐標!

你要得到x=1,2...對應(yīng)的y坐標,抱歉,FFT告訴你我做不到,你要求只能老老實實O(n2).
那么,什么樣的一組x坐標,FFT能O(nlogn)求出對應(yīng)的y坐標呢?
答案是x^n=1這個方程的所有根.

(下面所說的需要比較多的數(shù)學知識,請看不懂的同志自行百度惡補,這里有一篇參考文章:http://wenku.baidu.com/link?url=RWHc2uPEZWdbo0TTVM3z320hJidSR5hk-6YcYQbG6n9n2WzrDKcx1E4r_XS0ZtIqkLjgxODMq_Q8GlK-wGJgZCyHsd1aPhFjBC9Ute2yB0G)
哈?!這個方程的解不就只有1,如果n是偶數(shù)就是1和-1.我可是要n個x坐標喔!
注意了,我這里是說x^n=1的所有復(fù)數(shù)根,而不只是實數(shù)跟.
可能有些同學復(fù)數(shù)學得少,在這里就困惑了.x^n=1這方程一定有n個復(fù)數(shù)根嗎?
答案是一定的,一元n次方程有且只有n個復(fù)數(shù)根.這是代數(shù)基本定理.
并且,我還可以告訴你這n個復(fù)數(shù)根分別是什么,使用復(fù)數(shù)開方公式就可以了.

對于一個復(fù)數(shù)z,我們使用其三角函數(shù)形式表達,z=r(cosθ+isinθ)
那么z開n次方的所有根就是下圖:


啥?k∈Z,那這里不就有無限個復(fù)數(shù)?
并不是哦,注意k是在cos和sin里面的,是有周期的.
并且很容易知道,周期就是n,也就是說k=0和k=n的結(jié)果是一樣的.
就看式子可能還是很抽象,我們把這n個根畫在復(fù)平面上,就一目了然了.


這是算法導(dǎo)論上的圖片,它定義了一個符號.(下面我使用w[k,n]去表達這個符號)
w[i,n]其實是表示x^n=1這個方程的所有根中,令k=i的那一個(k就是剛才那個式子的k)
事實上w[k,n]的k是k次方的意思,而w[n]則是e^(2πi/n).
但這兩個定義是等價的,前者是我自己的理解,后者是算法導(dǎo)論上的定義.(詳情參考算法導(dǎo)論)
在這里,可以很清楚地看到x^8=1中真的有8個復(fù)數(shù)根,并且那個周期真的是8.

好,我們得到了n個x坐標,但是是復(fù)數(shù)的.
有人可能會問,這個x坐標是復(fù)數(shù)的,怎么對應(yīng)在xy平面上的某個點啊?
答案當然是不能對應(yīng),正確來說是不能在xy平面上對應(yīng).
這里的x坐標和y坐標,其實已經(jīng)是一個四維面上的兩條軸了,不能單單用xy平面去理解了.
不過,這里FFT并不關(guān)心你是實數(shù)還是復(fù)數(shù),它就是根據(jù)x坐標,算出y坐標而已.

那么,FFT算法要怎么O(nlogn)把這n個x坐標轉(zhuǎn)換成對應(yīng)的y坐標呢?
答案是分治,FFT算法是一個分治的算法,分而治之.
那么,兩個關(guān)鍵的問題就是,如何"分"?如何"治"?

考慮如何分,就要先說清楚原問題和子問題是什么,我們現(xiàn)在說清楚.
FFT算法,是已知一個n-1次的多項式函數(shù),求出n個x坐標對應(yīng)的y坐標,這n個x坐標是x^n=1的n個復(fù)數(shù)根.(注意是n-1次,不是n次,因為n-1次多項式有n項)
舉個例子,現(xiàn)在我有一個(8-1)次的多項式,我要求出x^8=1的8個根作為x坐標,對應(yīng)的y坐標.
那么子問題是,我有一個(4-1)次的多項式,我要求出x^4=1的4個根作為x坐標,對應(yīng)的y坐標.

那么,如何從原問題轉(zhuǎn)變成子問題?
首先,我們需要幾個挺容易理解的公式.
公式1:w[k,2n]^2=w[k,n] 公式2:w[dk,dn]=w[k,n] 公式3:w[j+k,n]=w[(j+k)%n,n] 公式4:w[-k,n]=1/w[k,n]
現(xiàn)在,以n=4為例子,y=a0+a1*x+a2*x^2+a3*x^3.
我們的問題是,要求出x^4=1的4個根作為x坐標,對應(yīng)的y坐標.
我們按照x的指數(shù)的奇偶性,把y分成兩部分ya和yb.
ya=a0+a2*x^2
yb=a1*x+a3*x^3=x(a1+a3*x^2)
然后我們定義y1,y2
y1=a0+a2*x
y2=a1+a3*x

然后,我們得到了子問題,我們已知y1和y2這兩個(2-1)次函數(shù).
我要求出x^2=1的2個根作為x坐標,對應(yīng)的y1和y2.

為什么這樣就是子問題了,或者說為什么要選這樣的y1和y2作為子問題的多項式?
因為首先它符合子問題的要求,其次是我們有了這兩個子問題的答案,能給出原問題的答案.
也就是,它能幫我們實現(xiàn)分而治之的"治".

怎么治?
我們從子問題得到的答案是y1(w[0,2]),y1(w[1,2]),y2(w[0,2])和y2(w[1,2]).
而我們要求的答案是y(w[0,4]),y(w[1,4]),y(w[2,4]),y(w[3,4]).
在這里,上面的幾條公式就有用了.

我們看y和y1y2的關(guān)系.
y(x)=y1(x^2)+x*y2(x^2)
那么,我們把w[0,4]代入x里,得到:
y(w[0,4])=y1(w[0,4]^2)+w[0,4]*y2(w[0,4]^2)
根據(jù)公式1,w[0,4]^2=w[0,2],所以:
y(w[0,4])=y1(w[0,2])+w[0,4]*y2(w[0,2])
我們也代入w[1,4],得到:
y(w[1,4])=y1(w[1,2])+w[1,4]*y2(w[1,2])
那w[2,4]呢?
y(w[2,4])=y1(w[2,2])+w[2,4]*y2(w[2,2])
w[2,2]是啥,根據(jù)公式三的周期性質(zhì),w[0,2]=w[2,2],所以:
y(w[2,4])=y1(w[0,2])+w[2,4]*y2(w[0,2])
同樣道理寫出w[3,4]:
y(w[3,4])=y1(w[1,2])+w[3,4]*y2(w[1,2])

寫在一起就是:
y(w[0,4])=y1(w[0,2])+w[0,4]*y2(w[0,2])
y(w[1,4])=y1(w[1,2])+w[1,4]*y2(w[1,2])
y(w[2,4])=y1(w[2,2])+w[2,4]*y2(w[2,2])
y(w[3,4])=y1(w[1,2])+w[3,4]*y2(w[1,2])

有沒有發(fā)現(xiàn),求y(w[0,4]),y(w[1,4]),y(w[2,4]),y(w[3,4])的式子中,除了w[0,4],w[1,4],w[2,4],w[3,4],其他項都是子問題返回的答案,y1(w[0,2]),y1(w[1,2]),y2(w[0,2])和y2(w[1,2]).
那剩下的問題是怎么求w[0,4],w[1,4],w[2,4],w[3,4]了.
剛開始,我以為這些根是用某種特殊形式表現(xiàn)的.
結(jié)果發(fā)現(xiàn),它就直接用浮點復(fù)數(shù)表現(xiàn)了!
也就是說,即使你想實現(xiàn)整數(shù)多項式相乘,使用FFT得到的結(jié)果是一個浮點數(shù)多項式,你需要四舍五入去得到正確的整數(shù)!

那w[0,4],w[1,4],w[2,4],w[3,4]到底怎么求?你可以直接用開方的那個式子!
那式子中要開方的復(fù)數(shù)是1,所以r=1,θ=0即可.因此:
w[k,n]=cos(2kπ/n)+isin(2kπ/n)
將k和n都代進去算的話就得到:
w[0,4]=1
w[1,4]=i
w[2,4]=-1
w[3,4]=-i
這里n=4的情況下,這些復(fù)數(shù)的實部和虛部都是整數(shù).
但從8開始就不是整數(shù)的了,甚至大部分都是無理數(shù),因此只能用浮點數(shù)去表示.

還有一個小問題是遞歸的終點,終點就是n=1的時候.
這個時候,我們要求的是w[0,1]其實就是1作為x坐標,對應(yīng)的y坐標.
在這里我們的多項式只有一項,也就是只有常數(shù),因此我們返回這個常數(shù)就行.

到這里,我們真的實現(xiàn)了從頻域到時域的轉(zhuǎn)換.
不難知道這個算法復(fù)雜度是O(nlogn).(f(n)=2*f(n/2)+O(n),和歸并排序類似.)
我們用python去實現(xiàn)一下,這個優(yōu)美的FFT算法.

代碼:

運行結(jié)果:

代碼解釋:

多項式我們使用列表來儲存,從前到后分別是0次項,1次項,2次項...
子問題返回的n個y坐標也是用列表儲存,從前到后分別是y0,y1,y2...

python中復(fù)數(shù)能直接表示和運算,無需特殊的類.
并且,虛數(shù)i在python中用j表示.

第6,7行是遞歸終點的特殊處理.
第9到13行是構(gòu)建y1和y2.
第15,16行是遞歸調(diào)用,得到子問題答案.
第19行到第23行是計算并返回原問題答案.
第20行和第21行的計算公式看起來比較復(fù)雜,但根據(jù)前文推導(dǎo)就發(fā)現(xiàn)也不是太難.

我們的樣例是f(x)=2+3x+x^2+2x^3,也就是頻域是[2,3,1,2]
從運行結(jié)果可以看到,FFT得到的時域是一堆浮點數(shù).
我們整理后得到的時域是[8,1+i,-2,1-i]
如果你將x^4=1的四個根1,i,-1,-i代入多項式,得到的結(jié)果也就是這四個值.
這說明,算法正確!

我們已經(jīng)能實現(xiàn)從頻域到時域了,那我們要怎么做才能實現(xiàn)時域到頻域?
數(shù)學家發(fā)現(xiàn),從頻域到時域和從時域到頻域的實現(xiàn)幾乎一樣!
這要使用線性代數(shù)的知識去解釋.

從頻域到時域,其實就是頻域的n個系數(shù)組成一個列,然后乘一個n*n矩陣,得到另一個列,而這個列就是y坐標組成的列.如下圖(出自算法導(dǎo)論):


為什么中間的n*n矩陣是這個樣子?你用矩陣乘法算一算某個行乘列a就懂了,這其實就是把w[k,n]代進了多項式后的式子.
為什么我們要把從頻域到時域理解成矩陣乘法呢?因為矩陣乘法有一個特殊的性質(zhì).
如果y=M*a,矩陣M是可逆的,那么a=M^(-1)*y.
它的意思就是,如果M是可逆矩陣,那么我們拿M的逆去乘時域y,就能得到頻域a.
也就是,它能實現(xiàn)從時域到頻域.

那M可逆嗎?這個矩陣其實是有名字的,叫范德蒙矩陣.
嗯,在線性代數(shù)里可以證明,它是可逆的.也就是說,有戲!
那M的逆是什么呢?這個線性代數(shù)也幫我們算好了,如下圖:


有了逆矩陣有什么用呢?矩陣乘法復(fù)雜度可是O(n3)的.
誰說用矩陣乘法了?我們還是用FFT.
我們把得到的時域y理解為新的一組系數(shù),把要求的頻域a理解為新的一組y坐標,然后新選取的一組x坐標是w[-1,n],w[-2,n]...
你就會發(fā)現(xiàn),這樣FFT轉(zhuǎn)換的含義剛好和這個矩陣乘法式子吻合.并且選取的這組x坐標,也能像剛才那組一樣分而治之(你可以像上面一樣試一下.)
也就是說,我們可以依樣畫葫蘆,使用FFT實現(xiàn)從時域到頻域,只是選取的x坐標換了,并且最后得到的結(jié)果要除以一個n.
w[-1,n],w[-2,n]怎么求...公式四w[-k,n]=1/w[k,n],完.

代碼:

運行結(jié)果:

代碼解釋:

我們在原來的FFT函數(shù)里加上inverse參數(shù),inverse為true的話,表示我們的FFT要實現(xiàn)從時域到頻域.
然后在第21行加上一句,如果inverse為true則w=1/w,即現(xiàn)在的w表示w[-i,n]而不是w[i,n].

樣例還是使用原來的,我們先將頻域轉(zhuǎn)為時域得到y(tǒng),再從時域轉(zhuǎn)為頻域得到new_fx.
由于FFT函數(shù)里沒有幫我們將結(jié)果除以n,所以我們在第30行自己實現(xiàn)了.

整理一下運行結(jié)果后可以看到,new_fx=[2,3,1,2].
也就是說,new_fx=fx!我們的FFT轉(zhuǎn)換是對的!
我們實現(xiàn)了從頻域到時域,也實現(xiàn)了從時域到頻域.

FFT算法實現(xiàn)多項式乘法:

最后我們來用FFT算法實現(xiàn)多項式乘法,在這里我們要特殊處理一些東西.
一個是多項式的長度,不是2^k要變成2^k.
另一個是積的時域的問題,假如兩個要乘的多項式都是有n項,我們在FFT轉(zhuǎn)換后得到n個點.
這n個點的y坐標相乘得到目標多項式的n個點,但是目標多項式理論上最多是有2n-1項的(如果是整數(shù)乘法則可以有2n位).
所以,這n個點會造成信息的丟失,它并不能構(gòu)成目標的時域.
怎么辦?我們一開始求2n個點就好了! 也就是把要乘的那兩個多項式的長度增長到2n,增長的部分的系數(shù)全部填0,OK~!
那假如兩個要乘的多項式的項數(shù)不一樣呢?簡單,把短的拉成和長的一樣.
(PS:這里的項數(shù)的意思其實是最高次數(shù)+1,例如x^2+1我們要說成3項.)

代碼:

運行結(jié)果:

代碼解釋:

上圖代碼調(diào)用的fft函數(shù)就是前一份代碼里的fft函數(shù). 代碼里的27行到34行,都是在處理長度問題,包括長度統(tǒng)一,找符合條件的最小的2的冪,以及擴大兩倍. 36,37行是第一個箭頭,實現(xiàn)頻域到時域. 39行則是第二個箭頭,求出目標多項式的時域. 41,42行是第三個箭頭,實現(xiàn)時域到頻域.
我們用(1+2x+3x^2)*(4+5x)作為樣例,手算一下答案是4+13x+22x^2+15x^3. 那運行結(jié)果呢,一堆很難看的浮點數(shù). 我們整理一下,得出結(jié)果就是4+13x+22x^2+15x^3. 雖然結(jié)果有8項,但很容易發(fā)現(xiàn)后面4項都是0.

FFT實現(xiàn)多項式乘法,在實際的操作上有許多種版本,我這里只是其中的一種,并且是效率不高的一種. 但是,只要思路掌握了,無論什么樣的版本都是一個套路.

總結(jié)

以上是生活随笔為你收集整理的FFT:快速傅里叶变换与高精度乘法的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。