Информатика / Цифровая обработка сигналов

Периодограмма спектральной плотности мощности

Периодограмма переводит квадраты модулей спектра в плотность мощности по частоте. Масштаб задают частота дискретизации и сумма квадратов отсчётов окна.

Опубликовано: Обновлено:

Формула

$$P_2[k]=\frac{|X_w[k]|^2}{f_s\sum_{n=0}^{N-1}|w[n]|^2},\quad X_w[k]=\sum_{n=0}^{N-1}w[n]x[n]e^{-2\pi i kn/N}$$
Схема Односторонняя периодограмма примера
Плотность, В²/Гц0120 Гц2 Гц4 ГцПлотность

Постоянная составляющая остаётся на нуле, чередование отсчётов попадает на частоту Найквиста 4 Гц. Среднее из блока не удалялось.

Обозначения

$P₂[k]$
двусторонняя оценка плотности в k-й ячейке, квадрат ед. сигнала/Гц
$Xw[k]$
коэффициент ДПФ сигнала после окна, ед. сигнала
$fₛ$
частота дискретизации, Гц
$w[n]$
вес окна, —
$x[n]$
отсчёт до умножения на окно, ед. сигнала
$N$
длина блока и ДПФ, отсчёты
n, k
индексы отсчёта и частоты, —

Условия применения

Здесь ДПФ имеет длину N и прямую нормировку без множителя. Отсчёты равномерны, энергия окна положительна. Если удаляется среднее или тренд, x[n] в формуле обозначает отсчёты после такой обработки.

Ограничения

Один блок даёт статистически изменчивую оценку плотности. Утечка из-за конечной длины и выбор окна влияют на форму пиков. Сумма оценок по частоте для произвольного окна отражает нормированный взвешенный квадрат сигнала.

Подробное объяснение

Сначала сигнал умножают на окно, затем вычисляют ДПФ. Квадрат модуля Xw измеряется в квадрате единиц исходного сигнала. Деление на fₛ и энергию окна добавляет в знаменатель герцы и приводит результат к плотности мощности.

Для вещественных отсчётов можно представить все данные на неотрицательных частотах. Ячейки парных частот складывают. Поэтому внутренние ячейки удваиваются, а нулевую ячейку и частоту Найквиста при чётном N оставляют без удвоения.

Частотный шаг равен fₛ/N. Умножение высоты ячейки на этот шаг даёт её вклад в средний квадрат. Из равенства Парсеваля следует, что сумма по всему спектру равна сумме квадратов взвешенных отсчётов, делённой на сумму квадратов окна.

При прямоугольном окне сумма квадратов окна равна N, а интегральная оценка совпадает с обычным средним квадратом блока. При другом окне отсчёты получают разные веса. Чтобы оценить разброс относительно среднего, сначала удаляют постоянную составляющую, если это уместно для задачи.

Как пользоваться формулой

  1. Запишите fₛ, длину блока и принятое окно; решите, нужно ли сохранить среднее значение записи.
  2. Умножьте отсчёты на окно, найдите комплексные коэффициенты ДПФ и квадраты их модулей.
  3. Разделите каждый квадрат на произведение частоты дискретизации и суммы квадратов весов окна.
  4. Для одностороннего спектра вещественной записи сложите парные частоты, сохранив одиночные ячейки.
  5. Для оценки мощности полосы умножьте значения плотности внутри полосы на частотный шаг и сложите их.

Историческая справка

Артур Шустер в работе 1898 года исследовал скрытые периодичности на примере предполагаемого 26-дневного цикла метеорологических явлений. Классическая периодограмма использовала частотный анализ для поиска таких компонент в ряду измерений.

В цифровой обработке частоты задаются сеткой ДПФ, а оконное взвешивание меняет форму спектральных пиков и масштаб оценки. Документация SciPy различает плотность спектра и спектр квадратов амплитуд, поскольку эти результаты имеют разные единицы и разные делители. Нормировка на частоту дискретизации и энергию окна относится к первому варианту. Масштаб в В²/Гц позволяет сравнивать оценки для разной частотной сетки при суммировании мощности по полосе.

Историческая линия формулы

Историческое название периодограммы связано с работой Шустера о скрытых периодичностях. Современная формула на странице описывает плотность для равномерно дискретизированной записи и следует соглашениям оконного спектрального анализа, изложенным в документации SciPy.

Пример

Даны отсчёты напряжения [3, 1, 3, 1] В при fₛ=8 Гц. Возьмём прямоугольное окно и сохраним среднее. Найдём одностороннюю периодограмму. ДПФ равно [8, 0, 4, 0] В. Знаменатель fₛN=32 Гц. Нулевая ячейка даёт 64/32=2 В²/Гц, внутренняя нулевая, а Найквист даёт 16/32=0,5 В²/Гц. Значения стоят на частотах 0, 2 и 4 Гц. Шаг равен 2 Гц. Сумма вкладов равна (2+0+0,5)·2=5 В². Проверка по отсчётам: (9+1+9+1)/4=5 В². Отдельно видны 4 В² постоянной составляющей и 1 В² чередования. Суммарное значение имеет единицы В² после умножения на герцы.

Частая ошибка

Удваивают частоту Найквиста как обычную пару. Применяют сумму весов вместо суммы их квадратов и получают другую шкалу. Сравнивают плотность одной ячейки с мощностью полосы без учёта Δf. Воспроизводят пример библиотекой с удалением среднего, хотя в примере оно сохранено.

Практика

Задачи с решением

Пара частот вещественного сигнала

Условие. x=[1, 0, −1, 0] В, fₛ=8 Гц, окно прямоугольное. Найдите плотность в односторонней ячейке 2 Гц.

Решение. Коэффициенты на частотах ±2 Гц имеют модули 2 В. Двусторонняя плотность каждой 4/32=0,125 В²/Гц. Складываем пару.

Ответ. 0,25 В²/Гц.

Нормировка по энергии окна

Условие. fₛ=10 Гц, веса [0, 1, 1, 0], модуль Xw[1]=6 В для вещественного сигнала. Найдите одностороннюю плотность в этой ячейке.

Решение. Энергия окна 2. Двусторонняя оценка 36/(10·2)=1,8 В²/Гц. При N=4 индекс 1 парный, поэтому удваиваем.

Ответ. 3,6 В²/Гц.

Дополнительные источники

  • https://doi.org/10.1029/TM003i001p00013
  • https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.periodogram.html
  • https://docs.scipy.org/doc/scipy/tutorial/signal.html

Связанные формулы

Информатика

Равенство Парсеваля для энергии блока ДПФ

$\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2$

Равенство Парсеваля позволяет вычислить сумму квадратов отсчётов по всем коэффициентам ДПФ. Множитель 1/N применяется, если прямое преобразование записано без нормировки.

Информатика

Циклическая свёртка через дискретное преобразование Фурье

$\begin{aligned}y[n]&=\sum_{m=0}^{N-1}x[m]\,h[(n-m)\bmod N],\\Y[k]&=X[k]H[k].\end{aligned}$

Циклическая свёртка возвращает N значений, обращаясь к индексам по модулю N. Она равна обратному ДПФ произведения спектров, а дополнение нулями позволяет получить линейную свёртку.

Информатика

Спектральная плотность мощности по методу Уэлча

$\widehat S_{\mathrm{Welch}}[k]=\frac{1}{M}\sum_{j=0}^{M-1}P_j[k],\quad P_j[k]=\frac{|\operatorname{DFT}(w\,x_j)[k]|^2}{f_s\sum_n w[n]^2}$

Метод Уэлча разбивает запись на сегменты, строит оконную периодограмму каждого сегмента и усредняет значения плотности по частотам. Сегменты могут перекрываться.

Информатика

Дискретное преобразование Фурье и обратная формула

$X_k=\sum_{n=0}^{N-1}x_n e^{-2\pi i kn/N},\qquad x_n=\frac1N\sum_{k=0}^{N-1}X_k e^{2\pi i kn/N}$

Дискретное преобразование Фурье переводит конечную последовательность в комплексные частотные коэффициенты. Обратная сумма восстанавливает все отсчёты; знак экспоненты и нормировку нужно выбирать согласованно.

Информатика

Окно Ханна: коэффициенты и поправка амплитуды

$w_n=\frac12\left(1-\cos\frac{2\pi n}{N-1}\right),\quad G_c=\frac1N\sum_{n=0}^{N-1}w_n,\quad N>1$

Симметричное окно Ханна плавно уменьшает края записи до нуля. Оно снижает дальнюю спектральную утечку ценой расширения главного лепестка и изменения амплитуды, поэтому масштаб спектра требует отдельной поправки.