Охота на слабый импульс в Python: логарифмический спектр проявляет скрытый сигнал

В этой статье мы спрячем слабый прямоугольный импульс внутрь мощной синусоиды
и покажем, что на обычном линейном спектре он невидим, а стоит переключить ось
амплитуд в децибелы — и проявляется характерный sinc-узор боковых лепестков
импульса. Попутно сверим численный спектр с аналитическим и убедимся, что БПФ
не врёт.

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

Представим, что прибор записал сигнал, в котором «звучит» сильная несущая
частота, а где-то посередине записи проскочил короткий прямоугольный импульс
в 50 раз меньшей амплитуды:

где , Гц, , с,
а — центр записи. Исходные данные:

частота дискретизации Гц;
длительность записи с, отсчётов ;
разрешение по частоте Гц.

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

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

Зададим параметры и соберём несущую, импульс и их сумму.

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

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

f0 = 50.0                       # частота несущей, Гц
A  = 5.0                        # амплитуда несущей
eps   = 0.1                     # амплитуда импульса
T_imp = 0.2                     # длительность импульса, с
tc    = T / 2                   # центр импульса

f1 = A * np.sin(2 * np.pi * f0 * t)
f2 = eps * ((t >= tc - T_imp/2) & (t <= tc + T_imp/2)).astype(float)
sig = f1 + f2

Посмотрим фрагмент вокруг центра записи, где спрятан импульс.

Рисунок 1. Слева — несущая (синяя) и импульс (красный) порознь; справа — суммарный сигнал, в котором импульс практически не виден

На правом графике импульс амплитуды 0.1 тонет в синусоиде амплитуды 5: глаз
видит только несущую. Во временной области такая «стеганография» надёжна.
Посмотрим, справится ли спектр.

Шаг 2. Односторонний амплитудный спектр

Вычислим нормированный односторонний спектр — тот же приём с нормировкой
, что и в предыдущей статье:

def one_sided_spectrum(x):
    Y = np.fft.fft(x) / len(x)
    n = len(Y)
    freq = np.arange(n // 2) * Fs / n
    amp = np.abs(Y[:n // 2])
    amp[1:] *= 2
    return freq, amp

freq, amp_sig = one_sided_spectrum(sig)

Нарисуем спектр суммарного сигнала в линейном масштабе, ограничив ось до
150 Гц.

Шаг 3. Линейный масштаб: один пик правит балом

Рисунок 2. Слева — линейный спектр: одинокий пик на 50 Гц, вклад импульса не виден; справа — тот же спектр в дБ: проявились боковые лепестки sinc-функции

Левая панель показывает то, чего мы ожидали: острый пик на 50 Гц высотой 5.0,
а всё остальное прижато к нулю. Прямоугольный импульс «размазан» по многим
гармоникам, и каждая из них настолько мала по сравнению с несущей, что на
линейной шкале сливается с фоном.

Шаг 4. Логарифмический масштаб: охота удалась

Переведём амплитуды в децибелы:

где — страховка от логарифма нуля.

db = 20 * np.log10(amp_sig + 1e-12)

Правая панель рисунка 2 показывает результат: над шумовой полкой проступает
характерный «забор» боковых лепестков — это спектр прямоугольного импульса,
который в линейном масштабе был невидим. Нули sinc-функции приходятся на
частоты, кратные Гц, и между ними видны лепестки,
убывающие примерно как . Охота удалась: слабый сигнал обнаружен.

Шаг 5. Проверка: численный спектр против аналитического

Убедимся, что БПФ не врёт. Для синусоиды аналитический спектр — это
единственный пик амплитуды на частоте . Для прямоугольного
импульса непрерывное преобразование Фурье даёт

где , а деление на появляется
из-за нормировки одностороннего спектра.

# Аналитический спектр импульса
freq_th = freq[freq <= 150]
amp_theor = (eps * T_imp / T) * np.abs(np.sinc(freq_th * T_imp))

Рисунок 3. Слева — спектр синусоиды: красная точка аналитики точно совпадает с пиком БПФ; справа — спектр прямоугольного импульса: пунктирная аналитическая кривая sinc ложится поверх численного спектра

Обе панели подтверждают: численный спектр совпадает с теорией. Синусоида даёт
ровно один пик нужной высоты, а импульс воспроизводит форму sinc с правильными
нулями и лепестками.

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

Увеличьте амплитуду импульса до : на каком значении он
становится заметен уже на линейном спектре?
Сократите длительность импульса до с: нули sinc
разойдутся шире. Сколько лепестков помещается в полосе до 150 Гц?
Добавьте белый шум sig += 0.05 * np.random.normal(size=N) и посмотрите,
при каком уровне шума лепестки ещё различимы на логарифмическом спектре.

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

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

# Охота на слабый импульс: спектр в логарифмическом масштабе (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"
GREEN = "#2a8a4b"


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.0
T  = 4.0
N  = int(Fs * T)
t  = np.arange(N) / Fs

f0 = 50.0
A  = 5.0
eps   = 0.1
T_imp = 0.2
tc    = T / 2.0

# ---------- Сигналы ----------
f1  = A * np.sin(2 * np.pi * f0 * t)
f2  = eps * ((t >= tc - T_imp/2) & (t <= tc + T_imp/2)).astype(float)
sig = f1 + f2


def one_sided_spectrum(x):
    Y = np.fft.fft(x) / len(x)
    n = len(Y)
    freq = np.arange(n // 2) * Fs / n
    amp = np.abs(Y[:n // 2])
    amp[1:] *= 2
    return freq, amp


freq, amp_f1  = one_sided_spectrum(f1)
_,    amp_f2  = one_sided_spectrum(f2)
_,    amp_sig = one_sided_spectrum(sig)

idx  = freq <= 150
zoom = (t >= tc - 0.5) & (t <= tc + 0.5)

# ---------- Рисунок 1: временная область ----------
fig, axes = plt.subplots(1, 2)

axes[0].plot(t[zoom], f1[zoom], color=BLUE, linewidth=1.0, label="Несущая 50 Гц")
axes[0].plot(t[zoom], f2[zoom], color=RED,  linewidth=1.8, label="Импульс (x50)")
axes[0].set_xlabel("Время, с")
axes[0].set_ylabel("Амплитуда")
axes[0].legend(frameon=False, fontsize=9)
finish(axes[0], "Несущая и слабый импульс (фрагмент)")

axes[1].plot(t[zoom], sig[zoom], color=GREEN, linewidth=1.0)
axes[1].set_xlabel("Время, с")
axes[1].set_ylabel("Амплитуда")
finish(axes[1], "Суммарный сигнал: импульс не виден")

fig.tight_layout()
fig.savefig("weak_pulse_time.png")
plt.close(fig)

# ---------- Рисунок 2: линейный vs логарифмический спектр ----------
fig, axes = plt.subplots(1, 2)

axes[0].plot(freq[idx], amp_sig[idx], color=GREEN, linewidth=1.4)
axes[0].set_xlabel("Частота, Гц")
axes[0].set_ylabel("Амплитуда гармоник")
finish(axes[0], "Линейный масштаб: один пик доминирует")

db = 20 * np.log10(amp_sig[idx] + 1e-12)
axes[1].plot(freq[idx], db, color=GREEN, linewidth=1.4)
axes[1].set_xlabel("Частота, Гц")
axes[1].set_ylabel("Амплитуда, дБ")
finish(axes[1], "Логарифмический масштаб: проявился sinc-узор импульса")

fig.tight_layout()
fig.savefig("weak_pulse_spectrum.png")
plt.close(fig)

# ---------- Рисунок 3: численный vs аналитический ----------
fig, axes = plt.subplots(1, 2)

axes[0].plot(freq[idx], amp_f1[idx], color=BLUE, linewidth=1.4, label="БПФ")
axes[0].plot(f0, A, "o", color=RED, markersize=8, markeredgecolor="white",
             label=f"Аналитика: A={A}")
axes[0].axvline(f0, color=RED, linestyle="--", alpha=0.4)
axes[0].set_xlim(0, 100)
axes[0].set_xlabel("Частота, Гц")
axes[0].set_ylabel("Амплитуда гармоник")
axes[0].legend(frameon=False, fontsize=9)
finish(axes[0], "Синусоида: БПФ совпадает с аналитикой")

freq_th   = freq[idx]
amp_theor = (eps * T_imp / T) * np.abs(np.sinc(freq_th * T_imp))
axes[1].plot(freq[idx], amp_f2[idx], color=RED, linewidth=1.4, label="БПФ")
axes[1].plot(freq_th, amp_theor, "k--", linewidth=1.2, label="Аналитический sinc")
axes[1].set_xlim(0, 100)
axes[1].set_xlabel("Частота, Гц")
axes[1].set_ylabel("Амплитуда гармоник")
axes[1].legend(frameon=False, fontsize=9)
finish(axes[1], "Прямоугольный импульс: БПФ vs непрерывный спектр")

fig.tight_layout()
fig.savefig("weak_pulse_analytical.png")
plt.close(fig)

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