FFT:快速傅立葉變換與高精度乘法
阿新 • • 發佈：2019-02-07
相信不少人和我一樣,第一次看到傅立葉變換是在演算法書上實現快速高精度乘法的章節,可是又看也看不懂,百度之後更加雲裡霧裡.
今天,我要試圖用簡單但不一定正確的理解,探討快速傅立葉變換(FFT)和高精度乘法之間的關係.
傅立葉級數:
在討論FFT之間,我們要說清楚一下各種傅立葉變換.
傅立葉變換最早要追朔到傅立葉級數,傅立葉級數其實和冪級數是同一個玩意.
冪級數是說,我使用x^0,x^1,x^2,x^3....加上不同的係數,可以表達世界上任何一個函式(誇張手法).
而傅立葉級數則說,我使用sin0x,cos0x,sinx,cosx,sin2x,cos2x,sin3x,cos3x.....加上不同的係數,可以表達世界上任何一個周期函式(甚至非週期).
於是乎,傅立葉級數的係數,和某個周期函式式,形成了對應關係.
我知道傅立葉級數的所有係數,就能求出對應的函式式.
我知道某個函式式,我可以求出傅立葉級數的所有係數.
這是高等數學上學到的.
離散傅立葉變換:
什麼是離散傅立葉變換(DFT)?我用一個非常簡單的例子來說明. 問題1: 已知y=f(x)=2+3x+x^2,求x=0,1,2時y的值. 問題2: 已知一個二次函式f(x)經過點(0,2),(1,6)(2,12),求f(x)的函式表示式. 相信這兩個題目,99%的人都會做.但是,我們不是要做這兩個題目,而是要尋找它們的奧妙. 在問題1,我們知道f(x)的三個係數,求不同的x時y的值. 而問題2則反過來,我們知道三個不同的x時y的值,反過來求係數. 顯然f(x)=2+3x+x^2和二次函式f(x)經過(0,2),(1,6)(2,12),是在表達同一個函式,只是用了不同的形式
時域和頻域:
傅立葉變換有兩個很重要的術語,時域和頻域.
其實,在DFT裡,簡單地說,時域就是點座標,頻域就是係數.
這麼說好像不容易理解,畢竟時域和頻域是在訊號與系統那方面命名的,用在這裡有點奇怪.
不過,我們就這樣記住好了.時域是點座標,頻域是係數.
要想理解時域和頻域的命名含義,可以看有關傅立葉變換的書籍.
快速傅立葉變換:
快速傅立葉變換FFT又是什麼呢?我們先看離散傅立葉變換DFT.
既然說FFT是個演算法,談到演算法我們當然要談時間複雜度.要實現DFT,時間複雜度是很高的.
比如從頻域到時域(知道係數,給出不同的x座標,求y座標),我們有O(n)個座標要算,每個座標是O(n)個項,複雜度就是O(n2)了.
反過來,從時域到頻域(知道座標,求係數),我們要解一個n元一次方程組.用線性代數的高斯消元法,複雜度是O(n3).
那麼,FFT就是來降低這個複雜度的,它能把DFT的這兩個轉換的複雜度降到O(nlogn).
哇,真的有這麼厲害?!
是的,只是這裡有一個很大的限制.
什麼限制?下面再說.
多項式乘法與高精度乘法:
本文的第二個主角,高精度乘法終於出場了. 實現高精度乘法,其實和實現多項式乘法是幾乎一樣的. 這話怎說?其實,我們用某個進位制去表達一個數字,恰恰和多項式的表達一模一樣. 比如15386,它其實是1*10^4+5*10^3+3*10^2+8*10^1+6*10^0. 這和函式f(x)=x^4+5*x^3+3*x^2+8*x^1+6是一個套路. 那麼,兩個數字的乘法,其實就是兩個多項式的乘法. 例如12*43,我們可以理解為(x+2)(4x+3). (x+2)(4x+3)=4x^2+11x+6,所以12*43=4*10^2+11*10+6.稍有不同的是,數字乘法要進位,所以(4,11,6)要變成(5,1,6),也就是12*43=516. 好了,我們已經明白了多項式乘法和高精度乘法是蛇鼠一窩了. 接下來,我們看看多項式乘法的演算法複雜度.分析很簡單,兩個for迴圈,O(n2). 而FFT告訴你,我能O(nlogn)幫你搞定!
FFT實現多項式乘法:
FFT如何實現多項式乘法,我們先探討多項式乘法的本質.多項式乘法,其實就是給出你兩個多項式的頻域,求它們的積的頻域!
我們知道兩個多項式的表示式,其實就是知道頻域(係數),而我們要求的是它們的積的表示式,也就是求頻域. 我們把上面的結論收集起來,說明如何用FFT實現O(nlogn)多項式乘法. 1.FFT能O(nlogn)實現頻域到時域 2.FFT能O(nlong)實現時域到頻域 3.我們知道兩個多項式的頻域,要求它們的積的頻域. 只有這三個結論是不夠的,我們還需要一個很白痴但你可能想不到的結論: 4.我們知道兩個多項式的時域,能O(n)時間求出它們的積的時域. 這個結論看似很難,但講清楚之後你會發現很白痴. 知道兩個多項式的時域,就是知道了f(x)經過(x1,f1),(x2,f2)...(xn,fn),g(x)經過(x1,g1),(x2,g2)....(xn,gn) 設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)會經過(x1,f1*g1),(x2,f2*g2)....(xn,fn*gn).這些點,就是k(x)的時域啊! 這裡的複雜度不用說都知道是O(n),原來"時域相乘"是這麼簡單的事情,比"頻域相乘"簡單多了. 現在我們可以用FFT實現O(nlogn)多項式乘法了! 1.用FFT在O(nlogn)時間,把f(x)和g(x)的頻域轉變為時域 2.用O(n)時間,把f(x)和g(x)的時域變成k(x)的時域 3.用FFT在O(nlogn)時間,把k(x)的時域轉變為頻域
FFT演算法:
前面我們說了一大堆,目的是說明為什麼FFT能實現nlogn多項式乘法.現在,我們剩下的問題,也就是最核心的問題是,如何實現FFT演算法?
FFT演算法包括頻域到時域,以及時域到頻域.
我們在這裡著重討論頻域到時域,也就是上面那個圖的第一個箭頭.
現在我們知道的是兩個函式的頻域(n個多項式係數),我們要求它們的時域(n個點座標).
還有一點是,這個n一定要是2的冪.如果不是,你可以強制改成是.
比如n=6,你可以加上0*n^7+0*n^6,強行讓它有8個項.
首先有一個問題,如何選擇時域的一組x座標呢?
我們知道,時域的n個點並不是確定的,它們只需要保證x座標互不相同就能構成一個時域.
那麼,我們是不是可以隨便找n個x座標就可以了呢?像x=1,2,3...n這樣,行嗎?
答案是不行!前文說道FFT轉換有一個限制,這個限制就在這裡:我們不能任意選擇x座標!
你要得到x=1,2...對應的y座標,抱歉,FFT告訴你我做不到,你要求只能老老實實O(n2).
那麼,什麼樣的一組x座標,FFT能O(nlogn)求出對應的y座標呢?
答案是x^n=1這個方程的所有根.
(下面所說的需要比較多的數學知識,請看不懂的同志自行百度惡補,這裡有一篇參考文章:http://wenku.baidu.com/link?url=RWHc2uPEZWdbo0TTVM3z320hJidSR5hk-6YcYQbG6n9n2WzrDKcx1E4r_XS0ZtIqkLjgxODMq_Q8GlK-wGJgZCyHsd1aPhFjBC9Ute2yB0G)
哈?!這個方程的解不就只有1,如果n是偶數就是1和-1.我可是要n個x座標喔!
注意了,我這裡是說x^n=1的所有複數根,而不只是實數跟.
可能有些同學複數學得少,在這裡就困惑了.x^n=1這方程一定有n個複數根嗎?
答案是一定的,一元n次方程有且只有n個複數根.這是代數基本定理.
並且,我還可以告訴你這n個複數根分別是什麼,使用複數開方公式就可以了.
對於一個複數z,我們使用其三角函式形式表達,z=r(cosθ+isinθ)
那麼z開n次方的所有根就是下圖:
啥?k∈Z,那這裡不就有無限個複數?
並不是哦,注意k是在cos和sin裡面的,是有周期的.
並且很容易知道,週期就是n,也就是說k=0和k=n的結果是一樣的.
就看式子可能還是很抽象,我們把這n個根畫在複平面上,就一目瞭然了.
這是演算法導論上的圖片,它定義了一個符號
.(下面我使用w[k,n]去表達這個符號)w[i,n]其實是表示x^n=1這個方程的所有根中,令k=i的那一個(k就是剛才那個式子的k)
事實上w[k,n]的k是k次方的意思,而w[n]則是e^(2πi/n).
但這兩個定義是等價的,前者是我自己的理解,後者是演算法導論上的定義.(詳情參考演算法導論)
在這裡,可以很清楚地看到x^8=1中真的有8個複數根,並且那個週期真的是8.
好,我們得到了n個x座標,但是是複數的.
有人可能會問,這個x座標是複數的,怎麼對應在xy平面上的某個點啊?
答案當然是不能對應,正確來說是不能在xy平面上對應.
這裡的x座標和y座標,其實已經是一個四維面上的兩條軸了,不能單單用xy平面去理解了.
不過,這裡FFT並不關心你是實數還是複數,它就是根據x座標,算出y座標而已.
那麼,FFT演算法要怎麼O(nlogn)把這n個x座標轉換成對應的y座標呢?
答案是分治,FFT演算法是一個分治的演算法,分而治之.
那麼,兩個關鍵的問題就是,如何"分"?如何"治"?
考慮如何分,就要先說清楚原問題和子問題是什麼,我們現在說清楚.
FFT演算法,是已知一個n-1次的多項式函式,求出n個x座標對應的y座標,這n個x座標是x^n=1的n個複數根.(注意是n-1次,不是n次,因為n-1次多項式有n項)
舉個例子,現在我有一個(8-1)次的多項式,我要求出x^8=1的8個根作為x座標,對應的y座標.
那麼子問題是,我有一個(4-1)次的多項式,我要求出x^4=1的4個根作為x座標,對應的y座標.
那麼,如何從原問題轉變成子問題?
首先,我們需要幾個挺容易理解的公式.
公式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]
現在,以n=4為例子,y=a0+a1*x+a2*x^2+a3*x^3.
我們的問題是,要求出x^4=1的4個根作為x座標,對應的y座標.
我們按照x的指數的奇偶性,把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)次函式.
我要求出x^2=1的2個根作為x座標,對應的y1和y2.
為什麼這樣就是子問題了,或者說為什麼要選這樣的y1和y2作為子問題的多項式?
因為首先它符合子問題的要求,其次是我們有了這兩個子問題的答案,能給出原問題的答案.
也就是,它能幫我們實現分而治之的"治".
怎麼治?
我們從子問題得到的答案是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的關係.
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)
根據公式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]是啥,根據公式三的週期性質,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])
有沒有發現,求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]了.
剛開始,我以為這些根是用某種特殊形式表現的.
結果發現,它就直接用浮點複數表現了!
也就是說,即使你想實現整數多項式相乘,使用FFT得到的結果是一個浮點數多項式,你需要四捨五入去得到正確的整數!
那w[0,4],w[1,4],w[2,4],w[3,4]到底怎麼求?你可以直接用開方的那個式子!
那式子中要開方的複數是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的情況下,這些複數的實部和虛部都是整數.
但從8開始就不是整數的了,甚至大部分都是無理數,因此只能用浮點數去表示.
還有一個小問題是遞迴的終點,終點就是n=1的時候.
這個時候,我們要求的是w[0,1]其實就是1作為x座標,對應的y座標.
在這裡我們的多項式只有一項,也就是隻有常數,因此我們返回這個常數就行.
到這裡,我們真的實現了從頻域到時域的轉換.
不難知道這個演算法複雜度是O(nlogn).(f(n)=2*f(n/2)+O(n),和歸併排序類似.)
我們用python去實現一下,這個優美的FFT演算法.
程式碼:
執行結果:
程式碼解釋:
多項式我們使用列表來儲存,從前到後分別是0次項,1次項,2次項...子問題返回的n個y座標也是用列表儲存,從前到後分別是y0,y1,y2...
python中複數能直接表示和運算,無需特殊的類.
並且,虛數i在python中用j表示.
第6,7行是遞迴終點的特殊處理.
第9到13行是構建y1和y2.
第15,16行是遞迴呼叫,得到子問題答案.
第19行到第23行是計算並返回原問題答案.
第20行和第21行的計算公式看起來比較複雜,但根據前文推導就發現也不是太難.
我們的樣例是f(x)=2+3x+x^2+2x^3,也就是頻域是[2,3,1,2]
從執行結果可以看到,FFT得到的時域是一堆浮點數.
我們整理後得到的時域是[8,1+i,-2,1-i]
如果你將x^4=1的四個根1,i,-1,-i代入多項式,得到的結果也就是這四個值.
這說明,演算法正確!
我們已經能實現從頻域到時域了,那我們要怎麼做才能實現時域到頻域?
數學家發現,從頻域到時域和從時域到頻域的實現幾乎一樣!
這要使用線性代數的知識去解釋.
從頻域到時域,其實就是頻域的n個係數組成一個列,然後乘一個n*n矩陣,得到另一個列,而這個列就是y座標組成的列.如下圖(出自演算法導論):
為什麼中間的n*n矩陣是這個樣子?你用矩陣乘法算一算某個行乘列a就懂了,這其實就是把w[k,n]代進了多項式後的式子.
為什麼我們要把從頻域到時域理解成矩陣乘法呢?因為矩陣乘法有一個特殊的性質.
如果y=M*a,矩陣M是可逆的,那麼a=M^(-1)*y.
它的意思就是,如果M是可逆矩陣,那麼我們拿M的逆去乘時域y,就能得到頻域a.
也就是,它能實現從時域到頻域.
那M可逆嗎?這個矩陣其實是有名字的,叫範德蒙矩陣.
嗯,線上性代數裡可以證明,它是可逆的.也就是說,有戲!
那M的逆是什麼呢?這個線性代數也幫我們算好了,如下圖:
有了逆矩陣有什麼用呢?矩陣乘法複雜度可是O(n3)的.
誰說用矩陣乘法了?我們還是用FFT.
我們把得到的時域y理解為新的一組係數,把要求的頻域a理解為新的一組y座標,然後新選取的一組x座標是w[-1,n],w[-2,n]...
你就會發現,這樣FFT轉換的含義剛好和這個矩陣乘法式子吻合.並且選取的這組x座標,也能像剛才那組一樣分而治之(你可以像上面一樣試一下.)
也就是說,我們可以依樣畫葫蘆,使用FFT實現從時域到頻域,只是選取的x座標換了,並且最後得到的結果要除以一個n.
w[-1,n],w[-2,n]怎麼求...公式四w[-k,n]=1/w[k,n],完.
程式碼:
執行結果:
程式碼解釋:
我們在原來的FFT函式里加上inverse引數,inverse為true的話,表示我們的FFT要實現從時域到頻域.然後在第21行加上一句,如果inverse為true則w=1/w,即現在的w表示w[-i,n]而不是w[i,n].
樣例還是使用原來的,我們先將頻域轉為時域得到y,再從時域轉為頻域得到new_fx.
由於FFT函式裡沒有幫我們將結果除以n,所以我們在第30行自己實現了.
整理一下執行結果後可以看到,new_fx=[2,3,1,2].
也就是說,new_fx=fx!我們的FFT轉換是對的!
我們實現了從頻域到時域,也實現了從時域到頻域.
FFT演算法實現多項式乘法:
最後我們來用FFT演算法實現多項式乘法,在這裡我們要特殊處理一些東西.一個是多項式的長度,不是2^k要變成2^k.
另一個是積的時域的問題,假如兩個要乘的多項式都是有n項,我們在FFT轉換後得到n個點.
這n個點的y座標相乘得到目標多項式的n個點,但是目標多項式理論上最多是有2n-1項的(如果是整數乘法則可以有2n位).
所以,這n個點會造成資訊的丟失,它並不能構成目標的時域.
怎麼辦?我們一開始求2n個點就好了! 也就是把要乘的那兩個多項式的長度增長到2n,增長的部分的係數全部填0,OK~!
那假如兩個要乘的多項式的項數不一樣呢?簡單,把短的拉成和長的一樣.
(PS:這裡的項數的意思其實是最高次數+1,例如x^2+1我們要說成3項.)
程式碼:
執行結果:
程式碼解釋:
上圖程式碼呼叫的fft函式就是前一份程式碼裡的fft函式.
程式碼裡的27行到34行,都是在處理長度問題,包括長度統一,找符合條件的最小的2的冪,以及擴大兩倍.
36,37行是第一個箭頭,實現頻域到時域.
39行則是第二個箭頭,求出目標多項式的時域.
41,42行是第三個箭頭,實現時域到頻域.
我們用(1+2x+3x^2)*(4+5x)作為樣例,手算一下答案是4+13x+22x^2+15x^3.
那執行結果呢,一堆很難看的浮點數.
我們整理一下,得出結果就是4+13x+22x^2+15x^3.
雖然結果有8項,但很容易發現後面4項都是0.FFT實現多項式乘法,在實際的操作上有許多種版本,我這裡只是其中的一種,並且是效率不高的一種. 但是,只要思路掌握了,無論什麼樣的版本都是一個套路.