Дискретное преобразование Фурье (ДПФ)

Жемчужина курса: записываем точную формулу спектра и считаем его настоящим кодом на cmath.

ДПФ (дискретное преобразование Фурье) переводит N отсчётов сигнала в N комплексных коэффициентов: X[k] = sum(x[n]*exp(-2j*pi*k*n/N)). Каждый X[k] — амплитуда и фаза частоты с номером k.

Мы подошли к главному инструменту анализа. ДПФ — это та самая операция, что превращает осциллограмму в спектр. И, к счастью, в стандартной библиотеке Python есть модуль cmath с комплексной экспонентой — значит, ДПФ можно реализовать честно, без numpy, и посчитать спектр реального сигнала прямо здесь.

Формула и её смысл

Комплексная экспонента exp(-2j*pi*k*n/N) — это «вращающийся пробный синус» частоты k: её вещественная часть — косинус, мнимая — синус (с минусом). Умножая сигнал на неё и суммируя, мы за один присест измеряем и косинусную, и синусную составляющую частоты k. Результат X[k] — комплексное число: его модуль |X[k]| — амплитуда, аргумент — фаза.

import math, cmath

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

# Прямоугольный импульс: 4 единицы, затем 4 нуля
x = [1, 1, 1, 1, 0, 0, 0, 0]
X = dft(x)
for k, c in enumerate(X):
    re = round(c.real, 3) or 0.0
    im = round(c.imag, 3) or 0.0
    print(f"X[{k}] = {re:+.3f} {im:+.3f}j   |X|={abs(c):.3f}")

Вывод:

X[0] = +4.000 +0.000j   |X|=4.000
X[1] = +1.000 -2.414j   |X|=2.613
X[2] = +0.000 +0.000j   |X|=0.000
X[3] = +1.000 -0.414j   |X|=1.082
X[4] = +0.000 +0.000j   |X|=0.000
X[5] = +1.000 +0.414j   |X|=1.082
X[6] = +0.000 +0.000j   |X|=0.000
X[7] = +1.000 +2.414j   |X|=2.613

X[0] = 4 — это сумма всех отсчётов (постоянная составляющая, «уровень»). Остальные коэффициенты описывают, какие частоты нужны, чтобы собрать прямоугольник. Спектр посчитан настоящим кодом.

Спектр реального сигнала-смеси

Теперь главное: подадим в ДПФ смесь двух тонов и увидим, как алгоритм точно вытащит их частоты и амплитуды.

import math, cmath

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)]

N, fs = 16, 16.0
# тон 2 Гц (амплитуда 1.0) + тон 5 Гц (амплитуда 0.5)
sig = [math.sin(2 * math.pi * 2 * n / fs) + 0.5 * math.sin(2 * math.pi * 5 * n / fs)
       for n in range(N)]
X = dft(sig)

print("Гц | амплитуда")
for k in range(N // 2 + 1):
    amp = round(abs(X[k]) / (N / 2), 3) or 0.0
    bar = "#" * int(round(amp * 20))
    print(f"{k * fs / N:4.1f} | {bar} {amp}")

Вывод:

Гц | амплитуда
 0.0 |  0.0
 1.0 |  0.0
 2.0 | #################### 1.0
 3.0 |  0.0
 4.0 |  0.0
 5.0 | ########## 0.5
 6.0 |  0.0
 7.0 |  0.0
 8.0 |  0.0

Две линии: 2 Гц амплитудой 1.0 и 5 Гц амплитудой 0.5. ДПФ безошибочно восстановил рецепт сигнала из голых отсчётов. Это и есть «рентген», который мы обещали.

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

Каждый коэффициент X[k] соответствует частоте f = k*fs/N. Шаг по частоте между соседними коэффициентами — fs/N; это разрешение по частоте. Чтобы различить два близких тона, нужно достаточно много отсчётов N (длиннее запись — тоньше «гребёнка» частот). Прямое суммирование требует N операций на каждый из N коэффициентов — всего O(N²). Для 1000 отсчётов это миллион умножений, для миллиона отсчётов — триллион: слишком медленно. Поэтому на практике применяют быстрый алгоритм БПФ, которому посвящён следующий раздел. Но математически БПФ считает ровно то же ДПФ — просто умнее.

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

Забыть нормировку. Сырой |X[k]| зависит от N; для физической амплитуды тона делят на N/2 (а постоянную составляющую — на N).
Считать частоту коэффициента неправильно. Частота бина k равна k*fs/N, а не просто k.
Применять ДПФ «как есть» к длинным сигналам. Прямое ДПФ — O(N²); для реальных данных нужен БПФ.

Итог

ДПФ: X[k] = sum(x[n]*exp(-2j*pi*k*n/N)) — переводит отсчёты в спектр.
Реализуется на чистом Python через cmath.exp без сторонних библиотек.
|X[k]| — амплитуда, аргумент — фаза; частота бина — k*fs/N.
Прямое ДПФ медленное (O(N²)); на практике используют БПФ.