快速傅立葉轉換-下

2025/9/6 數學傅立葉分析訊號與系統快速計算

接續上篇文章快速傅立葉轉換-上，我們介紹了離散傅立業轉換的推導還有性質，這篇文章中我們要介紹重頭戲，Cooley–Tukey FFT 演算法。

🔸Cooley–Tukey FFT 演算法推導

現在我們DFT代表離散傅立葉算子，且假設是的倍數

離散訊號，DFT可寫成 $是偶數是奇數將的基數項與偶數項拆開與無關週期可除性是偶數是奇數$ 我們得到Cooley–Tukey FFT 演算法中第一個等式接著我們將重點放在上面，觀察當時改變變數，令只是一個變數，我們將改名成，得到：接著利用旋轉因子的性質： $週期性對陣性$ 得到所以Cooley–Tukey FFT 演算法可化簡成以下方程組透過以上推導，我們知道一個個點的傅立葉轉換，可以拆成兩個個點的傅立葉轉換組合出來，並且前 $項的公式，與後$ 的公式差一個負號，的範例(圖片引用自 1)

觀察偶數的這組DFT，繼續往下拆解

觀察偶數的偶數 的這組DFT，現在只剩下兩個元素了，是最簡單的格式

觀察以上圖案就像一隻蝴蝶，這就是蝴蝶圖(Butterfly diagram)的由來

🔸FFT演算法說明

在使用蝴蝶圖（Butterfly Diagram）來計算離散傅立葉轉換（DFT）時，我們首先需要對輸入序列進行一種特別的重新排序，這種排序稱為位反轉排列（Bit-Reversal Permutation）。

關於如何將一般數列轉換為 Bit-Reversal 排列，我們在此先不展開說明。現在假設序列長度為，其 Bit-Reversal 排列結果為：我們採用自底向上的（bottom-up）疊代方式，這種方式較為直觀

令 Bit-Reversal 排列後的初始序列為：每一層的蝴蝶圖使用展開長度，其中是當前的疊代層數。FFT一共會進行層的運算，在的情況下總共會做層

① 第層（）

這是最簡單的蝴蝶圖，每組僅有 2 個數值。此時旋轉因子為：我們對每兩個元素分成一組，每組蝴蝶運算，總共做了四組蝴蝶運算：

②第層（）

這一層每組蝴蝶包含 4 個元素，使用旋轉因子：將分為兩組，每組 4 個元素，分別進行蝴蝶運算：

第一組（）：

第二組（）：

③第層（）

這是最終一層，將整個向量視為一組蝴蝶圖進行處理。使用旋轉因子：對執行如下蝴蝶運算，產出最終的 FFT 輸出：

完成第層後，FFT 演算法的疊代過程也隨之結束。最終輸出即為，也就是序列的離散傅立葉轉換結果

🔸FFT演算法實作

用Python實作以長度是2的冪的FFT演算法

import numpy as np

def fft_butterfly(x):
    # 假設輸入 x 已經先做過 BIT Reverse 排列
    N = len(x)
    Layer = int(np.log2(N))
    x = np.array(x, dtype=complex)
    for layer in range(1, Layer + 1):
        m = 2**layer  # 每次蝴蝶圖翅膀展開的長度
        # 如果用 N=8 的範例就是從 m=2 慢慢展開到 m=4，組合出所求 y[k]
        half_m = m // 2
        omega_m = np.exp(-2j * np.pi / m)
        for k in range(0, N, m):
            # layer=1 時，取 k=0, 2, 4, 6，也就是每個最小蝴蝶 (展開長 2) 的左邊 index
            # layer=2 時，取 k=0, 4，也就是每個中蝴蝶 (展開長 4) 的左邊 index
            for j in range(half_m):
                t = omega_m**j * x[k + j + half_m]
                u = x[k + j]
                x[k + j] = u + t
                x[k + j + half_m] = u - t
    return x

我們來驗證當演算法的正確性

from numpy.fft import fft
import numpy as np
from numpy import linalg as LA

x = [complex(i, 0) for i in range(8)]
bit_reverse = [0, 4, 2, 6, 1, 5, 3, 7]
x_bit_reverse = [x[i] for i in bit_reverse]
y = fft_butterfly(x_bit_reverse)
print(y)
y2= fft(x)
print(y2)
print(LA.norm(y - y2))

執行結果如下：

[28.+0.j         -4.+9.65685425j -4.+4.j         -4.+1.65685425j
 -4.+0.j         -4.-1.65685425j -4.-4.j         -4.-9.65685425j]
[28.+0.j         -4.+9.65685425j -4.+4.j         -4.+1.65685425j
 -4.+0.j         -4.-1.65685425j -4.-4.j         -4.-9.65685425j]
2.5121479338940403e-15

結果顯示我們的自己做的FFT 與 numpy 內建的FFT結果一致

Bit-Reversal Permutation

現在，我們可以來探討這種特殊排列是如何產生的。許多網路上文章會直接告訴你：「只要將數字轉成二進位，再把位元順序整個反轉，這樣就得到Bit-Reversal排列了」，接著配個範例說明。雖然這樣的說法快速直觀，很難真正理解為什麼要這麼做。

事實上，這個排列並不是憑空出現，而是伴隨著蝴蝶圖中「偶數在上、奇數在下」的分組方式產生的。在蝴蝶圖演算法中，每一層的運算會將資料分成偶數索引與奇數索引兩組來運算。分組方式隨著層數遞迴深入，最終會將資料重新排序成一個特殊的順序，而這個順序就是所謂的 Bit-Reversal 排列。

📌我們用Top Down 的方式來觀察蝴蝶圖 N = 8 計算的方式

① 第層（）

為了得到N = 8 的離散傅立葉轉換，需要將

的數列分成 偶數索引

與奇數索引

，觀察這些指標的二進位表示法

索引類別	十進位	二進位
偶數	0	000
偶數	2	010
偶數	4	100
偶數	6	110
奇數	1	001
奇數	3	011
奇數	5	101
奇數	7	111

從上述表格發現，在層次，事實上是透過二進位的右邊數過來第1個bit來分辨 偶數索引與奇數索引

②第層（）

將時候的 偶數索引與奇數索引，在細拆成分別的「偶數索引」與「奇數索引」。從這裡開始，現在的奇數偶數的意義，不是一般意義的奇數偶數，而是分組後數列的奇數與偶數，舉例來說，當目前的索引是時，他對應的偶數索引應該是，奇數索引是。

第一組（的 偶數索引）：

索引類別	十進位	二進位
偶數	0	000
偶數	4	100
奇數	2	010
奇數	6	110

第二組（的奇數索引）：

索引類別	十進位	二進位
偶數	1	001
偶數	5	101
奇數	3	011
奇數	7	111

從上述表格發現，在層次，事實上是透過二進位的右邊數過來第2個bit來分辨 偶數索引與奇數索引

③ 第層（）

同理，進一步把層分好的資料，依照右數過來第 2 個 bit來再分組：

第一組（的 偶數索引）

索引類別	十進位	二進位
偶數	0	000
奇數	4	100

第二組（的 奇數索引）

索引類別	十進位	二進位
偶數	2	010
奇數	6	110

第三組（的 偶數索引）

索引類別	十進位	二進位
偶數	1	001
奇數	5	101

仔細觀察就會發現，每往下一層遞迴，其實就是依據二進位表示中「右邊數來第個 bit」進行分組，順序如下：

在時，看 右邊數來第 1 個 bit
在時，看 右邊數來第 2 個 bit
在時，看 右邊數來第 3 個 bit

以上觀察可以整理如下圖：

蝴蝶圖的奧祕就在於：每次蝴蝶圖運算進行「偶數在上、奇數在下」分組時，參照的就是該數字二進位表示的「右邊數來第個 bit」進行分組。

我們來看上圖中的例子：

Bit-Reversal 排列中的第 0 個元素，由「偶偶偶」給出，對應二進位，也就是十進位的
第 4 個元素，由「奇偶偶」給出，對應二進位，也就是十進位的
第 6 個元素，由「奇奇偶」給出，對應二進位，也就是十進位的
依此類推

📌 綜合以上說明，我們可以總結：當要找到Bit-Reversal 排列中的第個元素時，只需要將轉成二進位，接著將位元順序反轉，最後再轉回十進位，這個值就是在 Bit-Reversal 排列中對應的位置。

而這個特殊排列，其實正是蝴蝶圖層層分組結構自然產生的副產品，而非刻意安排的編號方式

用Python實作以長度是2的冪的Bit-Reversal 排列

import math
def bit_reverse_permutation(arr):
    """
    對輸入陣列 arr 做 bit-reversal permutation（位元反轉排列）
    """
    
    n = len(arr)
    result = []
    
    for i in range(n):
        # 計算 i 的 bit-reverse 索引
        rev = int('{:0{width}b}'.format(i, width=N)[::-1], 2)
        result.append(arr[rev])
    
    return result

再次驗證FFT的實作

from numpy.fft import fft
import numpy as np
from numpy import linalg as LA


x = [complex(i, 0) for i in range(16)]
x_bit_reverse = bit_reverse_permutation(x)
y = fft_butterfly(x_bit_reverse)
print(y)
y2= fft(x)
print(y2)
print(LA.norm(y - y2))

執行結果如下：

[120. +0.j          -8.+40.21871594j  -8.+19.3137085j   -8.+11.9728461j
  -8. +8.j          -8. +5.3454291j   -8. +3.3137085j   -8. +1.59129894j
  -8. +0.j          -8. -1.59129894j  -8. -3.3137085j   -8. -5.3454291j
  -8. -8.j          -8.-11.9728461j   -8.-19.3137085j   -8.-40.21871594j]
[120. +0.j          -8.+40.21871594j  -8.+19.3137085j   -8.+11.9728461j
  -8. +8.j          -8. +5.3454291j   -8. +3.3137085j   -8. +1.59129894j
  -8. +0.j          -8. -1.59129894j  -8. -3.3137085j   -8. -5.3454291j
  -8. -8.j          -8.-11.9728461j   -8.-19.3137085j   -8.-40.21871594j]
1.6687357476335705e-14

🔸時間複雜度分析

📌我們先討論直接計算DFT的時間複雜度

回想DFT計算公式以下是一個直接計算DFT的程式

import numpy as np
def DFT_slow(x):
    """
    x : 輸入訊號 (list 或 numpy array)，長度 N
    return : DFT 結果 (numpy array)，長度 N
    """
    x = np.asarray(x, dtype=complex)   # 確保輸入是複數型別
    N = x.shape[0]
    X = np.zeros(N, dtype=complex)     # 輸出 DFT 結果
    
    for k in range(N):  # 頻率索引
        X[k] = 0
        for n in range(N):  # 時間索引
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

接著我們分析這個演算法，假設每行程式碼的執行時間是 t

透過上圖的討論，直接計算DFT的演算法複雜度為

當很大的時候，主導項是所以我們可以說，這個直接計算DFT的演算法複雜度為，更精確的寫法應該是，因為是tightest bound。

📌FFT的時間複雜度

回想FFT的分解式 $是偶數是奇數$ 也就是說，假設我們已經分別算出了 $是偶數$ 與 $是奇數$ ，我們要得到完整的長度的 DFT，只需要將這兩組結果進行簡單的「合併」即可。這個組合過程就是經典的 蝶形運算 (butterfly operation)。

假設為 FFT 的時間複雜度：

首先，我們需要遞迴計算兩個長度的 DFT，這部分花費時間為
接著，當兩個子問題計算完成後，必須將它們合併成最終的長度的結果。在這裡，每一對偶數、奇數頻率分量都要經過一次「乘法」與「加減法」，如下圖 1：

總共會有 個輸出(的項數)需要處理，因此合併過程的運算次數是 正比於 的，而且這個就是tightest bound。

所以，綜合上述兩部分討論，的遞迴關係式可以寫作套用，的定義，將(3)改寫：依據 Master theorem (大師套公式法)

,
「合併」計算量= 與相同 大師套公式法 Case 2

所以：不嚴格一點也可以寫成：

用以下程式碼比較一下在，直接計算DFT與FFT速度的差異

from numpy.fft import fft
import numpy as np
from numpy import linalg as LA
import math
import time

N = 2**12
x = [complex(i, 0) for i in range(N)]

ans = fft(x)
start_fft = time.time()
x_bit_reverse = bit_reverse_permutation(x)
y1 = fft_butterfly(x_bit_reverse)
end_fft = time.time()

start_DFT_slow = time.time()
y2 = DFT_slow(x)
end_DFT_slow = time.time()

print(f"FFT butterfly time: {end_fft - start_fft:.6f} seconds")
print(f"DFT slow time: {end_DFT_slow - start_DFT_slow:.6f} seconds")
print(f"FFT butterfly Norm difference: {LA.norm(ans - y1)}")
print(f"DFT slow Norm difference: {LA.norm(ans - y2)}")

執行結果如下：

FFT butterfly time: 0.012702 seconds
DFT slow time: 10.643429 seconds
FFT butterfly Norm difference: 1.596628055678928e-07
DFT slow Norm difference: 8.190390264959864e-06

看的出來透過FFT演算法來加速，可以提升計算DFT的效率到1000倍以上

🔸FFT的應用：快速離散卷積運算

在上一篇快速傅立葉轉換-上中，我們介紹了 離散卷積定理。根據該定理，若要在時域中計算離散卷積，可以改為：先將兩個序列轉換到頻域，在頻域中進行逐點乘積，最後再透過 iDFT 將結果轉換回時域。

直接在時域計算兩個長度為的序列卷積，其時間複雜度為：然而，若改用 FFT，整體流程會大幅加速：

FFT 轉換：將兩個序列分別做 FFT，計算量為
頻域逐點相乘：僅需
iDFT 還原：透過與 FFT 相同的快速演算法完成，計算量為事實上，iDFT 也能透過 FFT 演算法加速，其原理與 FFT 推導幾乎相同，因此這裡筆者不再贅述，僅點出最終結論。

綜合上述三個步驟，總時間複雜度為由於的主導項是對數線性，最終可化簡為：因此，透過 FFT 加速運算的方法，離散卷積的計算量從原本的平方等級，有效降到。

📌在實作加速離散卷積之前，我們需要先基於蝴蝶圖(Butterfly diagram)來實作iFFT，先前我們已經有FFT的實作，可以利用DFT與iDFT中的共軛性質來直接實作iFFT

回想DFT與iDFT的公式我們對取任意序列，做共軛、再做 DFT：對上式等號兩邊同時取共軛：其中所以注意到以上結果已經很接近iDFT公式了，我們做符號代換，並且在最前面乘上

得到基於以上推導，我們可以用原先的fft_butterfly(x)，加上共軛的操作來實作iFFT

from numpy.fft import fft
from numpy.fft import ifft
import numpy as np
from numpy import linalg as LA
import random

def ifft_butterfly(x):
    """
    使用 fft_butterfly 實作 IFFT。
    """
    N = len(x)
    # 共軛
    x_conj = np.conjugate(x)
    # FFT butterfly
    y = fft_butterfly(x_conj)
    # 取共軛後乘上1/N
    return np.conjugate(y) / N

print("測試 IFFT 蝴蝶算法的正確性:")
print("="*50)

# 生成測試訊號
test_length = 8  # 使用 2^3 長度做測試
test_signal = [complex(random.randint(1, 10), random.randint(1, 10)) for _ in range(test_length)]

# 使用 numpy 的 FFT 產生頻域數據
freq_data = fft(test_signal)
print(f"Numpy FFT 產生的頻域數據: {freq_data}")

# 使用自製的 IFFT 蝴蝶算法
freq_data_br = bit_reverse_permutation(freq_data)
butterfly_result = ifft_butterfly(freq_data_br)

# 使用 numpy 的 IFFT 作為基準
numpy_result = ifft(freq_data)

# 計算誤差
error = LA.norm(np.array(butterfly_result) - numpy_result)

print(f"\n蝴蝶算法 IFFT 結果: {butterfly_result}")
print(f"Numpy IFFT 結果: {numpy_result}")
print(f"\n誤差分析:")
print(f"蝴蝶算法與 Numpy IFFT 的差異 (Norm): {error}")
print("="*50 + "\n")

執行結果如下：

==================================================
Numpy FFT 產生的頻域數據: [ 43.        +45.j         -13.94974747 -4.70710678j
 -10.         -4.j           9.70710678-10.46446609j
  -9.         -1.j          -4.05025253 -3.29289322j
 -16.         +4.j           8.29289322-17.53553391j]

蝴蝶算法 IFFT 結果: [ 1. +1.j  6. +5.j  5. +1.j  6. +5.j  1.+10.j  9. +8.j 10.+10.j  5. +5.j]
Numpy IFFT 結果: [ 1. +1.j  6. +5.j  5. +1.j  6. +5.j  1.+10.j  9. +8.j 10.+10.j  5. +5.j]

誤差分析:
蝴蝶算法與 Numpy IFFT 的差異 (Norm): 1.2560739669470201e-15
==================================================

📌有幾個程式實作上有個細節需要特別注意：

由於我們實作的 FFT 僅支援長度為的陣列，因此在計算 convolution 之前，必須先將兩個輸入序列zero padding至大於等於的最小的冪次長度。
因為我們有將兩個陣列都進行zero padding，也因為我們僅有實作長度的FFT 與 iFFT，所以再比較結果的時候，我們需要截斷後面多餘的

接著我們用以下程式，來比較一下FFT加速離散卷積運算的效果

from numpy.fft import fft
from numpy.fft import ifft
import numpy as np
from numpy import linalg as LA
import math
import time
import random


def discrete_convolution(f, g):
    """
    計算兩個一維向量 f 和 g 的離散卷積。
    回傳長度為 len(f) + len(g) - 1 的 numpy array。
    """
    f = np.asarray(f)
    g = np.asarray(g)
    N = len(f)
    M = len(g)
    result = np.zeros(N + M - 1, dtype=complex)
    for n in range(N + M - 1):
        for k in range(N):
            if 0 <= n - k < M:
                result[n] += f[k] * g[n - k]
    return result

# Parameters
N = 12  # N > 2 so 2**N > 2**2
signal_length = 2**N
print(f"Signal length: {signal_length}")

# Generate random input signals
f = [complex(random.randint(1, 100), 0) for _ in range(signal_length)]
g = [complex(random.randint(1, 100), 0) for _ in range(signal_length)]

# Zero-padding to next power of two (for convolution)
pad_length = 2 * signal_length
f_padded = f + [0] * (pad_length - len(f))
g_padded = g + [0] * (pad_length - len(g))
print(f"Zero-padded length: {pad_length}")

# FFT-based convolution
start_fft = time.time()
# Bit-reverse permutation for custom FFT
f_padded_br = bit_reverse_permutation(f_padded)
g_padded_br = bit_reverse_permutation(g_padded)
# FFT butterfly
F = fft_butterfly(f_padded_br)
G = fft_butterfly(g_padded_br)
FG = [F[i] * G[i] for i in range(len(F))]
# Bit-reverse permutation for inverse FFT
FG_br = bit_reverse_permutation(FG)
conv_fft = ifft_butterfly(FG_br)
end_fft = time.time()

# Direct convolution (slow)
start_conv = time.time()
conv_slow = discrete_convolution(f, g)
end_conv = time.time()

# Results and comparison
print(f"FFT butterfly convolution time: {end_fft - start_fft:.6f} seconds")
print(f"Direct convolution time: {end_conv - start_conv:.6f} seconds")
# Compare only the first len(conv_slow) elements
norm_diff = LA.norm(conv_fft[:len(conv_slow)] - conv_slow)
print(f"Norm difference between FFT and direct convolution: {norm_diff}")

執行結果如下：

Signal length: 4096
Zero-padded length: 8192
FFT butterfly convolution time: 0.072530 seconds
Direct convolution time: 1.808187 seconds
Norm difference between FFT and direct convolution: 1.0590499027122327e-05

透過FFT演算法來離散捲積，可以大幅提升計算離散捲積的速度，而且4096維的離散捲積，其結果僅的誤差，這些誤差來自浮點數的精度問題

🔸FFT的應用：數字乘法

回想小時候計算直式的乘法的過程

為了更清楚地看出直式乘法與離散卷積的關係，我們可以用兩個陣列來表示數字：

根據離散卷積的公式：觀察上圖右側紫色框起來的部分，可以發現這正是的離散卷積係數：也就是說，離散卷積的計算過程，其實正好對應到直式乘法中分別相乘與「每行加總」的步驟，如下圖所示

✨在實作上還有另一個需要注意的地方：我們在紙筆計算直式乘法時，數字的書寫方式是由高位數到低位數（從左到右）。例如，數字 1234 代表的是然而，若我們直接把它存成陣列 [1, 2, 3, 4]，就會出現一個問題：我們沒辦法第一時間就確定最前面的 1 是代表還是，因為陣列的索引是由小到大遞增的，並沒有天然對應到「位數大小」。

為了解決這個對齊問題，我們需要先將數字反轉，把最低位（個位數）對應到陣列索引。以 1234 為例，我們會把它表示成這樣一來，陣列的第個元素對應，第個元素對應，以此類推。如此一來，當我們進行卷積運算時，就能保證運算的索引順序與位數的意義一致，確保最終結果正確模擬直式乘法的加總過程。

以下是一個利用 FFT（快速傅立葉轉換） 來進行整數乘法的範例：

from numpy.fft import fft
from numpy.fft import ifft
import numpy as np
from numpy import linalg as LA
import math
import time
import random


def int_to_digits(n):
    """
    將整數轉換為數字陣列（最低位在前）。
    """
    return [int(d) for d in str(n)][::-1]

def digits_to_int(digits):
    """
    將數字陣列（最低位在前）轉回整數。
    """
    return int(''.join(str(d) for d in digits[::-1]))

def handle_carries(arr):
    """
    處理數字陣列中的進位，產生正確的整數數字。
    """
    arr = np.round(np.real(arr)).astype(int)
    carry = 0
    result = []
    for v in arr:
        total = v + carry
        result.append(total % 10)
        carry = total // 10
    while carry > 0:
        result.append(carry % 10)
        carry //= 10
    # Remove leading zeros
    while len(result) > 1 and result[-1] == 0:
        result.pop()
    return result


# 整數乘法範例
a = 1234
b = 5678
print(f"\n整數乘法範例: {a} * {b}")
a_digits = int_to_digits(a)
b_digits = int_to_digits(b)
n = len(a_digits) + len(b_digits) - 1
N = int(np.ceil(np.log2(n)))
pad_length = 2**N
# 補零到最近的 2 的次方
a_padded = a_digits + [0] * (pad_length - len(a_digits))
b_padded = b_digits + [0] * (pad_length - len(b_digits))

# 使用 FFT 蝴蝶算法做卷積（乘法）並計算時間
start_fft_mul = time.time()
a_br = bit_reverse_permutation(a_padded)
b_br = bit_reverse_permutation(b_padded)
Fa = fft_butterfly(a_br)
Fb = fft_butterfly(b_br)
Fc = [Fa[i] * Fb[i] for i in range(len(Fa))]
Fc_br = bit_reverse_permutation(Fc)
c_fft = ifft_butterfly(Fc_br)
c_digits_fft = handle_carries(c_fft)
product_fft = digits_to_int(c_digits_fft)
end_fft_mul = time.time()

# 基準：直接卷積並計算時間
start_slow_mul = time.time()
c_slow = discrete_convolution(a_digits, b_digits)
c_digits_slow = handle_carries(c_slow)
product_slow = digits_to_int(c_digits_slow)
end_slow_mul = time.time()

print(f"FFT 蝴蝶算法結果: {product_fft}")
print(f"FFT 蝴蝶算法乘法時間: {end_fft_mul - start_fft_mul:.6f} 秒")
print(f"直接卷積結果: {product_slow}")
print(f"直接卷積乘法時間: {end_slow_mul - start_slow_mul:.6f} 秒")
print(f"正確答案: {a * b}")
print(f"FFT 與直接卷積的數字差異 Norm: {LA.norm(np.array(c_digits_fft)[:len(c_digits_slow)] - np.array(c_digits_slow))}")

執行結果如下：

整數乘法範例: 1234 * 5678        
FFT 蝴蝶算法結果: 7006652        
FFT 蝴蝶算法乘法時間: 0.000105 秒
直接卷積結果: 7006652
直接卷積乘法時間: 0.000020 秒
正確答案: 7006652
FFT 與直接卷積的數字差異 Norm: 0.0

從這個小範例可以發現，FFT 加速的效果反而不如直接的暴力卷積。由於數字規模較小，FFT 本身的額外運算開銷大於加速效果，因此運行時間反而更長。

為了更清楚觀察兩種方法的差異，我們隨機生成兩個長度為的整數，並比較不同乘法方式的效能：

# 隨機生成長度為 2**N 的兩個大整數相乘範例
N = 11 # 可自行調整 N
num_digits = 2**N
# 產生隨機 digit array（每位 0~9）
rand_digits_a = [random.randint(0, 9) for _ in range(num_digits)]
rand_digits_b = [random.randint(0, 9) for _ in range(num_digits)]
# 轉成整數
rand_a = digits_to_int(rand_digits_a)
rand_b = digits_to_int(rand_digits_b)
print(f"\n隨機大整數乘法範例: 長度 {num_digits} 位")
print(f"A = {rand_a}")
print(f"B = {rand_b}")
start_python_built_in_multiplication = time.time()
res = rand_a * rand_b
end_python_built_in_multiplication = time.time()
print(f"Python 內建乘法結果: {res}")

# digit array 準備
a_digits = int_to_digits(rand_a)
b_digits = int_to_digits(rand_b)
n = len(a_digits) + len(b_digits) - 1
pad_length = 2**int(np.ceil(np.log2(n)))
a_padded = a_digits + [0] * (pad_length - len(a_digits))
b_padded = b_digits + [0] * (pad_length - len(b_digits))

# FFT 乘法
start_fft = time.time()
a_br = bit_reverse_permutation(a_padded)
b_br = bit_reverse_permutation(b_padded)
Fa = fft_butterfly(a_br)
Fb = fft_butterfly(b_br)
Fc = [Fa[i] * Fb[i] for i in range(len(Fa))]
Fc_br = bit_reverse_permutation(Fc)
c_fft = ifft_butterfly(Fc_br)
c_digits_fft = handle_carries(c_fft)
product_fft = digits_to_int(c_digits_fft)
end_fft = time.time()

# 直接卷積
start_conv = time.time()
c_slow = discrete_convolution(a_digits, b_digits)
c_digits_slow = handle_carries(c_slow)
product_slow = digits_to_int(c_digits_slow)
end_conv = time.time()

# 比較結果
print(f"\nFFT 蝴蝶算法結果: {product_fft}")
print(f"直接卷積結果: {product_slow}")

min_len = min(len(c_slow), len(c_fft))
print(f"FFT 蝴蝶算法與直接卷積Norm差異 : {LA.norm(np.array(c_slow[:min_len]) - np.array(c_fft[:min_len]))} (直接卷積)")

print(f"FFT 蝴蝶算法時間: {end_fft - start_fft:.6f} 秒")
print(f"直接卷積時間: {end_conv - start_conv:.6f} 秒")

print(f"Python 內建乘法時間: {end_python_built_in_multiplication - start_python_built_in_multiplication:.6f} 秒")

執行結果如下：

隨機大整數乘法範例: 長度 2048 位
A =  
#省略
B =  
#省略
Python 內建乘法結果:   
#省略
FFT 蝴蝶算法結果:
#省略
直接卷積結果:
#省略
FFT 蝴蝶算法與直接卷積Norm差異 : 5.7675031363963895e-08 (直接卷積)
FFT 蝴蝶算法時間: 0.033830 秒
直接卷積時間: 0.909085 秒
Python 內建乘法時間: 0.000016 秒

在這個例子中可以明顯看到：

FFT 蝴蝶算法 已經大幅快於直接卷積（速度差超過 20 倍）。
不過，因為 FFT 過程中涉及浮點數運算，會產生微小的精度誤差，導致最後需要額外處理進位與四捨五入，Norm 差異大約在等級。這也是為什麼 FFT 並不是專為大數精確運算設計的演算法。
事實上，若要進行更精確的大數乘法，有其他更合適的轉換與演算法。

最後值得一提的是，Python 的便利性：

在 Python 中，當整數超過一般機器整數範圍時，語言本身會自動切換到「大數運算」（arbitrary-precision integer）。
使用者不需要額外處理，直接寫 a * b 就能正確計算任意大整數。
這與 Java 或 C++ 不同，在那些語言裡，內建整數型別（例如 int、long）有固定大小，若要進行大數運算，需要額外引入 BigInteger（Java）或使用第三方大數庫（C++），開發體驗相對麻煩。

Reference

庫利-圖基快速傅立葉變換演算法

convolution - 演算法筆記

LOADING...

快速傅立葉轉換-下

🔸Cooley–Tukey FFT 演算法推導

🔸FFT演算法說明

① 第 層（）

②第 層（）

第一組（）：

第二組（）：

③第 層（）

🔸FFT演算法實作

用Python實作以長度是2的冪的FFT演算法

Bit-Reversal Permutation

以上觀察可以整理如下圖：

用Python實作以長度是2的冪的Bit-Reversal 排列

🔸時間複雜度分析

🔸FFT的應用：快速離散卷積運算

🔸FFT的應用：數字乘法

Reference

① 第層（）

②第層（）

③第層（）