學前須知:
作為一名巨弱的數學競賽生&高數愛好者,數論知識無疑是我在oi最擅長的領域(沒有之一)了。那么我來結合網上的現有資料,以及我的個人見解,書寫一篇關於快速傅里葉變換的博客吧。
關於FFT我大約半年前掌握了,現有些許生疏,而且最近學了數學中有關拓撲學的DFT,有了些新的見解,所以寫了這篇。
此博客的作用是為了讓不會的同學快速入門學習,以及在我本人寫的過程中提升自我,無任何商業目的。PS:本人原創意識薄弱可能會從網上找一些現有資料。
此博客重在為初學者提供技巧以及結論的運用,至於理由我會粗略介紹,我認為還是要先知其然,學到手后再考慮證明什么的吧
如果你實在看了還不會的話背板子就行了,不用太過糾結。還有關於關於第二板塊的數學知識大家沒學過的話跳過就行就是科普一下。
另外此篇是給有數學基礎高中生已經自學了高中數學的同學看的。如果對虛數三角函數等知識不熟悉可以自行學習,此論文中不再贅述了。
1.FFT的作用以及功能
相信大家都曾聽聞過FFT,它是一種可以優化高精度乘法的高逼格算法。
首先我們要知道FFT是由DFT以更高效,快速計算的方式得到的。在介紹FFT之前我會再講一下關於DTF(離散傅里葉變換)的數學意義,這是其他的博客都不介紹的內容。為什么介紹呢,這是因為可以讓你了解數學中離散傅里葉變換是如何把信號從時間域轉到頻率域的,正如我們的現在的計算機的DFT也是用的這種轉換。即把系數轉成點值(過會會介紹具體)。如果數學沒有到大學數學本科or研究生水平的這段可以跳過了,感興趣的可以了解下。
2.數學上的DFT(可以忽略)
關於DFT:
如果你看到了這里,首先我默認你會狄拉克函數的基礎內容
1.狄拉克函數
此函數數學表達式為:
其時域波是周期為Ts的單位脈沖序列,也稱之為理想抽樣函數。
由於滿足了周期性,我們可以用傅里葉級數求其頻譜,傅里葉級數系數為:
其傅里葉展開為:
可見其頻率分布值在的離散頻率點,頻率分量幅值均為1/Ts。
引進沖擊函數的概念后,一般的周期函數也可求傅里葉變換,如下(cn為周期函數的傅里葉級數)
用狄拉克梳狀函數表示離散信號和離散頻譜,將給信號分析和處理帶來極大的方便。
2.離散時間傅里葉(DTFT)
既然你都看到這里了 ,我又默認你會傅里葉變換(調和分析范疇,非戰斗人員請撤離)
既然我們會了處理連續系統的,那離散系統怎么處理呢?
進行采樣,頻率為fs,沖擊采樣序列為:
取樣后:
連續信號的傅里葉變換定義:
采樣后:
同理於:
由於狄拉克函數篩選性質可知:
連續信號的傅里葉變換如下:
3.離散傅里葉變換(DFT)
我們容易知道,在當前的世界,信號的頻譜是不連續的,所以在DTFT中我們的頻率也需要離散化處理,所需要的就是DFT了。即具有周期特性離散信號的傅里葉級數(就是將無限長的離散限號進行截短至N個采樣點,然后將這個N個采樣點進行周期延拓 )
DFT:
且
IDFT(逆變換):
且
這樣DFT就介紹完了,至於證明還是算了(畢竟這里不是重點)。其實FFT本質也是DFT,只不過是利用DFT內部蘊藏的規律的一種簡化算法罷了。
3.關於多項式
首先想想我們要解決的問題是什么----兩個多項式想乘
容易知道:這兩個多項式我們可以有不同的表示方法
1.系數表示法
十分好理解的最朴素的表達方法
如果我們用這種表示法暴力算的話,復雜度肯定是O(n^2)的,顯然TLE。
2.點值表示法
我們不妨設兩多項式分別為f(x)與g(x)
具體表示比如,我選一個x0分別代入f(x)與g(x),那么f(x0)的值也就是原多項式a1的第0項的值,g(x0)的值也就是原多項式a2的第零項的值。如果不理解的話我列一下圖表吧。
f[x]={(x0,f[x0]),(x1,f[x1]).....(xn,f(xn)]
g[x]={x0,g[x0]),(x1,g[x1])......(xn,g[xn])
也就是說,當我們去一個x時,經過兩種不同的函數關系,可以得到不同的y值。那么我們只需要y值對於相乘就好了。復雜度O(n)
4.關於oi中的DFT
對於一般的多項式來說呢,我們當然可以隨便去n個x值代入計算。
也就是暴力計算出每個x的值,然后再代入y中。顯然這樣轉換復雜度是o(n^2)
怎么辦呢?
考慮DFT了。
我們看到我寫的第二部分的數學中的DFT可以發現,其本質原理就是就是將無限長的離散限號進行截短至N個采樣點,然后將這個N個采樣點進行周期延拓 。如果不知道的話沒關系,其實就是去一些“神奇”的x代入。
他規定了點值表示中的n個x為n個模長為一的復數,如圖:
然后要把這個單位圓N等分,如圖當N=8時
這些是當n=8時n等分后要取的點。記第k個點代表的復數值是,可以知道
。
那么此時稱為n次單位根。而且每個單位根都可以解出來:
那么這些Wn0,Wn1.......Wn n-1就是我們要代入的x0,x1 .......xn-1。(這里的W就是歐米伽蛤)
5.DFT優化后的FFT-快速傅里葉變換
你會發現其實你這樣暴力計算這些單位根后復雜度還是O(n^2).....
通過單位根的性質 ,考慮分治優化。
考慮多項式A(x)=a0+a1x+a2x^2+a3x^3+.....+an-1^(n-1)
我們可以根據下標的奇數性分成兩組,注意是“下標”的。
可以發現兩邊非常相似。
令A0(x)=a0 x^0+a2x+a4x^w+......+a(n-2)x^(n/2-1)
令A1(x)=a1 x^0+a3x+a5x^2+.....+a(n-1)x^(n/2-1)
顯然可以得到A(x)=A0(x^2)+x*A1(x^2)
蝴蝶變換:當我們關於A1(x)與A2(x)在的點值求出后,可以O(n)求出A(x)在
的點值了。
我們在處理過程中不斷遞歸分治,知道n=1時return。可以知道復雜度是o(nlog2n)的。
於是我在網上找了一個朴素的FFT板子
#include<complex> #define cp complex<double> void fft(cp *a,int n,int inv) { if (n==1)return; int mid=n/2; static cp b[MAXN]; fo(i,0,mid-1)b[i]=a[i*2],b[i+mid]=a[i*2+1]; fo(i,0,n-1)a[i]=b[i]; fft(a,mid,inv),fft(a+mid,mid,inv); fo(i,0,mid-1) { cp x(cos(2*pi*i/n),inv*sin(2*pi*i/n)); b[i]=a[i]+x*a[i+mid],b[i+mid]=a[i]-x*a[i+mid]; } fo(i,0,n-1)a[i]=b[i]; }
嗯,TLE了。又雙叒叕優化。
6.迭代版的FFT
可以發現,我們在進行FFT時會把各個系數分別放在兩側。一個系數的原來位置后最終位置的規律:
考慮二進制。
容易發現每個位置分治后其實就是二進制翻轉過來的。所以我們把每個數放在最后的位置,然后不斷向上還原同時求出點值就可以了。
復雜度 根號O(nlog2n)
FFT板子:
void fft(cp *a,int n,int inv) { int bit=0; while ((1<<bit)<n)bit++; fo(i,0,n-1) { rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1)); if (i<rev[i])swap(a[i],a[rev[i]]);//不加這條if會交換兩次(就是沒交換) } for (int mid=1;mid<n;mid*=2)//mid是准備合並序列的長度的二分之一 { cp temp(cos(pi/mid),inv*sin(pi/mid));//單位根,pi的系數2已經約掉了 for (int i=0;i<n;i+=mid*2)//mid*2是准備合並序列的長度,i是合並到了哪一位 { cp omega(1,0); for (int j=0;j<mid;j++,omega*=temp)//只掃左半部分,得到右半部分的答案 { cp x=a[i+j],y=omega*a[i+j+mid]; a[i+j]=x+y,a[i+j+mid]=x-y;//蝴蝶變換 } } } }
補充一點:關於IFFT(快速傅里葉逆變換)
當我們把兩個多項式從系數表示法變成點值表示法相乘后,還要將結果從點值表示法轉化為系數表示法的操作就是IFFT。
結論是:
一個多項式在分治的過程中乘上單位根的共軛復數,分治完的每一項除以n 即為原多項式的每一項系數
FFT和IFFT是很類似的。
7.有關於這個論文的后言:
這個論文我改過INF次了。一方面是oi里的DFT中有些錯誤以及表達問題,另一方面對於數學中的DFT做了許許多多的刪除。
其實我還是一開始打算講講傅里葉變換和傅里葉級數的,,但是后面都刪了還是。畢竟此論文主要闡述的還是oi里的算法
書呢是菲赫金哥爾茨的微積分學教程關於傅里葉級數理論的篇章(強推菲磚)
希望大家可以理解快速傅里葉變換的基本思想與流程叭,代碼的話多碼幾次背過就好了。
喜歡的記得給個贊哦qwq