Дискретное преобразование Фурье на пальцах

Самое сердце цифровой обработки сигналов — и его можно посчитать в несколько строк на чистом Python.

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

Формула ДПФ

Пусть есть $N$ отсчётов сигнала $x_0,x_1,\dots,x_{N-1}$. ДПФ вычисляет $N$ комплексных коэффициентов $X_k$ по формуле:

$$X_k=\sum_{n=0}^{N-1} x_n\, e^{-i\,2\pi kn/N},\qquad k=0,1,\dots,N-1.$$

Каждый коэффициент $X_k$ отвечает за частоту $k$. Его модуль $|X_k|$ — амплитуда этой частоты в сигнале, аргумент — фаза. Множители $e^{-i2\pi kn/N}$ — это степени корня из единицы, с которыми мы познакомились раньше.

Реализация в несколько строк

ДПФ — это просто двойная сумма, его легко написать без всяких библиотек:

import cmath, math
def dft(x):
    N = len(x)
    X = []
    for k in range(N):
        s = 0
        for n in range(N):
            s += x[n] * cmath.exp(-2j * math.pi * k * n / N)
        X.append(s)
    return X

# Сигнал: чистая синусоида одного периода на 4 отсчёта
signal = [0, 1, 0, -1]
spectr = dft(signal)
for k, X in enumerate(spectr):
    print(f"X[{k}] = {complex(round(X.real,3), round(X.imag,3))}, |X|={round(abs(X),3)}")

Вывод:

X[0] = 0j, |X|=0.0
X[1] = -2j, |X|=2.0
X[2] = 0j, |X|=0.0
X[3] = (-0+2j), |X|=2.0

Спектр «зажёгся» только на частотах $k=1$ и $k=3$ — именно там сосредоточена энергия синусоиды. Нулевой коэффициент $X_0$ равен нулю, потому что среднее значение сигнала равно нулю.

Проверка на импульсе

Если подать «импульс» $[1,0,0,0]$, все частоты должны присутствовать одинаково — резкий всплеск содержит весь спектр:

import cmath, math
def dft(x):
    N = len(x)
    return [sum(x[n]*cmath.exp(-2j*math.pi*k*n/N) for n in range(N)) for k in range(N)]
impulse = [1, 0, 0, 0]
print([round(abs(X), 3) for X in dft(impulse)])

Вывод:

[1.0, 1.0, 1.0, 1.0]

Все амплитуды равны $1$: импульс содержит все частоты поровну — фундаментальный факт обработки сигналов.

Как работает под капотом

Интуиция формулы такая: $e^{-i2\pi kn/N}$ — это «эталонное» вращение с частотой $k$. Умножая сигнал на этот эталон и суммируя, мы измеряем, насколько сигнал «похож» на данную частоту. Если в сигнале есть составляющая частоты $k$, сумма накапливается и даёт большой $|X_k|$; если нет — слагаемые крутятся в разные стороны и взаимно гасятся. Наивное ДПФ требует $N^2$ операций; знаменитый алгоритм БПФ использует симметрию корней из единицы и снижает это до $N\log N$.

Частые ошибки

Забывать знак минус в экспоненте: прямое ДПФ использует $e^{-i\dots}$, обратное — $e^{+i\dots}$.
Ожидать только один пик для вещественной синусоиды. У вещественного сигнала спектр симметричен — пики появляются парами ($k$ и $N-k$).
Интерпретировать высокие $k$ как «очень высокие частоты». Из-за периодичности $k$ и $N-k$ — это одна и та же физическая частота.

Итог

ДПФ переводит сигнал из времени в частоты: $X_k=\sum_n x_n e^{-i2\pi kn/N}$.
Модуль $|X_k|$ — амплитуда частоты $k$, аргумент — фаза.
Алгоритм основан на корнях из единицы; БПФ ускоряет его до $N\log N$.