FFT: быстрое преобразование

DFT: $O(N^2)$. FFT: $O(N \log N)$. При N = 1 000 000 разница в 50 000 раз. Именно FFT сделала MP3, JPEG и 5G возможными - и всё это один алгоритм, который Гаусс набросал в 1805 году и забыл.

• Whisper (OpenAI): mel-spectrogram через STFT (FFT под капотом) как input representation - без FFT речевое распознавание в реальном времени невозможно
• 5G NR: OFDM использует 4096-point FFT для каждого символа каждые 71 мкс - 14 000 FFT в секунду на одну несущую
• JPEG/HEIC: 2D DCT (родственник FFT) компрессирует каждый блок 8x8 пикселей - стандарт, которому 35 лет и он до сих пор используется везде
• GPS: кросс-корреляция через FFT для поиска спутникового сигнала - холодный старт за секунды, не минуты

Алгоритм Кули-Тьюки: разделяй и властвуй

**FFT был изобретён в 1965 году Джеймсом Кули и Джоном Тьюки. Потом историки нашли тот же алгоритм в набросках Гаусса 1805 года.** 160 лет алгоритм ждал своего часа - пока компьютеры не стали достаточно мощными, чтобы его ценность стала закономерной. Ключевая идея: DFT размера N раскладывается на два DFT размера N/2 - для чётных и нечётных индексов входного сигнала. Это рекурсия, которая превращает $O(N^2)$ в $O(N \log N)$.

**DFT напрямую:** $X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-j2\pi kn/N}$. Для каждого из N значений k нужно N умножений - итого $O(N^2)$. При N = 1 000 000 это $10^{12}$ операций. FFT сводит это к $N \log_2 N \approx 2 \times 10^7$ операций - разница в 50 000 раз.

Разложение Кули-Тьюки использует тот факт, что $e^{-j2\pi kn/N}$ - периодическая функция. Комплексный множитель $W_N^k = e^{-j2\pi k/N}$ называется twiddle factor. DFT чётных индексов ($x[0], x[2], x[4], \ldots$) и нечётных ($x[1], x[3], x[5], \ldots$) вычисляются независимо, а потом объединяются через butterfly operation: $X[k] = E[k] + W_N^k \cdot O[k]$, где $E$ и $O$ - DFT чётной и нечётной частей.

Whisper (OpenAI) использует STFT - короткое окно FFT, скользящее по аудиосигналу. Из выхода STFT строится mel-spectrogram, который и является входом трансформера. Без FFT речевое распознавание в реальном времени было бы невозможным.

Radix-2 DIT: bit-reversal и butterfly

**Radix-2 DIT (decimation-in-time) - самая распространённая форма FFT.** Требование: N должно быть степенью 2 (512, 1024, 2048...). Рекурсия раскрывается в $\log_2 N$ стадий butterfly-операций. Каждая стадия обрабатывает все N элементов - итого $N \log_2 N$ умножений вместо $N^2$. Вход нужно переставить в bit-reversal порядке: индекс n записывается в двоичном виде, биты переворачиваются - это условие правильного in-place вычисления.

**Bit-reversal пример для N=8:** индексы 0-7 в двоичном - 000,001,010,011,100,101,110,111. После переворота бит: 000,100,010,110,001,101,011,111 - то есть порядок 0,4,2,6,1,5,3,7. Именно в таком порядке поступают элементы в первую стадию butterfly. numpy.fft.fft делает это внутри автоматически.

Butterfly operation на каждой стадии соединяет пары элементов: $a' = a + W \cdot b$, $b' = a - W \cdot b$. Это in-place: $a'$ и $b'$ записываются на те же позиции, что $a$ и $b$. После $\log_2 N$ стадий массив содержит готовые коэффициенты DFT. В музыкальных приложениях именно такая структура позволяет применять reverb, EQ и pitch-shift в реальном времени - FFT блока за 23 мс (1024 сэмпла при 44100 Гц) вычисляется за микросекунды.

numpy.fft.fft принимает массив любой длины, не обязательно степень 2 - внутри используется смешанный radix алгоритм (radix-4, radix-8 для N вида 2^a * 3^b * 5^c). Для длин вида 2^n он всё равно быстрее всего. FFTW автоматически выбирает оптимальный plan под конкретное N и аппаратуру.

O(N log N): откуда берётся скорость

**Рекуррентное соотношение FFT:** $T(N) = 2T(N/2) + O(N)$. По теореме Мастера (случай 2) это даёт $T(N) = O(N \log N)$. Интерпретация: задача размера N делится на две задачи размера N/2 (рекурсия) и объединяется за $O(N)$ (butterfly stage). При N = $2^{20} \approx 10^6$ наивный DFT требует $10^{12}$ операций, FFT - $2 \times 10^7$. На современном CPU это разница между 16 минутами и 50 миллисекундами.

**5G OFDM на практике:** каждый OFDM символ использует 4096-point IFFT на стороне передатчика и 4096-point FFT на приёмнике. Символы передаются каждые 71 мкс. Это 14 000 FFT/IFFT в секунду - только для одной несущей. Без $O(N \log N)$ сложности мобильная связь пятого поколения физически невозможна.

FFTW (Fastest Fourier Transform in the West) - библиотека, которая автоматически подбирает оптимальный алгоритм под конкретное N и аппаратуру. При первом вызове FFTW запускает бенчмарк различных factorization планов и сохраняет результат в WISDOM файл. scipy.fft использует FFTW под капотом на большинстве платформ. Для N = prime (простое число) используются алгоритмы Блюштейна или Рейдера - сложность остаётся $O(N \log N)$.

GPS-приёмники используют FFT для кросс-корреляции принятого сигнала с известными PRN-кодами спутников. Без FFT поиск спутников занимал бы секунды вместо миллисекунд - холодный старт GPS был бы недопустимо медленным.

numpy.fft, scipy.fft, cuFFT: практика

**1994 год. Fraunhofer Institute. MP3 готов. Под капотом - Modified DCT (MDCT), двоюродный брат FFT. 128 kbps против 1411 kbps CD - разница в 11 раз при почти неотличимом звуке.** На практике FFT появляется везде: numpy.fft.fft для общего случая, scipy.fft.rfft для вещественных сигналов (выход в N/2+1 точках вместо N - эксплуатирует эрмитову симметрию), cuFFT на GPU для больших батчей. Mel-spectrogram pipeline Whisper: librosa.stft (FFT) -> амплитуды -> mel filterbank -> log -> нормализация -> трансформер.

**scipy.fft vs numpy.fft:** scipy.fft.fft быстрее на больших N за счёт использования FFTW или pocketfft backend. rfft принимает N вещественных чисел и возвращает N//2+1 комплексных - вдвое меньше памяти и операций. irfft восстанавливает оригинальный сигнал. Для батчей: scipy.fft.fft принимает параметр workers=-1 для многопоточности.

JPEG использует 2D DCT (дискретное косинусное преобразование) на каждом блоке 8x8 пикселей. DCT - частный случай DFT для вещественных чётных сигналов. MRI reconstruction использует 2D FFT (k-space -> image space). CT-сканеры используют filtered back-projection с FFT для свёртки. FFT - универсальный инструмент обработки сигналов, не только аудио.

Ключевые идеи

• Алгоритм Кули-Тьюки разбивает DFT на два DFT вдвое меньшего размера (divide and conquer по чётным/нечётным индексам), что превращает $O(N^2)$ в $O(N \log N)$ - без потери точности
• Radix-2 DIT требует N = степень 2, использует bit-reversal permutation для in-place вычисления через $\log_2 N$ стадий butterfly операций
• При N = $2^{20}$ FFT в 50 000 раз быстрее наивного DFT; FFTW автоматически выбирает оптимальный план под аппаратуру
• scipy.fft.rfft для вещественных сигналов возвращает N//2+1 коэффициентов, эксплуатируя эрмитову симметрию DFT

Связанные темы

FFT строится на DFT и открывает путь к оконному анализу; алгоритмически это классический divide and conquer

• DFT: дискретное преобразование Фурье — FFT - быстрый алгоритм вычисления DFT; без понимания DFT butterfly operation теряет смысл
• Оконные функции и STFT — STFT = скользящий FFT с оконной функцией; spectral leakage без окна - следствие конечности FFT
• Асимптотическая сложность алгоритмов — O(N log N) FFT - классический пример применения теоремы Мастера к divide and conquer
• Научные вычисления: numpy и scipy — numpy.fft и scipy.fft - стандартные реализации FFT; понимание API и выбора rfft/fft критично для практики

Вопросы для размышления

• Алгоритм Гаусса 1805 года и алгоритм Кули-Тьюки 1965 года математически идентичны. Что изменилось за 160 лет - и почему это изменение сделало алгоритм практически значимым?
• OFDM в 5G использует 4096-point FFT каждые 71 мкс. Какие требования это предъявляет к латентности вычислений - и почему именно $O(N \log N)$ здесь критично, а не просто 'быстрый процессор'?
• scipy.fft.rfft возвращает вдвое меньше коэффициентов для вещественного сигнала. При каких условиях стоит использовать fft вместо rfft - и теряется ли информация при использовании rfft?

Связанные уроки

• dsp-04
• dsp-06
• alg-01-big-o
• sci-05
• alg-19-divide-conquer