Преобразование Фурье

Каждый сигнал — это сумма синусоид. Преобразование Фурье показывает, каких именно.

Тип: Реализация Язык: Python Предварительные требования: Фаза 1, уроки 01–04, 19 (комплексные числа) Время: ~90 минут

Цели обучения

• Реализовать DFT с нуля и проверить её результат против FFT Кули—Тьюки (Cooley-Tukey) со сложностью O(N log N)
• Интерпретировать частотные коэффициенты: извлекать амплитуду, фазу и спектр мощности сигнала
• Применять теорему о свёртке, чтобы выполнять свёртку умножением в FFT
• Связать частотную декомпозицию Фурье с позиционными кодировками трансформеров и свёрточными слоями CNN

Проблема

Аудиозапись — это последовательность измерений давления во времени. Цена акции — последовательность значений по дням. Изображение — сетка интенсивностей пикселей в пространстве. Всё это данные во временной области (или пространственной области). Вы видите значения, меняющиеся вдоль некоторого индекса.

Но многие закономерности невидимы во временной области. Этот аудиосигнал — чистый тон или аккорд? Есть ли у цены акции недельный цикл? Есть ли у изображения повторяющаяся текстура? Эти вопросы относятся к частотному содержанию, а временная область его скрывает.

Преобразование Фурье переводит данные из временной области в частотную. Оно берёт сигнал и раскладывает его на синусоиды разных частот. У каждой синусоиды есть амплитуда (насколько она сильна) и фаза (с какого положения она начинается). Преобразование Фурье сообщает и то, и другое.

Это важно для ML, потому что частотное мышление встречается повсюду. Свёрточные нейронные сети выполняют свёртку, а она является умножением в частотной области. Позиционные кодировки трансформеров используют частотную декомпозицию для представления позиции. Аудиомодели (распознавание речи, генерация музыки) работают со спектрограммами — частотными представлениями звука. Модели временных рядов ищут периодические паттерны. Понимание преобразования Фурье даёт вам словарь для работы со всем этим.

Концепция

Определение DFT

Для N отсчётов x[0], x[1], …, x[N-1] дискретное преобразование Фурье (Discrete Fourier Transform, DFT) создаёт N частотных коэффициентов X[0], X[1], …, X[N-1]:

X[k] = sum_{n=0}^{N-1} x[n] * e^(-2*pi*i*k*n/N)

for k = 0, 1, ..., N-1

Каждый X[k] — комплексное число. Его модуль |X[k]| задаёт амплитуду частоты k. Его фазовый угол angle(X[k]) задаёт сдвиг фазы этой частоты.

Ключевая идея: e^(-2*pi*i*k*n/N) — это вращающийся фазор на частоте k. DFT вычисляет корреляцию между сигналом и каждой из N равномерно расположенных частот. Если в сигнале есть энергия на частоте k, корреляция велика. Если нет — она близка к нулю.

Что означает каждый коэффициент

X[0]: DC-компонента. Это сумма всех отсчётов — величина, пропорциональная среднему. Она представляет постоянное смещение сигнала (нулевую частоту).

X[0] = sum_{n=0}^{N-1} x[n] * e^0 = sum of all samples

X[k] для 1 <= k <= N/2: положительные частоты. X[k] представляет частоту k циклов на N отсчётов. Чем больше k, тем выше частота (быстрее колебание).

X[N/2]: частота Найквиста. Это наивысшая частота, которую можно представить N отсчётами. Выше неё возникает алиасинг — высокие частоты маскируются под низкие.

X[k] для N/2 < k < N: отрицательные частоты. Для вещественных сигналов X[N-k] = conj(X[k]). Отрицательные частоты являются зеркальными отражениями положительных. Поэтому полезная информация находится в первых N/2 + 1 коэффициентах.

Обратное DFT

Обратное DFT восстанавливает исходный сигнал из его частотных коэффициентов:

x[n] = (1/N) * sum_{k=0}^{N-1} X[k] * e^(2*pi*i*k*n/N)

for n = 0, 1, ..., N-1

От прямого DFT оно отличается только двумя вещами: знак в показателе степени положителен (а не отрицателен), и есть множитель нормализации 1/N.

Обратное DFT даёт идеальное восстановление. Информация не теряется. Можно перейти из временной области в частотную и обратно без какой-либо ошибки. DFT — это смена базиса: оно выражает ту же информацию в другой системе координат.

FFT: как ускорить вычисление

Определённое выше DFT имеет сложность O(N^2): для каждого из N выходных коэффициентов вы суммируете N входных отсчётов. При N = 1 миллион это 10^12 операций.

Быстрое преобразование Фурье (Fast Fourier Transform, FFT) вычисляет тот же результат за O(N log N). При N = 1 миллион это около 20 миллионов операций вместо триллиона. Именно это делает частотный анализ практичным.

Алгоритм Кули—Тьюки (наиболее распространённый FFT) использует подход «разделяй и властвуй»:

1. Разделить сигнал на отсчёты с чётными и нечётными индексами.
2. Рекурсивно вычислить DFT каждой половины.
3. Объединить два DFT половинного размера с помощью «twiddle-факторов» e^(-2pii*k/N).

X[k] = E[k] + e^(-2*pi*i*k/N) * O[k]          for k = 0, ..., N/2 - 1
X[k + N/2] = E[k] - e^(-2*pi*i*k/N) * O[k]    for k = 0, ..., N/2 - 1

where E = DFT of even-indexed samples
      O = DFT of odd-indexed samples

Симметрия означает, что каждый уровень рекурсии выполняет O(N) работы, а уровней log2(N). Итого: O(N log N).

FFT требует, чтобы длина сигнала была степенью 2. На практике сигналы дополняют нулями до следующей степени 2.

Спектральный анализ

Спектр мощности — это |X[k]|^2, то есть квадрат модуля каждого частотного коэффициента. Он показывает, сколько энергии приходится на каждую частоту.

Спектр фаз — это angle(X[k]), то есть сдвиг фазы каждой частоты. В большинстве задач анализа вас интересует спектр мощности, а фазой можно пренебречь.

Power at frequency k:  P[k] = |X[k]|^2 = X[k].real^2 + X[k].imag^2
Phase at frequency k:  phi[k] = atan2(X[k].imag, X[k].real)

Частотное разрешение

Частотное разрешение DFT зависит от числа отсчётов N и частоты дискретизации fs.

Frequency of bin k:      f_k = k * fs / N
Frequency resolution:    delta_f = fs / N
Maximum frequency:       f_max = fs / 2  (Nyquist)

Чтобы различить две близко расположенные частоты, нужно больше отсчётов. Чтобы захватывать высокие частоты, нужна более высокая частота дискретизации.

Теорема о свёртке

Это один из важнейших результатов обработки сигналов, напрямую относящийся к CNN.

Свёртка во временной области равна поточечному умножению в частотной области.

x * h = IFFT(FFT(x) . FFT(h))

where * is convolution and . is element-wise multiplication

Почему это важно:

• Прямая свёртка двух сигналов длины N и M требует O(N*M) операций.
• Свёртка через FFT требует O(N log N): преобразовать оба сигнала, перемножить, преобразовать обратно.
• Для больших ядер свёртка через FFT заметно быстрее.
• Именно это происходит в свёрточных слоях с большими рецептивными полями.

Замечание: DFT вычисляет циклическую свёртку (сигнал заворачивается по кругу). Для линейной свёртки (без заворачивания) дополните оба сигнала нулями до длины N + M - 1 до вычисления.

Оконные функции

DFT предполагает, что сигнал периодичен: оно рассматривает N отсчётов как один период бесконечно повторяющегося сигнала. Если сигнал не начинается и не заканчивается на одном значении, на границе возникает разрыв, который проявляется как ложное высокочастотное содержимое. Это называется спектральной утечкой.

Оконная функция уменьшает утечку, сглаживая сигнал к нулю на обоих концах перед вычислением DFT.

Распространённые окна:

Окно	Форма	Ширина главного лепестка	Уровень боковых лепестков	Сценарий использования
Прямоугольное	Плоское (без окна)	Самая узкая	Самый высокий (-13 dB)	Когда сигнал точно периодичен на N отсчётах
Ханна	Приподнятый косинус	Умеренная	Низкий (-31 dB)	Спектральный анализ общего назначения
Хэмминга	Модифицированный косинус	Умеренная	Более низкий (-42 dB)	Обработка аудио, анализ речи
Блэкмана	Тройной косинус	Широкая	Очень низкий (-58 dB)	Когда подавление боковых лепестков критично

Hann window:    w[n] = 0.5 * (1 - cos(2*pi*n / (N-1)))
Hamming window: w[n] = 0.54 - 0.46 * cos(2*pi*n / (N-1))

Примените окно, поэлементно умножив его на сигнал до DFT: X = DFT(x * w).

Свойства DFT

Свойство	Временная область	Частотная область
Линейность	ax + by	aX + bY
Сдвиг по времени	x[n - k]	X[f] * e^(-2piifk/N)
Сдвиг по частоте	x[n] * e^(2piif0n/N)	X[f - f0]
Свёртка	x * h	X * H (поточечно)
Умножение	x * h (поточечно)	X * H (циклическая свёртка, масштабированная на 1/N)
Теорема Парсеваля	sum |x[n]|^2	(1/N) * sum |X[k]|^2
Сопряжённая симметрия (вещественный вход)	x[n] вещественно	X[k] = conj(X[N-k])

Теорема Парсеваля утверждает, что суммарная энергия одинакова в обеих областях. Энергия сохраняется при преобразовании.

Связь с позиционными кодировками

Исходный Transformer использует синусоидальные позиционные кодировки:

PE(pos, 2i)   = sin(pos / 10000^(2i/d_model))
PE(pos, 2i+1) = cos(pos / 10000^(2i/d_model))

Каждая пара измерений (2i, 2i+1) колеблется на своей частоте. Частоты располагаются геометрически — от высоких (измерения 0,1) к низким (последние измерения). Это даёт каждой позиции уникальный паттерн по всем частотным полосам — подобно тому, как коэффициенты Фурье однозначно определяют сигнал.

Ключевые свойства, которые это даёт:

• Уникальность: у двух позиций нет одинаковой кодировки.
• Ограниченные значения: sin и cos всегда находятся в [-1, 1].
• Относительная позиция: кодировка позиции p+k может быть выражена как линейная функция кодировки позиции p. Модель может научиться обращать внимание на относительные позиции.

Связь с CNN

Свёрточный слой применяет обучаемый фильтр (ядро) к входу, сдвигая его по сигналу или изображению. Математически это операция свёртки.

Согласно теореме о свёртке, это эквивалентно следующим действиям:

1. Вычислить FFT входа.
2. Вычислить FFT ядра.
3. Перемножить в частотной области.
4. Вычислить IFFT результата.

Стандартные реализации CNN используют прямую свёртку (она быстрее для небольших ядер 3x3). Но для больших ядер или глобальной свёртки подходы на основе FFT значительно быстрее. Некоторые архитектуры (такие как FNet) полностью заменяют attention на FFT, достигая конкурентной точности со сложностью O(N log N) вместо O(N^2).

Спектрограммы и кратковременное преобразование Фурье

Одно FFT даёт частотное содержимое всего сигнала, но ничего не сообщает о том, когда эти частоты возникают. Чирп (сигнал, частота которого растёт со временем) и аккорд (все частоты присутствуют одновременно) могут иметь одинаковый спектр модулей.

Кратковременное преобразование Фурье (Short-Time Fourier Transform, STFT) решает эту проблему, вычисляя FFT на перекрывающихся окнах сигнала. Результат — спектрограмма: двумерное представление, где по одной оси расположено время, а по другой — частота. Интенсивность в каждой точке показывает энергию на данной частоте в данный момент.

STFT procedure:
1. Choose a window size (e.g., 1024 samples)
2. Choose a hop size (e.g., 256 samples -- 75% overlap)
3. For each window position:
   a. Extract the windowed segment
   b. Apply a Hann/Hamming window
   c. Compute FFT
   d. Store the magnitude spectrum as one column of the spectrogram

Спектрограммы — стандартное входное представление для аудио ML-моделей. Модели распознавания речи (Whisper, DeepSpeech) работают с мел-спектрограммами — спектрограммами, в которых частоты отображены на шкалу mel, лучше соответствующую восприятию высоты звука человеком.

Алиасинг

Если сигнал содержит частоты выше fs/2 (частоты Найквиста), дискретизация с частотой fs создаст алиасированные копии. Сигнал 90 Hz, дискретизированный с частотой 100 Hz, выглядит точно так же, как сигнал 10 Hz. Только по отсчётам их невозможно различить.

Example:
  True signal: 90 Hz sine wave
  Sampling rate: 100 Hz
  Apparent frequency: 100 - 90 = 10 Hz

  The samples from the 90 Hz signal at 100 Hz sampling rate
  are identical to the samples from a 10 Hz signal.
  No amount of math can recover the original 90 Hz.

Поэтому аналого-цифровые преобразователи включают anti-aliasing-фильтры, которые удаляют частоты выше Найквиста перед дискретизацией. В ML алиасинг возникает при понижении разрешения карт признаков без корректной низкочастотной фильтрации; некоторые архитектуры решают это слоями pooling с подавлением алиасинга.

Дополнение нулями не повышает разрешение

Распространённое заблуждение: дополнение сигнала нулями перед FFT улучшает частотное разрешение. Это не так. Дополнение нулями интерполирует между существующими частотными бинами, создавая более гладкий на вид спектр. Но оно не может раскрыть частотные детали, которых не было в исходных отсчётах.

Настоящее частотное разрешение зависит только от времени наблюдения T = N / fs. Чтобы различить две частоты, разделённые delta_f, нужно не менее T = 1 / delta_f секунд данных. Никакое дополнение нулями не меняет этот фундаментальный предел.

fourier-synthesis

Реализуйте

Шаг 1: DFT с нуля

DFT со сложностью O(N^2) напрямую следует из определения.

import math

class Complex:
    ...

def dft(x):
    N = len(x)
    result = []
    for k in range(N):
        total = Complex(0, 0)
        for n in range(N):
            angle = -2 * math.pi * k * n / N
            w = Complex(math.cos(angle), math.sin(angle))
            xn = x[n] if isinstance(x[n], Complex) else Complex(x[n])
            total = total + xn * w
        result.append(total)
    return result

Шаг 2: Обратное DFT

Та же структура: положительный показатель степени и деление на N.

def idft(X):
    N = len(X)
    result = []
    for n in range(N):
        total = Complex(0, 0)
        for k in range(N):
            angle = 2 * math.pi * k * n / N
            w = Complex(math.cos(angle), math.sin(angle))
            total = total + X[k] * w
        result.append(Complex(total.real / N, total.imag / N))
    return result

Шаг 3: FFT (Кули—Тьюки)

Рекурсивный FFT требует длины, равной степени 2. Разделите на чётные и нечётные элементы, рекурсивно обработайте их и объедините с twiddle-факторами.

def fft(x):
    N = len(x)
    if N <= 1:
        return [x[0] if isinstance(x[0], Complex) else Complex(x[0])]
    if N % 2 != 0:
        return dft(x)

    even = fft([x[i] for i in range(0, N, 2)])
    odd = fft([x[i] for i in range(1, N, 2)])

    result = [Complex(0)] * N
    for k in range(N // 2):
        angle = -2 * math.pi * k / N
        twiddle = Complex(math.cos(angle), math.sin(angle))
        t = twiddle * odd[k]
        result[k] = even[k] + t
        result[k + N // 2] = even[k] - t
    return result

Шаг 4: Вспомогательные функции спектрального анализа

def power_spectrum(X):
    return [xk.real ** 2 + xk.imag ** 2 for xk in X]

def convolve_fft(x, h):
    N = len(x) + len(h) - 1
    padded_N = 1
    while padded_N < N:
        padded_N *= 2

    x_padded = x + [0.0] * (padded_N - len(x))
    h_padded = h + [0.0] * (padded_N - len(h))

    X = fft(x_padded)
    H = fft(h_padded)

    Y = [xk * hk for xk, hk in zip(X, H)]

    y = idft(Y)
    return [y[n].real for n in range(N)]

Используйте

Для практической работы используйте FFT из numpy: оно опирается на высокооптимизированные библиотеки C.

import numpy as np

signal = np.sin(2 * np.pi * 5 * np.arange(256) / 256)
spectrum = np.fft.fft(signal)
freqs = np.fft.fftfreq(256, d=1/256)

power = np.abs(spectrum) ** 2

positive_freqs = freqs[:len(freqs)//2]
positive_power = power[:len(power)//2]

Для оконных функций и более продвинутого спектрального анализа:

from scipy.signal import windows, stft

window = windows.hann(256)
windowed = signal * window
spectrum = np.fft.fft(windowed)

Для свёртки:

from scipy.signal import fftconvolve

result = fftconvolve(signal, kernel, mode='full')

Для спектрограмм:

from scipy.signal import stft

frequencies, times, Zxx = stft(signal, fs=sample_rate, nperseg=256)
spectrogram = np.abs(Zxx) ** 2

Матрица спектрограммы имеет форму (n_frequencies, n_time_frames). Каждый столбец — это спектр мощности в одном временном окне. Именно его аудио ML-модели используют как вход.

Доведите до результата

Запустите code/fourier.py, чтобы создать outputs/prompt-spectral-analyzer.md.

Упражнения

1. Идентификация чистого тона. Создайте сигнал с одной синусоидой неизвестной частоты (от 1 до 50 Hz), дискретизированный с частотой 128 Hz в течение 1 секунды. Используйте своё DFT, чтобы определить частоту. Убедитесь, что ответ совпадает. Теперь добавьте гауссов шум со стандартным отклонением 0.5 и повторите. Как шум влияет на спектр?

2. Проверка FFT против DFT. Сгенерируйте случайный сигнал длины 64. Вычислите и DFT (O(N^2)), и FFT. Проверьте, что все коэффициенты совпадают с точностью до 1e-10. Замерьте время обеих функций на сигналах длины 256, 512, 1024 и 2048. Постройте график отношения времени DFT ко времени FFT.

3. Доказательство теоремы о свёртке примером. Создайте сигнал x = [1, 2, 3, 4, 0, 0, 0, 0] и фильтр h = [1, 1, 1, 0, 0, 0, 0, 0]. Вычислите их циклическую свёртку напрямую (вложенным циклом). Затем вычислите её через FFT (преобразовать, перемножить, выполнить обратное преобразование). Убедитесь, что результаты совпадают. Теперь выполните линейную свёртку, правильно дополнив сигналы нулями.

4. Эффекты оконных функций. Создайте сигнал, равный сумме двух синусоид на 10 Hz и 12 Hz (очень близкие частоты). Дискретизируйте его с частотой 128 Hz в течение 1 секунды. Вычислите спектр мощности без окна, с окном Ханна и с окном Хэмминга. Какое окно позволяет легче всего различить два пика? Почему?

5. Анализ позиционной кодировки. Сгенерируйте синусоидальные позиционные кодировки для d_model = 128 и max_pos = 512. Для каждой пары позиций (p1, p2) вычислите скалярное произведение их кодировок. Покажите, что скалярное произведение зависит только от |p1 - p2|, а не от абсолютных позиций. Что происходит со скалярным произведением по мере увеличения расстояния?

Ключевые термины

Термин	Значение
DFT (дискретное преобразование Фурье)	Преобразует N отсчётов временной области в N коэффициентов частотной области. Каждый коэффициент — это корреляция с комплексной синусоидой на данной частоте
FFT (быстрое преобразование Фурье)	Алгоритм O(N log N) для вычисления DFT. Алгоритм Кули—Тьюки рекурсивно разделяет чётные и нечётные индексы
Обратное DFT	Восстанавливает сигнал временной области из частотных коэффициентов. Та же формула, что у DFT, но с перевёрнутым знаком показателя степени и масштабированием 1/N
Частотный бин	Каждый индекс k в выходе DFT представляет частоту k*fs/N Hz. «Бин» — это дискретная ячейка частоты
DC-компонента	X[0], коэффициент нулевой частоты. Пропорционален среднему сигнала
Частота Найквиста	fs/2, максимальная частота, представимая при частоте дискретизации fs. Частоты выше неё дают алиасинг
Спектр мощности	|X[k]|^2, квадрат модуля каждого частотного коэффициента. Показывает распределение энергии по частотам
Спектр фаз	angle(X[k]), сдвиг фазы каждой частотной компоненты. В анализе часто игнорируется
Спектральная утечка	Ложное частотное содержимое, вызванное рассмотрением непериодического сигнала как периодического. Уменьшается оконной функцией
Оконная функция	Сглаживающая функция (Ханна, Хэмминга, Блэкмана), применяемая перед DFT для уменьшения спектральной утечки
Twiddle-фактор	Комплексная экспонента e^(-2pii*k/N), используемая для объединения под-DFT в вычислении FFT «бабочкой»
Теорема о свёртке	Свёртка во временной области равна поточечному умножению в частотной области. Основа обработки сигналов и CNN
Циклическая свёртка	Свёртка, в которой сигнал заворачивается по кругу. Именно её естественным образом вычисляет DFT
Линейная свёртка	Стандартная свёртка без заворачивания. Достигается дополнением нулями перед DFT
Теорема Парсеваля	Суммарная энергия сохраняется при преобразовании Фурье. sum |x[n]|^2 = (1/N) sum |X[k]|^2
Алиасинг	Ситуация, когда частоты выше Найквиста выглядят как более низкие из-за недостаточной частоты дискретизации

Дополнительное чтение

• Cooley & Tukey: An Algorithm for the Machine Calculation of Complex Fourier Series (1965) — оригинальная статья об FFT, изменившая вычислительную технику
• 3Blue1Brown: But what is the Fourier Transform? — лучшее визуальное введение в преобразования Фурье
• Lee-Thorp et al.: FNet: Mixing Tokens with Fourier Transforms (2021) — замена self-attention на FFT в трансформерах
• Smith: The Scientist and Engineer’s Guide to Digital Signal Processing — бесплатный онлайн-учебник, подробно охватывающий FFT, оконные функции и спектральный анализ
• Vaswani et al.: Attention Is All You Need (2017) — синусоидальные позиционные кодировки, выведенные из частотной декомпозиции Фурье
• Radford et al.: Whisper (2022) — распознавание речи, использующее мел-спектрограммы как входное представление

---

Источник: The Fourier Transform Навигация: ← 01.19 — Комплексные числа для AI · ↑ Фаза 1 — Математические основы · Полный каталог · → 01.21 — Теория графов для машинного обучения