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

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

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

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

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

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

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

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

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

clear; clc;

fs = 1000;                 // частота дискретизации, Гц
T  = 1.0;                  // длительность сигнала, с
N  = fs * T;               // число отсчётов
t  = (0:N-1) / fs;         // ось времени, с

f1 = 5;    A1 = 1.0;       // первая гармоника: частота, Гц и амплитуда
f2 = 20;   A2 = 0.5;       // вторая гармоника
f3 = 100;  A3 = 0.2;       // третья гармоника

s = A1*sin(2*%pi*f1*t) + A2*sin(2*%pi*f2*t) + A3*sin(2*%pi*f3*t);

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

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

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

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

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

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

S = fft(s);                       // комплексный спектр (двусторонний)
f = (0:N-1) * fs / N;             // ось частот, Гц

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

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

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

amp_raw = abs(S);

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

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

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

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

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

nh   = N/2;
fpos = f(1:nh);
amp  = 2 * abs(S(1:nh)) / N;
amp(1) = amp(1) / 2;              // постоянную составляющую не удваиваем

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

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

idx = find(fpos <= 150);
k1 = find(abs(fpos - f1) < 0.5);
k2 = find(abs(fpos - f2) < 0.5);
k3 = find(abs(fpos - f3) < 0.5);

printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f1, amp(k1), A1);
printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f2, amp(k2), A2);
printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f3, amp(k3), 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*rand(1, N, "normal") и посмотрите, при каком
уровне шума пик 0.2 ещё различим.
Задайте Гц — частоту «между» бинами сетки — и убедитесь, что
пик «размазывается» по соседним частотам. Это и есть эффект утечки спектра.

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

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

// Амплитудный спектр сигнала в Scilab: разложение суммы синусоид по частотам
clear; clc;

// ---------- Параметры ----------
fs = 1000;                 // частота дискретизации, Гц
T  = 1.0;                  // длительность сигнала, с
N  = fs * T;               // число отсчётов
t  = (0:N-1) / fs;         // ось времени, с

f1 = 5;    A1 = 1.0;       // первая гармоника: частота, Гц и амплитуда
f2 = 20;   A2 = 0.5;       // вторая гармоника
f3 = 100;  A3 = 0.2;       // третья гармоника

// ---------- Сигнал ----------
s = A1*sin(2*%pi*f1*t) + A2*sin(2*%pi*f2*t) + A3*sin(2*%pi*f3*t);

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

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

// Рисунок 1: сигнал во временной области
fh1 = scf(1); clf(fh1);
fh1.figure_size = [900, 450];
plot(t, s, "b-", "thickness", 1.5);
title("Сигнал во временной области: сумма трёх синусоид");
xlabel("Время, с");
ylabel("Амплитуда");
gca().box = "on";
xgrid();
xs2png(fh1.figure_id, "fft_signal.png");

// Рисунок 2: «сырой» спектр без нормировки (промежуточный, несовершенный вид)
fh2 = scf(2); clf(fh2);
fh2.figure_size = [900, 450];
plot(f(1:N/2), amp_raw(1:N/2), "b-", "thickness", 2);
title("Модуль спектра без нормировки: амплитуды зависят от N");
xlabel("Частота, Гц");
ylabel("|S(k)|");
gca().box = "on";
xgrid();
xs2png(fh2.figure_id, "fft_spectrum_raw.png");

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

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

// Рисунок 3: нормированный амплитудный спектр с отмеченными пиками
fh3 = scf(3); clf(fh3);
fh3.figure_size = [900, 450];
h_sp = plot(fpos(idx), amp(idx), "b-", "thickness", 2);
// Отметим три ожидаемые гармоники (разрешение по частоте fs/N = 1 Гц)
k1 = find(abs(fpos - f1) < 0.5);
k2 = find(abs(fpos - f2) < 0.5);
k3 = find(abs(fpos - f3) < 0.5);
h_pk = plot([fpos(k1), fpos(k2), fpos(k3)], ..
            [amp(k1),  amp(k2),  amp(k3)], "r.", "markersize", 12);
legend([h_sp, h_pk], ["Амплитудный спектр"; "Гармоники 5, 20 и 100 Гц"], 1);
title("Амплитудный спектр сигнала (нормировка 2/N)");
xlabel("Частота, Гц");
ylabel("Амплитуда гармоник");
gca().box = "on";
xgrid();
xs2png(fh3.figure_id, "fft_spectrum.png");

// ---------- Проверка числами ----------
printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f1, amp(k1), A1);
printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f2, amp(k2), A2);
printf("Пик на %d Гц: амплитуда %.3f (ожидалось %.1f)\n", f3, amp(k3), A3);

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