FFT: thuật toán thay đổi thế giới

Mở đầu: cùng một câu trả lời, nhanh hơn hàng trăm lần

Bài 5 để lại một con số đáng lo: DFT trực tiếp cho N=4096 mẫu (một khung phổ âm thanh phổ biến) mất hàng trăm mili-giây để chạy trên chính máy của bạn — và với âm thanh thời gian thực cần xử lý hàng chục khung mỗi giây, đó là "chậm không thể chấp nhận được". Nhưng DFT và FFT cho ra ĐÚNG CÙNG MỘT kết quả — không phải một phép tính gần đúng nhanh hơn, mà là chính xác từng con số.

Bài này mổ xẻ FFT (Fast Fourier Transform) — thuật toán được nhiều nhà khoa học máy tính gọi thẳng là "quan trọng nhất thế kỷ 20" — từ ý tưởng chia-để-trị, qua cấu trúc butterfly, đến cài đặt lặp thật sự chạy được, VERIFY từng con số khớp DFT, rồi đo tốc độ thật trên máy của chính bạn.

1. Chia để trị: ý tưởng làm nên tốc độ

Quan sát mấu chốt: DFT $N$ điểm có thể tách thành 2 DFT $N/2$ điểm — một trên các mẫu CHẴN, một trên các mẫu LẺ — cộng thêm đúng $N$ phép "vá" bằng hệ số xoay gọi là twiddle factor ($W_N^k = e^{-j2\pi k/N}$). Áp dụng ý tưởng này ĐỆ QUY — tách $N/2$ thành $N/4$, rồi $N/8$... tới khi chỉ còn 1 mẫu (DFT của 1 điểm chính là bản thân nó) — biến độ phức tạp từ $O(N^2)$ thành $O(N \log N)$.

Càng tăng $N$, khoảng cách càng nới rộng — đây chính là lý do FFT không chỉ "nhanh hơn" mà thay đổi hẳn những gì khả thi: xử lý âm thanh/hình ảnh/tín hiệu vô tuyến thời gian thực chỉ tồn tại được nhờ khoảng cách này.

2. Butterfly & bit-reversal

Mỗi bước "vá" 2 kết quả con (chẵn/lẻ) lại với nhau có dạng cố định gọi là butterfly: 2 đầu vào, 2 đầu ra, đúng 1 phép nhân twiddle và 1 cặp cộng/trừ:

$$X_{even}[k] + W_N^k \cdot X_{odd}[k] \quad \text{ và } \quad X_{even}[k] - W_N^k \cdot X_{odd}[k]$$

Khi chuyển từ đệ quy sang cài đặt LẶP tại chỗ (Mục 3), mảng đầu vào phải được sắp xếp lại THEO THỨ TỰ ĐẢO BIT trước khi chạy các stage butterfly. Đây KHÔNG phải phép thuật — là hệ quả trực tiếp của việc tách chẵn/lẻ liên tục: mỗi lần tách là một bit quyết định "đi vào nửa chẵn hay nửa lẻ", lặp lại $\log_2 N$ lần từ bit thấp nhất tạo ra đúng thứ tự đảo bit đó. Với $N=8$ (3 bit): mẫu thứ 1 (001) chuyển tới vị trí 4 (100), mẫu thứ 3 (011) chuyển tới vị trí 6 (110) — verify bằng số thật ở Mục 4.

3. Cài đặt radix-2 in-place: từ đệ quy sang 3 vòng lặp

Đệ quy dễ hiểu nhưng tốn ngăn xếp gọi hàm. Cài đặt thực dụng thay bằng ĐÚNG 3 vòng lặp lồng nhau chạy trên 1 mảng, không đệ quy — sau khi đã sắp xếp lại theo bit-reversal (Mục 2):

for (let stage = 1; stage <= numBits; stage++) {
  const m = 1 << stage;        // kich thuoc nhom butterfly o stage nay
  const halfM = m >> 1;
  const angleStep = (-2 * Math.PI) / m;
  for (let start = 0; start < N; start += m) {      // tung nhom
    for (let k = 0; k < halfM; k++) {                // tung cap butterfly
      const angle = angleStep * k;                    // twiddle factor
      // ... nhan twiddle, cong/tru, ghi de tai cho (in-place) ...
    }
  }
}

VERIFY bắt buộc (kỷ luật xuyên suốt series): FFT phải cho ra ĐÚNG kết quả DFT trực tiếp, từng con số một, sai số dưới $10^{-10}$ — không phải "gần đúng", mà chính xác tới sai số làm tròn dấu phẩy động. Số đo thật từ self-test (node dsp-core.js): với $N=8$ và $N=16$ (tín hiệu sine thật), sai số lớn nhất giữa fft() và dft() đều dưới $10^{-10}$.

4. FFT nhanh đến đâu trên máy BẠN?

Lý thuyết là một chuyện — đo thật trên chính máy bạn là chuyện khác. Bấm nút bên dưới để chạy DFT trực tiếp và FFT radix-2 trên cùng một tín hiệu $N=8192$ mẫu, cùng import từ DSPJS thật:

5. Thực hành: Butterfly Diagram tương tác

Demo dưới đây chạy đúng thuật toán 3-vòng-lặp của Mục 3 trên $N=8$ mẫu — bấm "Stage tiếp theo" để xem dữ liệu chảy qua từng lớp butterfly, rồi so khớp kết quả cuối với fft() thật của DSPJS:

import { fft, bitReverse } from './dsp-core.js';

const N = 8, numBits = 3;
const inputSignal = [1, 2, 3, 4, 5, 6, 7, 8];

// Sap xep lai theo bit-reversal - dung HET logic that cua fft()
let re = new Array(N), im = new Array(N);
for (let i = 0; i < N; i++) {
  const j = bitReverse(i, numBits);
  re[j] = inputSignal[i];
  im[j] = 0;
}

// Chay tung stage khi bam nut - dung 3 vong lap giong het Muc 3
function runStage(stage) { /* ... giong dsp-core.js fft() ... */ }

// Sau stage cuoi: so sanh voi fft() THAT de xac nhan dung
const real = fft(inputSignal);
// verify: max(|real[i].re - re[i]|, |real[i].im - im[i]|) < 1e-9 voi moi i

Tóm lược

• ✅ FFT chia-để-trị: DFT $N$ điểm = 2 DFT $N/2$ điểm (chẵn/lẻ) + $N$ phép vá twiddle, đệ quy tới đáy → $O(N \log N)$ thay vì $O(N^2)$.
• ✅ Verified: $N=1024$ nhanh hơn khoảng 102 lần (1.048.576 vs 5.120 phép nhân); $N=4096$ nhanh hơn khoảng 341 lần.
• ✅ Butterfly: 2 vào, 2 ra, 1 twiddle, 1 cặp cộng/trừ — khối xây dựng lặp lại của mọi stage.
• ✅ Cài đặt lặp tại chỗ cần bit-reversal TRƯỚC, rồi đúng 3 vòng lặp (stage → nhóm → twiddle) — không đệ quy.
• ✅ Verified bắt buộc: FFT ≡ DFT từng con số, sai số <$10^{-10}$ ($N=8$ và $N=16$). FFT ném lỗi ngay nếu $N$ không phải luỹ thừa 2.
• ✅ Zero-padding nội suy phổ mượt hơn nhưng KHÔNG thêm thông tin tần số mới.
• ✅ Nền tảng của MP3/JPEG/5G-OFDM/MRI — và fast convolution (trả lời cổ tức Bài 4).

Trắc nghiệm ôn tập