傅里葉變換

本文轉載自查看原文 2021-09-23 10:18 733

1. 連續傅立葉變換（Continuous Fourier Transform）

對於時域連續函數 $x(t)$ ，它的傅立葉正變換（FT）定義為

$\mathcal F[ x(t)] = X(\omega) = \int_{-\infty}^\infty x(t) e^{-i\omega t}dt$ （用角頻率 $\omega$ 表示) 或者 $\mathcal F[ x(t)] = X(f) = \int_{-\infty}^\infty x(t) e^{-i2\pi f t}dt$ （用頻率 $f$ 表示, $\omega = 2\pi f$ )

傅立葉逆變換（inverse FT）定義為

$\mathcal F ^{-1}[ X(\omega)] = x(t) = \frac{1}{2\pi} \int_{-\infty}^\infty X(\omega) e^{+i\omega t}d\omega = \frac{1}{2\pi} \int_{-\infty}^\infty X(f) e^{+i2\pi f t}df$

2. 離散傅立葉變換（Discrete Fourier Transform, DFT）

如果有等時間間距 $\Delta t$ 的時域信號N個： $x[n],\text{where } n=0,\ldots,N-1$ ,它的離散傅立葉變換（DFT）定義為

$\mathcal F[x[n]]=X(e^{i\omega}) = X[k] = \sum_{n=0}^{N-1}x[n] e^{-i\frac{2\pi}{N}nk}, \;\;(k=0,\ldots,N-1)\\$

這里 $\omega = \frac{2\pi}{N} k$ . 相應的逆變換定義為

$\mathcal F^{-1}[X[k]]=X(e^{i\omega}) = x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k] e^{i\frac{2\pi}{N}nk}, \;\;(n=0,\ldots,N-1)\\$

3. 快速傅立葉變換（Fast Fourier Transform, FFT）

快速傅立葉變換FFT是一種更加高效的計算離散傅立葉變換的算法，有着 $O(N log_2 N)$ 的計算復雜度，比原始的DFT $O(N^2)$ 計算復雜度有更好的可擴展性。Cooley-Tuckey算法是一種常見的FFT實現方式，但是需要限制輸入數據點數 $N=2^m$ （2的乘方）形式。

如果你的數據數量不是2的幾次冪，可以通過在后面加0（zero padding）的方式湊成下一個2的乘方個數據，添加的0不會影響頻域譜，只會增加表觀數據密度，即多出來的頻域點其實是有效數據的線性內插值 Linear interpolation，如下圖所示：

對於zero padding證明如下：

如果記原始時域信號的離散傅立葉變換為 $F(e^{i\omega}) = \sum_{n=0}^{M-1} f[n]e^{-1\omega n}$ ，其中時域數據數為 $M$ . 我們添加0到后面，記zero padding的信號為 $g[n]=f[n]+zeros, (M\le n \le N-1)$ , N 是2的乘方。所以 $g[n], (n=0,\ldots,N-1)$ 的離散傅立葉變換為 $G(e^{i\omega}) = \sum_{n=0}^{N-1} g[n]e^{-i\omega n}= \sum_{n=0}^{M-1} f[n]e^{-i\omega n}+ \sum_{n=M}^{N-1} 0\times e^{-i\omega n}=F(e^{i\omega}) \qquad \square\\$

離散傅立葉變換得到的頻域信號是以 $\Delta f=\frac{1}{N\Delta t}$ 為間距的，最大頻率稱做Nyquist頻率，即為 $f_{Nyquist}=f_{max} = \frac{1}{2\Delta t}$ ，它只跟時域信號的步長有關。經過FFT變換的頻域信號對應的下列頻率

$\begin{align} f=\frac{\omega}{2\pi} = & 0, \Delta f, 2\Delta f, \ldots, (N/2-1)\Delta f,\\ & \pm (N/2)\Delta f, -(N/2-1)\Delta f, -(N/2-2)\Delta f, \ldots, -\Delta f. \end{align}\\$

也就是說，前一半是正的頻率，后一半是負的頻率，中間點是正負最大頻率。一般畫圖的話，需要把后一半挪到前面。

連續傅里葉變換

一般情況下，若“傅里葉變換”一詞的前面未加任何限定語，則指的是“ 連續傅里葉變換”。“連續傅里葉變換”將平方可積的函數

表示成復指數函數的積分形式：

上式其實表示的是連續傅里葉變換的逆變換，即將時間域的函數表示為頻率域的函數

的積分。反過來，其正變換恰好是將頻率域的函數

表示為時間域的函數

的積分形式。一般可稱函數

為原函數，而稱函數

為傅里葉變換的像函數，原函數和像函數構成一個傅里葉變換對（transform pair）。

當

為奇函數（或偶函數）時，其余弦（或正弦）分量為零，而可以稱這時的變換為余弦變換（或正弦變換）。

離散傅里葉變換

為了在科學計算和數字信號處理等領域使用計算機進行傅里葉變換，必須將函數定義在離散點上而非連續域內，且須滿足有限性或周期性條件。這種情況下，序列

的離散傅里葉變換（discrete Fourier transform, DFT）為。

其逆變換為

直接使用DFT的定義計算的計算復雜度為

，而快速傅里葉變換（fast Fourier transform, FFT）可以將復雜度改進為

。計算復雜度的降低以及數字電路計算能力的發展使得DFT成為在信號處理領域十分實用且重要的方法。

如果寫成一般形式，用 $x[n]$ 替換原信號，用 $N$ 替換40，得到

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})$

問題來了，雖然貌似聯系很緊密，但這怎么跟DFT的公式長得不一樣。。。DFT的公式應該是這樣的：

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

就算用歐拉公式展開，我們得到的是

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})-j\sum_{n=0}^{N-1}x[n]sin(2\pi\frac{kn}{N})$

作者：肖暢
鏈接：https://www.zhihu.com/question/21314374/answer/542909849
來源：知乎
著作權歸作者所有。商業轉載請聯系作者獲得授權，非商業轉載請注明出處。

覺得現有答案都是從連續傅里葉變換引出離散的，各位大神功底很扎實啊。但問題是很多同學就是因為對那些繁瑣的 $\int_{-\infty}^{\infty}$ 和 $e^{jwt}$ 感到暈厥，才轉到沒有積分，只有 $\sum$ ，看着比較簡單的DFT來試試。所以DFT這樣明顯更簡單的一個東西為什么一定要從連續傅里葉變換出發來講呢？換個思路不行嗎？

本文直接從DFT出發，不需要連續傅里葉變換的知識。

DFT做了什么事？

對一段有限長的離散信號，找出它含有的各個頻率的正弦波分量。（這是句大家都知道的廢話）

但如果換個單獨針對DFT的解釋就會好理解很多：

在長度為N的離散信號中，針對k=(0,1,2...)，分別找出在長度N內振動k個周期的三角波分量的權值

怎么做呢？舉個例子，大家可以打開MATLAB一起做

拋開時間拋開采樣頻率，只看點的個數，我們對某個余弦信號在兩個周期內采樣了40次：

$x[n]=cos(2\pi \frac{n}{20}), n=0,1,\cdots,39$

n = 0:39; y = cos(2*pi*(2*n)/40); stem(n,y);

現在我們想通過DFT知道它在40次采樣時間內振動了幾個周期，算法是怎么處理的呢？其實很暴力：

事先選取40個長度為40個點的基信號，它們分別長這樣：

第一個，40次采樣內振動0個周期： $cos(2\pi\frac{0n}{40})$ ，即常值

40次采樣內振動0個周期

第二個，40次采樣內振動1個周期： $cos(2\pi \frac{n}{40})$

40次采樣內振動1個周期

第三個，40次采樣內振動2個周期： $cos(2\pi \frac{2n}{40})$

40次采樣內振動2個周期

第四個，40次采樣內振動3個周期： $cos(2\pi \frac{3n}{40})$

40次采樣內振動3個周期

第五個，40次采樣內振動4個周期： $cos(2\pi\frac{4n}{40})$

40次采樣內振動4個周期

……

以此類推一直到40個采樣內振動39個周期。

接下來，對於上述40個基信號，判斷它們跟原信號的相關程度。代數里判斷兩個長度相同的離散函數x[k],y[k]的相關程度，是用把它們在同一點的函數值相乘，再對得到的積求和的方法（就是向量求內積的方法），即下面的公式：

$correlation(x,y) = \sum_k{x[k]y[k]}$

這個值越大，x[k]和y[k]長得越像。於是DFT把

$cos(2\pi\frac{0n}{40})$ 到 $cos(2\pi\frac{39n}{40})$ 這40個基函數都跟當前函數 $cos(2\pi\frac{2n}{40})$ 比較了一下，發現

$cos(2\pi \frac{2n}{40})$ 和 $cos(2\pi \frac{38n}{40})$ 跟它長得最像！

這很顯然，因為 $cos(2\pi\frac{38n}{40})=cos(2\pi n-2\pi\frac{38n}{40})=cos(2\pi\frac{2n}{40})$

於是 $X_2=X_{38}=\sum_{n=0}^{39}cos(2\pi\frac{2n}{40})cos(2\pi\frac{2n}{40})=20$

下面，如果我們把這40次每次比較的correlation值記下來，就得到了原信號在每個頻率上的分量大小，比如，我要看它跟40次采樣中振動8次的那個基信號相似程度，就是：

把這40個correlation值記下來，就得到了原信號的頻域X：

stem(n, abs(fft(y)));

可以看到 $X_2 = X_{38} = 20$ ，其他值都是0

如果寫成一般形式，用 $x[n]$ 替換原信號，用 $N$ 替換40，得到

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})$

問題來了，雖然貌似聯系很緊密，但這怎么跟DFT的公式長得不一樣。。。DFT的公式應該是這樣的：

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

就算用歐拉公式展開，我們得到的是

$X_k=\sum_{n=0}^{N-1}x[n]cos(2\pi\frac{kn}{N})-j\sum_{n=0}^{N-1}x[n]sin(2\pi\frac{kn}{N})$

這又是為什么？

這是因為，對於一個信號，如果只跟余弦函數比較，會損失一些信息，比如相位。

如果我們的原信號有一些相位偏移， $x=cos(2\pi \frac{2n}{40} + \frac{\pi}{4})$

如果對這個函數同樣按照上面的方法計算頻域，結果會有一些不一樣：

n = 0:39; x = cos(2 * pi * (2 * n)/40 + pi/4); stem(n, real(fft(x)));

雖然對應頻率還是篩選出來了，但是權值大小不一樣了，如果計算一下得知

$X_2=X_{38}=\sum_{n=0}^{39}cos(2\pi\frac{2n}{40}+\frac{\pi}{4})cos(2\pi\frac{2n}{40})=10\sqrt{2}$

在相位偏移前，這個值是20，是原來的 $\frac{\sqrt{2}}{2}$ 倍

注意如果我們再找一個信號 $y$ ，沒有相位偏移，但是直接把幅值砍到 $\frac{\sqrt{2}}{2}$ ，即：

$y=\frac{\sqrt{2}}{2}cos(2\pi \frac{2n}{40})$

那么這個信號的

$Y_2=Y_{38}=\sum_{n=0}^{39}\frac{\sqrt{2}}{2}\cos(2\pi\frac{2n}{40})\cos(2\pi\frac{2n}{40})=10\sqrt{2}$

跟 $x$ 信號的結果一模一樣，這樣我們就無法通過頻域恢復信號了，或者說我們損失了信息

解決方法是另選一組以正弦函數（實際上選了負正弦）為基准的“基信號”，即 $-\sin(2\pi \frac{0n}{40})$ 到 $-\sin(2\pi \frac{39n}{40})$ ，計算另一組原信號與正弦基的相關系數，這兩組系數一起作為DFT的最終結果。而復數只是一個工具，用來方便地同時存儲兩組計算結果。當然它還有一個好處就是能夠比較直觀地表現出模和相位。

選負正弦還是正弦做基信號其實無所謂，只是最后的結果算出來的相位反一下而已，幅值是一樣的。如果一個信號跟某頻率余弦和負正弦的相關系數分別為a,b，那么代表這個信號差不多型如

$a\cos(\frac{kn}{N})-b\sin(\frac{kn}{N})$

根據高中數學可以求得其模為 $\sqrt{a^2+b^2}$ ，（相對余弦的）相位為 $\arctan(\frac{b}{a})$ ，這與復數 $a+bi$ 的模和相位是相同的，因此DFT的公式

$X_k=\sum_{n=0}^{N-1}x[n]e^{-2\pi j\frac{kn}{N}}$

相當於同時把x[n]做了跟N個余弦基、N個負正弦基的比對，結果用N個復數存儲。如果想要看某個頻率的相位和模，就看對應復數的相位和模。

我們再來看看上面有相位偏移的那個例子：

原信號： $x[n]=cos(2\pi\frac{2n}{40}+\frac{\pi}{4})$

與余弦比對： $X_2 =\sum_{n=0}^{39} cos(2\pi\frac{2n}{40}+\frac{\pi}{4})cos(2\pi\frac{2n}{40})=10\sqrt{2}$

與負正弦比對： $X_2 =\sum_{n=0}^{39} cos(2\pi\frac{2n}{40}+\frac{\pi}{4})sin(-2\pi\frac{2n}{40})=10\sqrt{2}$

在40個點內振動兩個周期這個頻率上，其DFT的結果為 $10\sqrt{2}+10\sqrt{2}j$

其模為20，與相位偏移前相同。

其相位為 $\frac{\pi}{4}$ ，也沒有問題。

如果將原信號變為 $x[n]=sin(2\pi\frac{2n}{40}+\frac{\pi}{4})$ ，會求得該頻率DFT結果為 $10\sqrt{2}-10\sqrt{2}j$ ，求得其相位為 $-\frac{\pi}{4}$ 。因此，根據DFT結果求得的相位是相對余弦信號的相位。

對於實際信號

上面做的分析沒有考慮時間和采樣頻率，但只要給定一個采樣頻率，時域的橫軸就可以換算成時間，頻域的橫軸就可以換算成頻率。

DFT所做的只是選取了在給定長度內振動了整數個周期的正弦和余弦波作為基。

上面的信號是40個采樣，如果給一個采樣頻率是100Hz，那么信號長度就是0.4s，原信號在40個采樣內振動了兩個周期，可以算出其頻率為5Hz

對於頻域，每個“在40個采樣內振動了k個周期”的基信號的實際頻率為 $\frac{100}{40}k=2.5k$ ，也就是說頻域圖中的每一個點代表2.5Hz，這個系統的頻域分辨率是2.5Hz，所以我們從 $X_2=20$ 得知原信號的實際頻率為 $2\times2.5=5Hz$

FFT(Fast Fourier Transform)是什么

FFT是 Cooley & Tuket 兩人1965年提出的快速計算DFT的算法。這背后還有個故事，美國和蘇聯1963年簽了個核試驗禁令，互相約定大家都不搞核試驗了。但是美國不放心啊，怕毛子說一套做一套，肯尼迪就請了一堆科學家開會，說想搞一套不用去蘇聯檢查就能探測到核試驗的設備。美國在蘇聯周圍放了一圈聲波探測儀，但是遇到個問題——DFT算得太慢了。正巧 Tuket 和IBM的一個叫 Richard Garwin 的出席了那次會議。會上Garwin就想起來Tuket和Cooley好像合作過一個類似的算法，就讓Tuket去找Cooley，但是隱藏了真實目的，告訴Cooley這個算法是為了探測氦3晶體的自旋周期……（也是醉了orz）。兩人的算法最后於1965年發表，極大提升了DFT的計算速度，甚至被后世譽為“20世紀最偉大的工程算法”。

從DFT的公式可以看出其算法復雜度是 $O(n^2)$ 的，兩人的FFT算法表明，在N為合數的情況下 $N=n_1n_2$ ，原長度為N的DFT可以分解為兩個長度分別為 $n_1$ 和 $n_2$ 的DFT！

當然算法不停止於此，如果 $n_1$ 和 $n_2$ 還可以繼續分解，就可以分解為更多小DFT來計算。

使用了算法中的分治（divide and conquer）策略。

而由於DFT的長度其實是可以任意選擇的，因此通常在使用FFT時，都選擇2的整數次冪作為FFT長度，這樣可以一直分解到N/2個長度為2的DFT，復雜度直接降到 $O(nlogn)$

更詳細的解答，請移步我做的視頻：https://www.bilibili.com/video/av44600709

想從理論層面更深入了解DFT的，可以看我另外一個視頻：https://www.bilibili.com/video/BV1a

0x00 寫在前面

為了讓更多人能夠看到這個教程，希望大家收藏之前，也要點贊哦！！！蟹蟹大家的認可和鼓勵。

傅里葉變換

快速傅里葉變換（Fast Fourier Transform，FFT）是一種可在 $O(nlogn)$ 時間內完成的離散傅里葉變換（Discrete Fourier transform，DFT）算法。

在算法競賽中的運用主要是用來加速多項式的乘法。

考慮到兩個多項式 $A(x),B(x)$ 的乘積 $C(x)$ ，假設 $A(x)$ 的項數為 $n$ ，其系數構成的 $n$ 維向量為 $(a_0,a_1,a_2,...,a_{n-1})$ ， $B(x)$ 的項數為 $m$ ，其系數構成的 $m$ 維向量為 $(b_0,b_1,b_2,...,b_{m-1})$ 。

我們要求 $C(x)$ 的系數構成的 $n+m-1$ 維的向量，先考慮朴素做法。

可以用這段代碼表示：

for ( int i = 0 ; i < n ; ++ i ) for ( int j = 0 ; j < m ; ++ j ) { c [i + j] += a [i] * b [j] ; }

思路非常清晰，其時間復雜度是 $O(n^2)$ 的。

所以我們來學習快速傅里葉變換。

0x01 關於多項式

多項式有兩種表示方法，系數表達法與點值表達法

多項式的系數表示法

設多項式 $A(x)$ 為一個 $n-1$ 次的多項式，顯然，所有項的系數組成的系數向量 $(a_0,a_1,a_2,...,a_{n-1})$ 唯一確定了這個多項式。

$A(x)=\sum_{i=0}^{n-1}a_i\cdot x^i$

多項式的點值表示法

將一組互不相同的 $(x_0,x_1,x_2,...,x_n)$ （叫插值節點）分別帶入 $A(x)$ ，得到 $n$ 個取值 $(y_0,y_1,y_2,...,y_n)$ .

其中

$y_i=\sum_{j=0}^{n-1}a_j\cdot x_i^j$

定理:

一個 $n-1$ 次多項式在 $n$ 個不同點的取值唯一確定了該多項式。

證明：

假設命題不成立，存在兩個不同的 $n-1$ 次多項式 $A(x),B(x)$ ，滿足對於任何 $i \in [0,\ n - 1]$ ，有 $A(x_i) = B(x_i)$ 。
令 $C(x) = A(x) - B(x)$ ，則 $C(x)$ 也是一個 $n-1$ 次多項式。對於任何 $i \in [0,\ n - 1]$ ，都有 $C(x_i) = 0$ 。
即 $C(x)$ 有 $n$ 個根，這與代數基本定理（一個 $n-1$ 次多項式在復數域上有且僅有 $n-1$ 個根）相矛盾，故 $C(x)$ 並不是一個 $n-1$ 次多項式，推到矛盾。
原命題成立，證畢。

如果我們按照定義求一個多項式的點值表示，時間復雜度為 $O(n^2)$

已知多項式的點值表示，求其系數表示，可以使用插值。朴素的插值算法時間復雜度為 $O(n^2)$ 。

關於多項式的乘法

已知在一組插值節點 $(x_0,x_1,x_2,...,x_n)$ 中 $A(x),B(x)$ (假設個多項式的項數相同，沒有的視為 $0$ )的點值向量分別為 $(y_{a0},y_{a1},y_{a_2},...,y_{an}),(y_{b0},y_{b1},y_{b_2},...,y_{bn})$ ，那么 $C(x)=A(x)\cdot B(x)$ 的點值表達式可以在 $O(n)$ 的時間內求出，為 $(y_{a0}\cdot y_{b0},y_{a1}\cdot y_{b1},y_{a_2}\cdot y_{b2},...,y_{an}\cdot y_{bn})$ 。

因為 $C(x)$ 的項數為 $A(x),B(x)$ 的項數之和。

設 $A(x),B(x)$ 分別有 $n,m$ 項所以我們帶入的插值節點有至少有 $n+m$ 個。

如果我們能快速通過點值表式求出系數表示，那么就搭起了它們之間的一座橋了。

這也是快速傅里葉變換的基本思路，由系數表達式到點值表達式到結果的點值表達式再到結果的系數表達式。

0x02 關於復數的基本了解

我們把形如 $a+bi$ 這樣的數叫做復數，復數集合用 $C$ 來表示。其中 $a$ 稱為實部 $(real\;part)$ ， $b$ 稱為虛部 $(imaginary\;part)$ , $i$ 為虛數單位，指滿足 $x^2=-1$ 的一個解 $\sqrt{-1}$ ；此外，對於這樣對復數開偶次冪的數叫做虛數 $( imaginary\;number)$ .

每一個復數 $a+bi$ 都對應了一個平面上的向量 $(a,b)$ 我們把這樣的平面稱為復平面 $(complex\;plane)$ ，它是由水平的實軸與垂直的虛軸建立起來的復數的幾何表示。

故每一個復數唯一對應了一個復平面上的向量，每一個復平面上的向量也唯一對應了一個復數。其中 $0$ 既被認為是實數，也被認為是虛數。

其中復數 $z=a+bi$ 的模長 $\left| z \right|$ 定義為 $z$ 在復平面的距離到原點的距離， $\left| z \right|=\sqrt{a^2+b^2}$ 。幅角 $\theta$ 為實軸的正半軸正方向（逆時針）旋轉到 $z$ 的有向角度。

由於虛數無法比較大小。復數之間的大小關系只存在等於與不等於兩種關系，兩個復數相等當且僅當實部虛部對應相等。對於虛部為 $0$ 的復數之間是可以比較大小的，相當於實數之間的比較。

復數之間的運算滿足結合律，交換律和分配律。

由此定義復數之間的運算法則：

$(a+bi)+(c+di)=(a+c)+(b+d)i$

$(a+bi)-(c+di)=(a-c)+(b-d)i$ $(a+bi)\cdot(c+di)=ac+adi+bci+bdi^2=(ac-bd)+(ad+cb)i$

$\frac{a+bi}{c+di}=\frac{(a+bi)\cdot(c-di)}{(c+di)\cdot(c-di)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2}=\frac{(ac+bd)}{c^2+d^2}+\frac{(bc-ad)i}{c^2+d^2}$

復數運算的加法滿足平行四邊形法則，乘法滿足幅角相加，模長相乘。

對於一個復數 $z=a+bi$ ,它的共軛復數是 $z^{'}=a-bi$ , $z^{'}$ 稱為 $z$ 的復共軛 $(complex\;conjugate)$ .

共軛復數有一些性質

$z\cdot z^{'}=a^2+b^2$

$\left| z \right| = \left | z^{'} \right|$

0x03 復數中的單位根

復平面中的單位圓

其中 $\bar{OA}$ 單位根，表示為 $e^{i\theta}$ ，可知 $e^{i\theta}=cos\theta+i\cdot sin\theta$

（順便一提著名的歐拉幅角公式 $e^{i\pi}+1=0$ 其實是由定義來的...）

將單位圓等分成 $n$ 個部分（以單位圓與實軸正半軸的交點一個等分點），以原點為起點，圓的這 $n$ 個 $n$ 等分點為終點，作出 $n$ 個向量。

其中幅角為正且最小的向量稱為 $n$ 次單位向量，記為 $\omega_{n}^{1}$ 。

（有沒有大佬幫我補張圖啊，畫不來）

其余的 $n-1$ 個向量分別為 $\omega_{n}^{2},\omega_{n}^{3},......,\omega_{n}^{n}$ ，它們可以由復數之間的乘法得來 $w_{n}^{k}=w_{n}^{k-1}\cdot w_{n}^{1}\ (2 \leq k \leq n)$ 。

容易看出 $w_{n}^{n}=w_{n}^{0}=1$ 。

對於 $w_{n}^{k}$ ,它事實上就是 $e^{2\pi\cdot \frac{k}{n}i}$ 。

所以 $\omega_{n}^{k}=cos(2\pi\cdot\frac{k}{n})+i\cdot sin(2\pi\cdot\frac{k}{n})$

關於單位根有兩個性質

性質一（又稱為折半引理）：

$\omega_{2n}^{2k}=\omega_{n}^{k}$

證明一：

由幾何意義，這兩者表示的向量終點是一樣的。

證明二：

由計算的公式：

$\omega_{2n}^{2k}=cos(2\pi\frac{2k}{2n})+i\cdot sin(2\pi\frac{2k}{2n})=cos(2\pi\frac{k}{n})+i\cdot sin(2\pi\frac{k}{n})=\omega_{n}^{k}$

其實由此我們可以引申出

$\omega_{mn}^{mk}=\omega_{n}^{k}$

性質二（又稱為消去引理）

$\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$

證明一：

由幾何意義，這兩者表示的向量終點是相反的，左邊較右邊在單位圓上多轉了半圈。

證明二：

由計算的公式：

$\omega_{n}^{k+\frac{n}{2}}=cos(2\pi\frac{k+\frac{n}{2}}{n})+i\cdot sin(2\pi\frac{k+\frac{n}{2}}{n})=cos(2\pi\frac{k}{n}+\pi)+i\cdot sin(2\pi\frac{k}{n}+\pi)=-cos(2\pi\frac{k}{n})-i\cdot sin(2\pi\frac{k}{n})=-\omega_{n}^{k}$

最后一步由三角恆等變換得到。

0x04 離散傅里葉變換（Discrete Fourier Transform）

首先我們單獨考慮一個 $n$ 項（ $n=2^x$ ）的多項式 $A(x)$ ，其系數向量為 $(a_0,a_1,a_2,...,a_{n-1})$ 。我們將 $n$ 次單位根的 $0$ ~ $n-1$ 次冪分別帶入 $A(x)$ 得到其點值向量 $(A(w_n^{0}),A(w_n^{1}),A(w_n^{2}),...,A(w_n^{n-1}))$ 。

這個過程稱為離散傅里葉變換（Discrete Fourier Transform）。

如果朴素帶入，時間復雜度也是 $O(n^2)$ 的。

所以我們必須要利用到單位根 $\omega$ 的特殊性質。

對於 $A(x)=a_0+a_1\cdot x^1+a_2\cdot x^2+a_3\cdot x^3+...+a_{n-1}\cdot x^{n-1}$

考慮將其按照奇偶分組

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+(a_1\cdot x^1+a_3\cdot x^3+a_{5}\cdot x^{5}+...+a_{n-1}\cdot x^{n-1})$

$A(x)=(a_0+a_2\cdot x^2+a_{4}\cdot x^{4}...+a_{n-2}\cdot x^{n-2})+x\cdot(a_1+a_3\cdot x^2+a_{5}\cdot x^{4}+...+a_{n-1}\cdot x^{n-2})$

令

$A1(x)=(a_0+a_2\cdot x+a_{4}\cdot x^{2}...+a_{n-2}\cdot x^{\frac{n-2}{2}})$

$A2(x)=(a_1+a_3\cdot x+a_{5}\cdot x^{2}...+a_{n-1}\cdot x^{\frac{n-2}{2}})$

則可得到

$A(x)=A1(x^2)+x\cdot A2(x^2)$

分類討論

設 $0\leq k\leq \frac{n}{2}-1$ ， $k\in Z$

$A(\omega_{n}^{k})=A1(\omega_{n}^{2k})+\omega_{n}^{k}\cdot A2(\omega_{n}^{2k})$

由上文提到的折半引理

$A(\omega_{n}^{k})=A1(\omega_{\frac{n}{2}}^{k})+\omega_{n}^{k}\cdot A2(\omega_{\frac{n}{2}}^{k})$

對於 $\frac{n}{2}\leq k+\frac{n}{2}\leq n-1$

$A(\omega_{n}^{k+\frac{n}{2}})=A1(\omega_{n}^{2k+n})+\omega_{n}^{k+\frac{n}{2}}\cdot A2(\omega_{n}^{2k+n})$

其中 $\omega_{n}^{2k+n}=\omega_{n}^{2k}\cdot \omega_{n}^{n}=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}$

由消去引理 $\omega_{n}^{k+\frac{n}{2}}=-\omega_{n}^{k}$

故 $A(\omega_{n}^{k+\frac{n}{2}})=A1(\omega_{\frac{n}{2}}^{k})-\omega_{n}^{k}\cdot A2(\omega_{\frac{n}{2}}^{k})$

注意， $k$ 與 $k+\frac{n}{2}$ 取遍了 $[0,n-1]$ 中的 $n$ 個整數，保證了可以由這 $n$ 個點值反推解出系數（上文已證明）。

於是我們可以知道

如果已知了 $A1(x),A2(x)$ 分別在 $\omega_{\frac{n}{2}}^{0},\omega_{\frac{n}{2}}^{1},...,\omega_{\frac{n}{2}}^{\frac{n}{2}-1},$ 的取值，可以在 $O(n)$ 的時間內求出 $A(x)$ 的取值。

而 $A1(x),A2(x)$ 都是 $A(x)$ 一半的規模，顯然可以轉化為子問題遞歸求解。

時間復雜度：

$T(n)=2T(\frac{n}{2})+O(n)=O(nlogn)$

0x05 離散傅里葉反變換（Inverse Discrete Fourier Transform）

使用快速傅里葉變換將點值表示的多項式轉化為系數表示，這個過程叫做離散傅里葉反變換（Inverse Discrete Fourier Transform）。

即由 $n$ 維點值向量 $(A(x_0),A(x_1),...,A(x_{n-1}))$ 推出 $n$ 維系數向量 $(a_0,a_1,...,a_{n-1})$ 。

設 $(d_0,d_1,...,d_{n-1})$ 為 $(a_0,a_1,...,a_{n-1})$ 得到的離散傅里葉變換的結果。

我們構造一個多項式 $F(x)=d_0+d_1\cdot x+d_2\cdot x^2+...+d_{n-1}\cdot x^{n-1}$

設向量 $(c_0,c_1,...,c_{n-1})$ 中

$c_k$ 為 $F(x)$ 在 $x=\omega_{n}^{-k}$ 的點值表示

即 $c_k=\sum_{i=0}^{n-1}d_i\cdot(\omega_{n}^{-k})^i$ ，

我們考慮對 $d_i$ 進行還原

於是

$c_k=\sum_{i=0}^{n-1}[\sum_{j=0}^{n-1}a_j\cdot(\omega_{n}^{i})^j]\cdot(\omega_{n}^{-k})^i$

由和式的性質

$c_k=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{i})^j\cdot(\omega_{n}^{-k})^i$

$=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$

令 $S(j,k)=\sum_{i=0}^{n-1}(\omega_{n}^{i})^{j-k}$

對其進行化簡

設 $j-k=\delta$

則 $S(j,k)=\omega_{n}^{0}+\omega_{n}^{\delta}+\omega_{n}^{2\delta}+...+\omega_{n}^{(n-1)\delta}$

其公比為 $\omega_{n}^{\delta}$

當 $\omega_{n}^{\delta}=1$ 即 $\delta=0$ 時

$S(j,k)=n$ 此時 $\delta=0\Rightarrow j-k=0 \Rightarrow j=k$

當 $\omega_{n}^{\delta}\ne1$ 即 $\delta\ne0$ 時

由等比數列求和公式

$S(j,k)=\frac{\omega_{n}^{0}[(\omega_{n}^{\delta})^{n}-1]}{\omega_{n}^{\delta}-1}=\frac{1 [(\omega_{n}^{n})^{\delta}-1]}{\omega_{n}^{\delta}-1}=\frac{[(1)^{\delta}-1]}{\omega_{n}^{\delta}-1}=\frac{0}{\omega_{n}^{\delta}-1}=0$ ,此時 $j\ne k$ .

所以

$S(j,k)=[j=k]\cdot n$

將 $S(j,k)$ 帶入原式

$c_k=\sum_{j=0}^{n-1}a_j\cdot S(j,k)=\sum_{j=0}^{n-1}a_j\cdot [j=k]\cdot n=a_k\cdot n$

所以 $a_k=\frac{c_k}{n}$ .

其中 $a_k$ 為原多項式 $A(x)$ 的系數向量 $(a_0,a_1,...,a_n)$ 中的 $a_k$ .

由此得到：

對於多項式 $A(x)$ 由插值節點 $(\omega_{n}^{0},\omega_{n}^{1},\omega_{n}^{2},...,\omega_{n}^{n-1})$ 做離散傅里葉變換得到的點值向量 $(d_0,d_1,...,d_{n-1})$ 。我們將 $(\omega_{n}^{0},\omega_{n}^{-1},\omega_{n}^{-2},...,\omega_{n}^{-(n-1)})$ 作為插值節點， $(d_0,d_1,...,d_{n-1})$ 作為系數向量，做一次離散傅里葉變換得到的向量每一項都除以 $n$ 之后得到的 $(\frac{c_0}{n},\frac{c_1}{n},...,\frac{c_{n-1}}{n})$ 就是多項式的系數向量 $(a_0,a_1,...,a_{n-1})$ 。

注意到 $\omega_{n}^{-k}$ 是 $\omega_{n}^{k}$ 的共軛復數。

這個過程稱為離散傅里葉反變換。

0x06 關於FFT在C++的實現

首先要解決復數運算的問題，我們可以使用C++STL自帶的 $std :: complex < T >$ 依照精度要求 $T$ 一般為 $double,long\;double$ 。

也可以自己封裝，下面是我封裝的復數類。

struct Complex { double r, i ; Complex ( ) { } Complex ( double r, double i ) : r ( r ), i ( i ) { } inline void real ( const double& x ) { r = x ; } inline double real ( ) { return r ; } inline Complex operator + ( const Complex& rhs ) const { return Complex ( r + rhs.r, i + rhs.i ) ; } inline Complex operator - ( const Complex& rhs ) const { return Complex ( r - rhs.r, i - rhs.i ) ; } inline Complex operator * ( const Complex& rhs ) const { return Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator /= ( const double& x ) { r /= x, i /= x ; } inline void operator *= ( const Complex& rhs ) { *this = Complex ( r * rhs.r - i * rhs.i, r * rhs.i + i * rhs.r ) ; } inline void operator += ( const Complex& rhs ) { r += rhs.r, i += rhs.i ; } inline Complex conj ( ) { return Complex ( r, -i ) ; } } ;

我們由上面的分析可以得到這個遞歸的寫法。

bool inverse = false ; inline Complex omega ( const int& n, const int& k ) { if ( ! inverse ) return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ) ; return Complex ( cos ( 2 * PI / n * k ), sin ( 2 * PI / n * k ) ).conj ( ) ; } inline void fft ( Complex *a, const int& n ) { if ( n == 1 ) return ; static Complex buf [N] ; const int m = n >> 1 ; for ( int i = 0 ; i < m ; ++ i ) { buf [i] = a [i << 1] ; buf [i + m] = a [i << 1 | 1] ; } memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ; Complex *a1 = a, *a2 = a + m; fft ( a1, m ) ; fft ( a2, m ) ; for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega ( n, i ) ; buf [i] = a1 [i] + t * a2 [i] ; buf [i + m] = a1 [i] - t * a2 [i] ; } memcpy ( a, buf, sizeof ( Complex ) * ( n + 1 ) ) ; }

但是這樣的 $FFT$ 要用到輔助數組，並且常數比較大。

能不能優化呢？

我們把每一次分組的情況推演出來

遞歸分類的每一層

觀察到每一個位置的數其實都是原來位置上的數的二進制后 $log_{2}n$ 位 $reverse$ 了一下。

於是我們可以想，先將原數組調整成最底層的位置（很好調整吧）。

然后從倒數第二層由底向上計算。

這就是我們一般用來實現 $FFT$ 的 $Cooley-Tukey$ 算法。

考慮怎么合並？

在 $Cooley-Tukey$ 算法中，合並操作被稱作是蝴蝶操作。

慮合並兩個子問題的過程，這一層有 $n$ 項需要處理。假設 $A_1(\omega_{ \frac{n}{2} } ^ k)$ 和 $A_2(\omega_{ \frac{n}{2} } ^ k)$ 分別存在 $a(k)$ 和 $a(\frac{n}{2} + k)$ 中， $A(\omega_n ^ {k})$ 和 $A(\omega_n ^ {k + \frac{n}{2} })$ 將要被存放在 $buf(k)$ 和 $buf(\frac{n}{2} + k)$ 中，合並的單位操作可表示為

$buf(k):=a(k)+\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

$buf(k+\frac{n}{2}):=a(k)-\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

只要將合並順序換一下，再加入一個臨時變量，合並過程就可以在原數組中進行。

令 $t:=\omega_{n}^{k}\cdot a(k+\frac{n}{2})$

合並過程如下：

$a(k+\frac{n}{2}):=a(k)-t$

$a(k):=a(k)+t$ 。

至此，我們可以給出 $Cooley-Tukey$ 算法的實現。

struct FastFourierTransform { Complex omega [N], omegaInverse [N] ; void init ( const int& n ) { for ( int i = 0 ; i < n ; ++ i ) { omega [i] = Complex ( cos ( 2 * PI / n * i), sin ( 2 * PI / n * i ) ) ; omegaInverse [i] = omega [i].conj ( ) ; } } void transform ( Complex *a, const int& n, const Complex* omega ) { for ( int i = 0, j = 0 ; i < n ; ++ i ) { if ( i > j ) std :: swap ( a [i], a [j] ) ; for( int l = n >> 1 ; ( j ^= l ) < l ; l >>= 1 ) ; } for ( int l = 2 ; l <= n ; l <<= 1 ) { int m = l / 2; for ( Complex *p = a ; p != a + n ; p += l ) { for ( int i = 0 ; i < m ; ++ i ) { Complex t = omega [n / l * i] * p [m + i] ; p [m + i] = p [i] - t ; p [i] += t ; } } } } void dft ( Complex *a, const int& n ) { transform ( a, n, omega ) ; } void idft ( Complex *a, const int& n ) { transform ( a, n, omegaInverse ) ; for ( int i = 0 ; i < n ; ++ i ) a [i] /= n ; } } fft ;