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

Главный инструмент электрофизиолога — спектр: сколько в сигнале какой частоты. Посчитаем его на чистом Python.

Спектр мощности показывает, как энергия сигнала распределена по частотам; мощность в полосе — сумма по её частотам.

Идея преобразования Фурье

Любой сигнал можно представить как сумму синусоид. Чтобы узнать, сколько в сигнале частоты $f$, его умножают на синус и косинус этой частоты и усредняют:

$$ a_f = \frac{1}{N}\sum_{n} x_n \cos(2\pi f t_n), \quad b_f = \frac{1}{N}\sum_{n} x_n \sin(2\pi f t_n) $$

Мощность на частоте $f$ равна $P_f = a_f^2 + b_f^2$. Если в сигнале есть эта частота, синусоиды «срезонируют» и сумма будет большой; если нет — близкой к нулю.

Мощность по полосам

Чтобы оценить, например, альфа-ритм, складывают $P_f$ по всем частотам полосы $8\le f \le 12$ Гц. Отношение мощностей полос (скажем, $\alpha/\beta$) — частый признак в BCI и нейробиоуправлении.

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

Соберём сигнал из чистого альфа (10 Гц, сильный) и бета (20 Гц, слабее) плюс шум, затем посчитаем мощность в полосах:

import math, random
random.seed(7)
fs, N = 128, 256
t = [i/fs for i in range(N)]
sig = [2.0*math.sin(2*math.pi*10*x) + 1.0*math.sin(2*math.pi*20*x)
       + 0.5*random.gauss(0,1) for x in t]

def power_at(f):
    re = sum(sig[i]*math.cos(2*math.pi*f*t[i]) for i in range(N)) / N
    im = sum(sig[i]*math.sin(2*math.pi*f*t[i]) for i in range(N)) / N
    return re*re + im*im

def band_power(lo, hi):
    return sum(power_at(f) for f in range(lo, hi+1))

alpha = band_power(8, 12)
beta  = band_power(13, 30)
print("Мощность alpha (8-12 Гц):", round(alpha, 3))
print("Мощность beta  (13-30 Гц):", round(beta, 3))
print("Отношение alpha/beta:", round(alpha/beta, 2))

Вывод:

Мощность alpha (8-12 Гц): 0.964
Мощность beta  (13-30 Гц): 0.269
Отношение alpha/beta: 3.58

Мы заложили сильную альфу и слабую бету — и спектр это подтвердил: мощность альфы примерно в 3.6 раза больше. Формула спектра сошлась с тем, что мы синтезировали.

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

• Считать, что один отсчёт несёт частоту. Частота — свойство последовательности отсчётов.
• Анализировать слишком короткое окно: на 0.1 с нельзя различить близкие частоты (разрешение $\approx 1/T$).
• Забывать про окно (Хэннинга и т.п.) на реальных данных — иначе края дают ложные частоты (растекание спектра).

Итог

• Спектр: умножаем сигнал на синус/косинус частоты и усредняем; $P_f = a_f^2 + b_f^2$.
• Мощность полосы — сумма $P_f$ по её частотам.
• Отношения мощностей полос (α/β) — рабочий признак для BCI.