Амплитудный спектр сигнала в Python: БПФ раскладывает клубок синусоид на три пика

В этой статье мы построим амплитудный спектр сигнала из трёх синусоид на Python:
вычислим быстрое преобразование Фурье (БПФ) через NumPy, поймаем ошибку
нормировки на «сыром» спектре и исправим её так, чтобы пики спектра показали
истинные амплитуды гармоник — 1.0, 0.5 и 0.2.

Что будем анализировать

Пусть прибор снимает сигнал, в котором одновременно «звучат» три гармоники:

где , и — частоты в герцах, а коэффициенты перед синусами — их
амплитуды. Исходные данные:

частота дискретизации Гц — тысяча отсчётов в секунду;
длительность записи с, значит, отсчётов ;
разрешение по частоте Гц — этого хватит, чтобы
различить наши гармоники.

Задача: по одному лишь записанному сигналу узнать, из каких частот он
состоит и какова амплитуда каждой.

Шаг 1. Генерируем сигнал

Зададим параметры и соберём сигнал как сумму трёх синусоид.

import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

plt.rcParams.update({"figure.figsize": (9, 4.5), "axes.grid": True})

fs = 1000                 # частота дискретизации, Гц
T = 1.0                   # длительность сигнала, с
N = int(fs * T)           # число отсчётов
t = np.arange(N) / fs     # ось времени, с

f1, f2, f3 = 5, 20, 100   # частоты гармоник, Гц
A1, A2, A3 = 1.0, 0.5, 0.2

s = A1 * np.sin(2 * np.pi * f1 * t) + \
    A2 * np.sin(2 * np.pi * f2 * t) + \
    A3 * np.sin(2 * np.pi * f3 * t)

Построим сигнал во временной области и сохраним рисунок в PNG.

Рисунок 1. Сигнал во временной области: сумма трёх синусоид

По временной картине состав сигнала не угадать: кривая выглядит как клубок
проводов, в котором перемешаны все три частоты. Глаз видит только период
около 0.2 с — это самая сильная гармоника 5 Гц, — а слабая сотня герц прячется
в мелкой ряби. Именно для таких случаев и существует спектр.

Шаг 2. БПФ и ось частот

Преобразование Фурье раскладывает сигнал по гармоникам. В дискретном виде оно
вычисляется формулой

где — отсчёты сигнала, а — комплексные коэффициенты спектра.
Каждому индексу соответствует своя частота . В NumPy
всё это делает одна функция:

S = np.fft.fft(s)                 # комплексный спектр (двусторонний)
f = np.arange(N) * fs / N         # ось частот, Гц

Спектр S — двусторонний: пики стоят и на положительных, и на отрицательных
частотах, а для вещественного сигнала две половины зеркально симметричны.

Шаг 3. Первая попытка: спектр без нормировки

Самый простой способ посмотреть спектр — нарисовать модуль :

amp_raw = np.abs(S)

Рисунок 2. Модуль спектра без нормировки: амплитуды зависят от N

Пики стоят там, где надо, — на 5, 20 и 100 Гц, — но высоты у них странные:
500, 250 и 100 вместо ожидаемых 1.0, 0.5 и 0.2. Дело в том, что БПФ «суммирует»
отсчётов, и высота пика пропорциональна для гармоники единичной
амплитуды. Такой спектр годится для поиска частот, но не для измерения
амплитуд: удвоим длину записи — и «амплитуды» на графике тоже удвоятся.
Исправим нормировку.

Шаг 4. Нормировка 2/N: физический амплитудный спектр

Чтобы пики показывали истинные амплитуды, поделим модуль спектра на и
оставим только положительные частоты, удвоив их, — ведь половина энергии
гармоники ушла в зеркальное отрицательное плечо:

где — индекс пика на частоте -й гармоники. Постоянную составляющую
() удваивать не нужно: у неё нет «зеркала».

nh = N // 2
fpos = f[:nh]
amp = 2 * np.abs(S[:nh]) / N
amp[0] /= 2                       # постоянную составляющую не удваиваем

Шаг 5. Финальный спектр и проверка числами

Ограничим ось до 150 Гц, отметим найденные пики и сверим их высоты с
заданными амплитудами.

idx = fpos <= 150
k1 = np.argmin(np.abs(fpos - f1))
k2 = np.argmin(np.abs(fpos - f2))
k3 = np.argmin(np.abs(fpos - f3))

print(f"Пик на {f1} Гц: амплитуда {amp[k1]:.3f} (ожидалось {A1})")
print(f"Пик на {f2} Гц: амплитуда {amp[k2]:.3f} (ожидалось {A2})")
print(f"Пик на {f3} Гц: амплитуда {amp[k3]:.3f} (ожидалось {A3})")

Консоль отвечает:

Пик на 5 Гц: амплитуда 1.000 (ожидалось 1.0)
Пик на 20 Гц: амплитуда 0.500 (ожидалось 0.5)
Пик на 100 Гц: амплитуда 0.200 (ожидалось 0.2)

Рисунок 3. Амплитудный спектр сигнала (нормировка 2/N)

Клубок распутан: три острых пика ровно на заданных частотах, а их высоты
совпадают с амплитудами гармоник до третьего знака. Спектр сделал то, что не
смогла временная область, — разделил составляющие и измерил каждую.

Что попробовать самостоятельно

Уменьшите до 0.5 с: разрешение станет Гц. Проверьте, что
произойдёт с пиком на 5 Гц, если задать частоту 6 Гц.
Добавьте шум s = s + 0.1 * np.random.normal(size=N) и посмотрите, при
каком уровне шума пик 0.2 ещё различим.
Задайте Гц — частоту «между» бинами сетки — и убедитесь, что
пик «размазывается» по соседним частотам. Это и есть эффект утечки спектра.

Полный код скрипта

Весь листинг fft_base.py целиком — от генерации сигнала до сохранения
графиков:

# Амплитудный спектр сигнала в Python: разложение суммы синусоид по частотам
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

# ---------- Оформление ----------
plt.rcParams.update({
    "figure.figsize": (9, 4.2),
    "figure.dpi": 160,
    "axes.edgecolor": "#c9ccc4",
    "axes.linewidth": 1.0,
    "axes.grid": True,
    "axes.axisbelow": True,
    "grid.color": "#e7e9e3",
    "axes.titlesize": 13,
    "axes.titleweight": "bold",
    "axes.labelsize": 11,
    "axes.labelcolor": "#4f5a53",
    "xtick.labelsize": 10,
    "ytick.labelsize": 10,
    "xtick.color": "#4f5a53",
    "ytick.color": "#4f5a53",
})
BLUE = "#1e5bc4"                 # линия графика
RED = "#d93a2b"                  # акцентные маркеры

def finish(ax, title):
    # убираем лишние рамки и выравниваем заголовок по левому краю
    ax.set_title(title, loc="left", pad=12)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)

# ---------- Параметры ----------
fs = 1000                 # частота дискретизации, Гц
T = 1.0                   # длительность сигнала, с
N = int(fs * T)           # число отсчётов
t = np.arange(N) / fs     # ось времени, с

f1, f2, f3 = 5, 20, 100   # частоты гармоник, Гц
A1, A2, A3 = 1.0, 0.5, 0.2

# ---------- Сигнал ----------
s = A1 * np.sin(2 * np.pi * f1 * t) + \
    A2 * np.sin(2 * np.pi * f2 * t) + \
    A3 * np.sin(2 * np.pi * f3 * t)

# ---------- БПФ ----------
S = np.fft.fft(s)                    # комплексный спектр (двусторонний)
f = np.arange(N) * fs / N            # ось частот, Гц

# Первая попытка: модуль спектра без нормировки
amp_raw = np.abs(S)

# Рисунок 1: сигнал во временной области
fig, ax = plt.subplots()
ax.plot(t, s, color=BLUE, linewidth=1.1)
ax.set_xlabel("Время, с")
ax.set_ylabel("Амплитуда")
finish(ax, "Сигнал во временной области: сумма трёх синусоид")
fig.tight_layout()
fig.savefig("fft_signal.png")
plt.close(fig)

# Рисунок 2: «сырой» спектр без нормировки (промежуточный, несовершенный вид)
fig, ax = plt.subplots()
ax.plot(f[:N // 2], amp_raw[:N // 2], color=BLUE, linewidth=1.8)
ax.set_xlabel("Частота, Гц")
ax.set_ylabel("|S(k)|")
finish(ax, "Модуль спектра без нормировки: амплитуды зависят от N")
fig.tight_layout()
fig.savefig("fft_spectrum_raw.png")
plt.close(fig)

# ---------- Нормировка: физический амплитудный спектр ----------
# Односторонний спектр: положительные частоты, амплитуды удвоены,
# чтобы учесть симметричный вклад отрицательных частот
nh = N // 2
fpos = f[:nh]
amp = 2 * np.abs(S[:nh]) / N
amp[0] /= 2                          # постоянную составляющую не удваиваем

# Ограничим частотную ось до 150 Гц — выше значимых компонент нет
idx = fpos <= 150

# Рисунок 3: нормированный амплитудный спектр с отмеченными пиками
fig, ax = plt.subplots()
ax.plot(fpos[idx], amp[idx], color=BLUE, linewidth=1.8, label="Амплитудный спектр")
k1 = np.argmin(np.abs(fpos - f1))
k2 = np.argmin(np.abs(fpos - f2))
k3 = np.argmin(np.abs(fpos - f3))
ax.plot([fpos[k1], fpos[k2], fpos[k3]], [amp[k1], amp[k2], amp[k3]],
        "o", color=RED, markersize=7, markeredgecolor="white",
        label="Гармоники 5, 20 и 100 Гц")
ax.set_ylim(top=1.15)
ax.legend(frameon=False, fontsize=10)
ax.set_xlabel("Частота, Гц")
ax.set_ylabel("Амплитуда гармоник")
finish(ax, "Амплитудный спектр сигнала (нормировка 2/N)")
fig.tight_layout()
fig.savefig("fft_spectrum.png")
plt.close(fig)

# ---------- Проверка числами ----------
print(f"Пик на {f1} Гц: амплитуда {amp[k1]:.3f} (ожидалось {A1})")
print(f"Пик на {f2} Гц: амплитуда {amp[k2]:.3f} (ожидалось {A2})")
print(f"Пик на {f3} Гц: амплитуда {amp[k3]:.3f} (ожидалось {A3})")

Итог: мы прошли путь от «клубка» во временной области до спектра с тремя
пиками точных амплитуд и убедились, что нормировка делает спектр
физическим. Та же задача, решённая в Scilab, разобрана
в парной статье — полезно сравнить
оба инструмента.