Fourier transform
「傅立葉轉換」是雙射函數，輸入輸出都是一串複數，可以是離散數列或者連續函數，各有對應名稱。混淆視聽罷了。
輸入丨輸出丨名稱
一一十一一十一一一一一一一一一一一一一一一一一一一
離散丨離散丨discrete Fourier transform
離散丨連續丨discrete-time Fourier transform
連續丨離散丨Fourier series
連續丨連續丨Fourier transform
離散到離散的傅立葉轉換，輸入和輸出都是一串複數數列。
電腦做運算，數值皆離散。本文介紹離散版本。
連續到連續的傅立葉轉換，輸入和輸出都是一個ℝ ⇨ ℂ函數。
連續版本是離散版本的推廣：輸入輸出無限密無限長。
discrete Fourier transform物理意義
N個複弦波，頻率是0倍到N-1倍，分別是e𝑖⋅(2π/N)⋅0⋅t、e𝑖⋅(2π/N)⋅1⋅t、……、e𝑖⋅(2π/N)⋅(N-1)⋅t。寫成代數是e𝑖⋅(2π/N)⋅f⋅t。（複弦波很難畫，故圖例為實弦波。）
輸入數列與一個波，靠左對齊。N個對應位置，相除後求和，得到一個輸出數值。可以簡單想做：輸入數列除以波，求比例。
輸入數列，分別除以N個波，得到N個輸出數值，形成輸出數列。這就是傅立葉轉換。
正向傅立葉轉換：一個複雜的波，拆解成N個平穩的波，頻率是0倍到N-1倍，振幅與相位是N個輸出數值的強度與相位。
逆向傅立葉轉換：N個平穩的波，頻率是0倍到N-1倍，分別乘上振幅、添上相位，疊加成一個複雜的波。
discrete Fourier transform數學公式
正向傅立葉轉換
N-1
y[f] = ∑ { x[t] ÷ e𝑖⋅(2π/N)⋅f⋅t } ÷ √N
t=0
N-1
= ∑ { x[t] ⋅ e-𝑖⋅(2π/N)⋅f⋅t } ÷ √N
t=0
逆向傅立葉轉換（反函數）
N-1
x[t] = ∑ { y[f] ⋅ e𝑖⋅(2π/N)⋅f⋅t } ÷ √N
f=0
為了加快計算速度，正向傅立葉轉換經常改成不除以√N，逆向傅立葉轉換經常改成多除以√N。
N-1
y[f] = ∑ { x[t] ÷ e𝑖⋅(2π/N)⋅f⋅t }
t=0
N-1
x[t] = ∑ { y[f] ⋅ e𝑖⋅(2π/N)⋅f⋅t } ÷ N
f=0
discrete Fourier transform是線性函數
傅立葉轉換是線性函數，恰是正規正交矩陣。
ω = e𝑖⋅2π/N
⎡y0 ⎤ ⎡ ω-0⋅0 ω-0⋅1 ω-0⋅2 .. ω-0⋅(N-1) ⎤ ⎡x0 ⎤
⎢y1 ⎥ ⎢ ω-1⋅0 ω-1⋅1 ω-1⋅2 .. ω-1⋅(N-1) ⎥ ⎢x1 ⎥
⎢y2 ⎥ = ⎢ ω-2⋅0 ω-2⋅1 ω-2⋅2 .. ω-2⋅(N-1) ⎥ ⎢x2 ⎥
⎢: ⎥ ⎢ : : : : ⎥ ⎢: ⎥
⎣yN-1⎦ ⎣ ω-(N-1)⋅0 ω-(N-1)⋅1 ω-(N-1)⋅2 .. ω-(N-1)⋅(N-1) ⎦ ⎣xN-1⎦
⎡x0 ⎤ ⎡ ω0⋅0 ω0⋅1 ω0⋅2 .. ω0⋅(N-1) ⎤ ⎡y0 ⎤
⎢x1 ⎥ 1 ⎢ ω1⋅0 ω1⋅1 ω1⋅2 .. ω1⋅(N-1) ⎥ ⎢y1 ⎥
⎢x2 ⎥ = ——— ⎢ ω2⋅0 ω2⋅1 ω2⋅2 .. ω2⋅(N-1) ⎥ ⎢y2 ⎥
⎢: ⎥ N ⎢ : : : : ⎥ ⎢: ⎥
⎣xN-1⎦ ⎣ ω(N-1)⋅0 ω(N-1)⋅1 ω(N-1)⋅2 .. ω(N-1)⋅(N-1) ⎦ ⎣yN-1⎦
物理教科書規定角速度ω = 2πf。數學教科書規定第一個N次單位根ω = e𝑖⋅2π/N。符號恰巧一樣。此處是後者。
複弦波，變成離散數列，可以畫成這樣子：
傅立葉轉換的矩陣，可以畫成這樣子：
演算法（公式解）
依照公式實作，時間複雜度O(N²)。
演算法（Horner's rule）
依照公式實作，點積視作多項式，時間複雜度O(N²)。
演算法（Cooley–Tukey algorithm）
時間複雜度優於O(N²)的傅立葉轉換演算法，老人家稱作「快速傅立葉轉換fast Fourier transform, FFT」。
這裡介紹最經典的快速傅立葉轉換。公式的偶數項與奇數項分開整理，採用dynamic programming，時間複雜度O(NlogN)。
由於必須剛好對半分，所以N必須剛好是2的次方。當N不是2的次方，可在輸入數列末端補零。其功效是輸出數列的波形，重新取樣。
逆向轉換的演算法如出一轍，此處省略。
FFT
(x₀ x₁ x₂ x₃ x₄ x₅ x₆ x₇) ----> (y₀ y₁ y₂ y₃ y₄ y₅ y₆ y₇)
N = 8, ω = e-𝑖⋅2π/N 注意到ω放入了負號，讓下面的數學式子比較簡潔
y₀ = x₀ω⁰ + x₁ω⁰ + x₂ω⁰ + x₃ω⁰ + x₄ω⁰ + x₅ω⁰ + x₆ω⁰ + x₇ω⁰
= (x₀ω⁰ + x₂ω⁰ + x₄ω⁰ + x₆ω⁰) + (x₁ω⁰ + x₃ω⁰ + x₅ω⁰ + x₇ω⁰)
= (x₀ω⁰ + x₂ω⁰ + x₄ω⁰ + x₆ω⁰) + ω⁰ ⋅ (x₁ω⁰ + x₃ω⁰ + x₅ω⁰ + x₇ω⁰)
= (x₀ x₂ x₄ x₆)轉換結果第0項 + ω⁰ ⋅ (x₁ x₃ x₅ x₇)轉換結果第0項
= y偶0 + ω⁰ ⋅ y奇0
y₁ = x₀ω⁰ + x₁ω¹ + x₂ω² + x₃ω³ + x₄ω⁴ + x₅ω⁵ + x₆ω⁶ + x₇ω⁷
= (x₀ω⁰ + x₂ω² + x₄ω⁴ + x₆ω⁶) + (x₁ω¹ + x₃ω³ + x₅ω⁵ + x₇ω⁷)
= (x₀ω⁰ + x₂ω² + x₄ω⁴ + x₆ω⁶) + ω¹ ⋅ (x₁ω⁰ + x₃ω² + x₅ω⁴ + x₇ω⁶)
= (x₀υ⁰ + x₂υ¹ + x₄υ² + x₆υ³) + ω¹ ⋅ (x₁υ⁰ + x₃υ¹ + x₅υ² + x₇υ³)
= (x₀ x₂ x₄ x₆)轉換結果第1項 + ω¹ ⋅ (x₁ x₃ x₅ x₇)轉換結果第1項
= y偶1 + ω¹ ⋅ y奇1
y₂ = x₀ω⁰ + x₁ω² + x₂ω⁴ + x₃ω⁶ + x₄ω⁸ + x₅ω¹⁰ + x₆ω¹² + x₇ω¹⁴
= (x₀ω⁰ + x₂ω⁴ + x₄ω⁸ + x₆ω¹²) + (x₁ω² + x₃ω⁶ + x₅ω¹⁰ + x₇ω¹⁴)
= (x₀ω⁰ + x₂ω⁴ + x₄ω⁸ + x₆ω¹²) + ω² ⋅ (x₁ω⁰ + x₃ω⁴ + x₅ω⁸ + x₇ω¹²)
= (x₀υ⁰ + x₂υ² + x₄υ⁴ + x₆υ⁶ ) + ω² ⋅ (x₁υ⁰ + x₃υ² + x₅υ⁴ + x₇υ⁶ )
= (x₀ x₂ x₄ x₆)轉換結果第2項 + ω² ⋅ (x₁ x₃ x₅ x₇)轉換結果第2項
= y偶2 + ω² ⋅ y奇2
y₃ = x₀ω⁰ + x₁ω³ + x₂ω⁶ + x₃ω⁹ + x₄ω¹² + x₅ω¹⁵ + x₆ω¹⁸ + x₇ω²¹
= (x₀ω⁰ + x₂ω⁶ + x₄ω¹² + x₆ω¹⁸) + (x₁ω³ + x₃ω⁹ + x₅ω¹⁵ + x₇ω²¹)
= (x₀ω⁰ + x₂ω⁶ + x₄ω¹² + x₆ω¹⁸) + ω³ ⋅ (x₁ω⁰ + x₃ω⁶ + x₅ω¹² + x₇ω¹⁸)
= (x₀υ⁰ + x₂υ³ + x₄υ⁶ + x₆υ⁹ ) + ω³ ⋅ (x₁υ⁰ + x₃υ³ + x₅υ⁶ + x₇υ⁹ )
= (x₀ x₂ x₄ x₆)轉換結果第3項 + ω³ ⋅ (x₁ x₃ x₅ x₇)轉換結果第3項
= y偶3 + ω³ ⋅ y奇3
注意到 ω⁸ = 1
y₄ = x₀ω⁰ + x₁ω⁴ + x₂ω⁸ + x₃ω¹² + x₄ω¹⁶ + x₅ω²⁰ + x₆ω²⁴ + x₇ω²⁸
= (x₀ω⁰ + x₂ω⁸ + x₄ω¹⁶ + x₆ω²⁴) + (x₁ω⁴ + x₃ω¹² + x₅ω²⁰ + x₇ω²⁸)
= (x₀ω⁰ + x₂ω⁸ + x₄ω¹⁶ + x₆ω²⁴) + ω⁴ ⋅ (x₁ω⁰ + x₃ω⁸ + x₅ω¹⁶ + x₇ω²⁴)
= (x₀ω⁰ + x₂ω⁰ + x₄ω⁰ + x₆ω⁰ ) + ω⁴ ⋅ (x₁ω⁰ + x₃ω⁰ + x₅ω⁰ + x₇ω⁰ )
= (x₀ x₂ x₄ x₆)轉換結果第0項 + ω⁴ ⋅ (x₁ x₃ x₅ x₇)轉換結果第0項
= y偶0 + ω⁴ ⋅ y奇0
y₅ y₆ y₇以此類推
y₀ = y偶0 + y奇0 ⋅ ω⁰
y₁ = y偶1 + y奇1 ⋅ ω¹
y₂ = y偶2 + y奇2 ⋅ ω²
y₃ = y偶3 + y奇3 ⋅ ω³
y₄ = y偶0 + y奇0 ⋅ ω⁴
y₅ = y偶1 + y奇1 ⋅ ω⁵
y₆ = y偶2 + y奇2 ⋅ ω⁶
y₇ = y偶3 + y奇3 ⋅ ω⁷
觀察DP的遞推過程，偶數項與奇數項分開處理，索引值不連續，不易取值。預先重新排列陣列元素，符合遞推過程，減少cache miss；還可以重複使用記憶體、節省空間。
如何重新排列呢？索引值的二進位表示法，高低位數顛倒之後，恰是正確結果！
重新排列的時間複雜度是O(N)。假設高低位數顛倒是O(1)。
Hartley transform
哈特利轉換是雙射函數，輸入和輸出都是一串實數。
哈特利轉換與傅立葉轉換如出一轍，只少了虛數𝑖而已。
傅立葉轉換：
2πft 2πft -𝑖2πft/N
cos ———— - 𝑖 ⋅ sin ———— = e
N N
哈特利轉換：
2πft 2πft 2πft
cos ———— + sin ———— = cas ————
N N N
另一個哈特利轉換，比較沒人用：
2πft 2πft 2πft
cos ———— - sin ———— = cis ————
N N N
傅立葉轉換：
N-1
y[f] = ∑ { x[t] ÷ e𝑖⋅(2π/N)⋅f⋅t } ÷ √N
t=0
哈特利轉換：
N-1
y[f] = ∑ { x[t] ⋅ cas((2π/N)⋅f⋅t) } ÷ √N
t=0
由於哈特利轉換與傅立葉轉換的公式幾乎相同，所以兩者的演算法也是一一對應。這裡介紹的方法也是運用divide-and-conquer method。不一樣的是奇數項的處理方式，提出常數的步驟變複雜了。
N-1
∑ { x[t] ⋅ cas((2π/N)⋅f⋅t) }
t=1,3,5,...
N/2-1
= ∑ { x[2t+1] ⋅ cas((2π/N)⋅f⋅(2t+1)) }
t=0,1,2,...
N/2-1
= ∑ { x[2t+1] ⋅ ( cas((2π/N)⋅f⋅2t) ⋅ cos((2π/N)⋅f⋅1)
t=0,1,2,... + cas(-(2π/N)⋅f⋅2t) ⋅ sin((2π/N)⋅f⋅1) ) }
N/2-1
= ∑ { x[2t+1] ⋅ ( cas((2π/(N/2))⋅f⋅t) ⋅ cos((2π/N)⋅f⋅1)
t=0,1,2,... + cas(-(2π/(N/2))⋅f⋅t) ⋅ sin((2π/N)⋅f⋅1) ) }
N/2-1
= ∑ { x[2t+1] ⋅ cas( (2π/(N/2))⋅f⋅t) } ⋅ cos((2π/N)⋅f⋅1)
t=0,1,2,...
N/2-1
+ ∑ { x[2t+1] ⋅ cas(-(2π/(N/2))⋅f⋅t) } ⋅ sin((2π/N)⋅f⋅1)
t=0,1,2,...
= y奇[f] ⋅ cos((2π/N)⋅f⋅1) + y奇[-f] ⋅ sin((2π/N)⋅f⋅1)
θ = 2π / N
y0 = y偶0 + y奇0 ⋅ cos0θ + y奇0 ⋅ sin0θ
y1 = y偶1 + y奇1 ⋅ cos1θ + y奇3 ⋅ sin1θ
y2 = y偶2 + y奇2 ⋅ cos2θ + y奇2 ⋅ sin2θ
y3 = y偶3 + y奇3 ⋅ cos3θ + y奇1 ⋅ sin3θ
y4 = y偶0 + y奇0 ⋅ cos4θ + y奇0 ⋅ sin4θ
y5 = y偶1 + y奇1 ⋅ cos5θ + y奇3 ⋅ sin5θ
y6 = y偶2 + y奇2 ⋅ cos6θ + y奇2 ⋅ sin6θ
y7 = y偶3 + y奇3 ⋅ cos7θ + y奇1 ⋅ sin7θ
哈特利轉換的輸出，可以調整成傅立葉轉換的輸出，O(N)：
實數運算比複數運算簡單，哈特利轉換比傅立葉轉換迅速。
聲音處理，訊號都是實數，習慣採用哈特利轉換，再把結果調整成傅立葉轉換。
2D discrete Fourier transform
傅立葉轉換可以推廣到高維度。
二維傅立葉轉換，輸入輸出都是一個N×N複數方陣。輸入方陣，分別除以N×N種二維複弦波，得到N×N個輸出數值，形成輸出方陣。（由於二維複弦波很難畫，以下改畫二維實弦波。）
依照公式實作，時間複雜度O(N⁴)。高速演算法是每一橫條各自傅立葉轉換，然後每一直條各自傅立葉轉換，時間複雜度O(NNlogN + NNlogN) = O(N²logN)。
