Документация Engee

periodogram

Оценка спектральной плотности мощности методом периодограммы.

Библиотека

EngeePhased

Синтаксис

Вызов функции

  • pxx,f = periodogram(x,out=:data) — возвращает оценку спектральной плотности мощности (СПМ) сигнала x, вычисленную методом периодограммы. Столбцы x рассматриваются как независимые каналы. Также функция возвращает частоты f в рад/отсчет.

  • pxx,f = periodogram(x,win,out=:data) — также использует win — оконную функцию, применяемую к сигналу.

  • pxx,f = periodogram(x,win,nfft,out=:data) — также использует nfft — число точек в дискретном преобразовании Фурье (ДПФ).

  • pxx,f = periodogram(x,win,freq,out=:data) — также использует freq — нормализованные или циклические частоты.

    Можно указать либо freq как вектор, либо nfft как целое положительное число. Одновременное использование этих аргументов приведет к ошибке.

  • pxx,f,pxxc = periodogram(___,Name=Value,out=:data) — задает дополнительные параметры с помощью аргументов типа «имя-значение» для любого из предыдущих вариантов синтаксиса. Также функция возвращает границы доверительного интервала pxxc, если для оценки спектральной плотности мощности указан уровень ConfidenceLevel=p.

  • periodogram(___) — строит график оценки СПМ или спектра мощности. Для этого синтаксиса можно дополнительно передавать параметры для отображения итогового графика. Например, если вызвать periodogram(x, color = :red), то график будет красного цвета. Если указан аргумент ConfidenceLevel, то для каждого столбца доверительные интервалы также будут отображены пунктирными линиями.

Аргументы win, nfft и freq можно задать как аргументы типа «имя-значение». Например, вызов periodogram(x, nfft = 512) аналогичен вызову periodogram(x, fill(1, size(x, 1)), 512). Это позволяет опустить некоторые аргументы, для которых есть реализация по умолчанию.

Аргументы

Входные аргументы

# x — входной сигнал
вектор | матрица

Details

Входной сигнал, заданный как вектор или матрица. Если x — матрица, функция periodogram рассматривает ее столбцы как независимые каналы.

Типы данных

Float32, Float64
Поддержка комплексных чисел

Да

# win — окно
fill(1, size(x, 1)) (по умолчанию) | вещественный или комплексный вектор

Details

Окно, заданное как вектор той же длины, что и входной сигнал. Если аргумент win не задан, функция periodogram использует прямоугольное окно.

Стандартное прямоугольное окно характеризуется уровнем подавления боковых лепестков, равным 13.3 дБ; это может привести к маскировке спектральных составляющих, уровень которых ниже данного значения (относительно пикового значения спектра). Выбор других типов окон позволяет найти компромисс между разрешающей способностью (например, при использовании прямоугольного окна) и степенью подавления боковых лепестков (например, при использовании окна Ханна).
Типы данных

Float32, Float64

# nfft — количество точек в ДПФ
max(256, nextpow(2, size(x, 1))) (по умолчанию) | целое положительное число

Details

Количество точек в ДПФ, заданное как целый положительный скаляр.

Если аргумент nfft не задан, то функция использует значение , где . Символы обозначают округление вверх, а length(win).

Типы данных

Float32, Float64

# freq — нормализованные или циклические частоты
вещественный вектор

Details

Значения частот, заданные как вещественный вектор. Тип задаваемых частот зависит от значения аргумента fs:

  • Если значение входного аргумента fs не задано, то частоты по умолчанию считаются нормализованными.

  • Если значение входного аргумента fs указано, то частоты считаются циклическими.

Можно задать либо freq как вектор, либо nfft как целое положительное число. Одновременное использование этих аргументов приведет к ошибке.

Входные аргументы «имя-значение»

Укажите необязательные пары аргументов в виде Name=Value, где Name — имя аргумента, а Value — соответствующее значение. Аргументы типа «имя-значение» должны располагаться после других аргументов, но порядок пар не имеет значения. Можно указать несколько пар «имя-значение».

# fs — частота дискретизации
NaN (по умолчанию) | положительный скаляр

Details

Частота дискретизации, заданная как положительный скаляр. Частота дискретизации — это количество отсчетов в единицу времени. Если единицей времени являются секунды, то частота дискретизации выражается в Гц.

# freqrange — частотный диапазон оценки СПМ
"onesided" | "twosided"

Details

Частотный диапазон для оценки СПМ, заданный как "onesided" или "twosided".

Для сигналов с вещественными значениями по умолчанию используется значение "onesided". Для комплексных сигналов по умолчанию используется значение "twosided", а указание значения "onesided" приводит к ошибке.

  • Значение "onesided" — возвращает одностороннюю периодограмму вещественного входного сигнала.

    • Если значение аргумента nfft четное, то аргумент pxx имеет nfft/2 + 1 строк и вычисляется на интервале [0, π] радиан на отсчет.

    • Если значение аргумента nfft нечетное, то pxx имеет (nfft + 1)/2 строк и вычисляется на интервале [0, π) радиан на отсчет. Если вы укажете fs, то интервалы будут соответственно [0, fs/2] циклов на единицу времени и [0, fs/2) циклов на единицу времени.

  • Значение "twosided" — возвращает двустороннюю периодограмму вещественного или комплексного сигнала. Аргумент pxx имеет nfft строк и вычисляется на интервале [0, 2π) радиан на отсчет. Если вы укажете аргумент fs, то интервал составит [0, fs) циклов на единицу времени.

# centered — является ли частотный диапазон центрированным
false (по умолчанию) | true

Details

Включение центрирования частотного диапазона, заданное как false или true.

# spectrumtype — масштабирование спектра мощности
"psd" (по умолчанию) | "power"

Details

Масштабирование спектра мощности, заданное одним из следующих значений:

  • "psd" — функция periodogram возвращает спектральную плотность мощности;

  • "power" — функция масштабирует каждую оценку СПМ на эквивалентную шумовую полосу частот окна и возвращает оценку мощности на каждой частоте.

В таблице показано соотношение масштабирования между оценкой СПМ и оценкой спектра мощности, возвращаемой в pxx, при заданных входном сигнале x, векторе окна win, количестве точек в ДПФ nfft и частоте дискретизации fs.

Частота дискретизации Соотношение масштабирования

fs задан

 
psd = periodogram(x, win, nfft, fs=Fs, spectrumtype="psd")
pow = periodogram(x, win, nfft, fs=Fs, spectrumtype="power")

Тогда pow эквивалентно psd * enbw(win, Fs), где win — это вектор окна.

fs не задан

 
psd = periodogram(x, win, nfft, spectrumtype="psd")
pow = periodogram(x, win, nfft, spectrumtype="power")

Тогда pow эквивалентно psd * enbw(win, 2 * π).

# ConfidenceLevel — доверительный интервал для оценки СПМ
NaN (по умолчанию) | вещественное число в диапазоне от 0 до 1

Details

Вероятность покрытия для оценки СПМ, заданная как скаляр в диапазоне (0, 1).

Если задано ConfidenceLevel = p, функция periodogram возвращает в выходном аргументе pxxc нижнюю и верхнюю границы доверительного интервала с уровнем p × 100% для истинной СПМ.

# out — тип выходных данных
:plot (по умолчанию) | :data

Details

Тип выходных данных:

  • :plot — функция возвращает график;

  • :data — функция возвращает данные.

Выходные аргументы

# pxx — оценка СПМ или спектра мощности
вектор | матрица

Details

Оценка СПМ или спектра мощности, возвращаемая в виде вещественного неотрицательного вектора-столбца или матрицы.

  • Каждый столбец pxx представляет собой оценку СПМ или спектр мощности соответствующего столбца x в зависимости от того, как задан аргумент spectrumtype.

  • Единицы измерения оценки СПМ — квадраты величины входного сигнала на единицу частоты. Например, если входные данные x заданы в вольтах, частота дискретизации fs — в герцах, а сопротивление равно 1 Ом, то оценка СПМ будет выражена в Вт/Гц.

  • Единицы измерения спектра мощности — квадраты величины входного сигнала. Например, если входные данные x заданы в вольтах, а сопротивление равно 1 Ом, то оценка спектра мощности будет выражена в ваттах.

# f — частоты, связанные с оценкой СПМ
вектор

Details

Частоты, связанные с оценкой СПМ, возвращаемые в виде вещественного вектора-столбца.

  • Если аргумент fs задан, f содержит циклические частоты в Гц.

  • Если аргумент fs не задан, f содержит нормализованные частоты в рад/отсчет.

# pxxc — доверительные интервалы
матрица

Details

Границы доверительных интервалов, возвращаемые в виде вещественной матрицы.

  • Матрица pxxc имеет столько же строк, сколько pxx.

  • Матрица pxxc имеет вдвое больше столбцов, чем pxx.

    • Нечетные столбцы содержат нижние границы доверительных интервалов.

    • Четные столбцы содержат верхние границы доверительных интервалов.

    Таким образом, pxxc[m, 2 * n − 1] — это нижняя граница доверительного интервала, а pxxc[m, 2 * n] — верхняя граница доверительного интервала, соответствующая оценке pxx[m, n].

  • Вероятность покрытия доверительных интервалов определяется значением аргумента ConfidenceLevel.

Зависимости

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

Типы данных

Float32, Float64

Примеры

Периодограмма со стандартными параметрами

Details

Получим периодограмму входного сигнала, представляющего собой дискретную синусоиду с угловой частотой рад/отсчет, к которой добавлен белый шум с распределением .

Сформируем входной сигнал. Длина сигнала составляет 320 отсчетов. Вычислим периодограмму, используя прямоугольное окно и длину ДПФ, заданные по умолчанию. Длина ДПФ выбирается равной ближайшей степени двойки, превышающей длину сигнала, то есть 512 точкам. Поскольку сигнал является вещественным и имеет четную длину, периодограмма получается односторонней и содержит 512/2 + 1 точку.

import EngeePhased.Functions: periodogram

n = collect(0:319);
x = cos.(pi ./ 4 .* n) .+ randn(size(n));
periodogram(x)

periodogram 1

Модифицированная периодограмма с окном Хэмминга

Details

Получим модифицированную периодограмму входного сигнала, представляющего собой дискретную синусоиду с угловой частотой рад/отсчет, к которой добавлен белый шум с распределением .

Сформируем входной сигнал. Длина сигнала составляет 320 отсчетов. Вычислим периодограмму, используя прямоугольное окно и длину ДПФ, заданные по умолчанию. Вычислим модифицированную периодограмму, используя окно Хэмминга и длину ДПФ по умолчанию. Длина ДПФ выбирается равной ближайшей степени двойки, превышающей длину сигнала, то есть 512 точкам. Поскольку сигнал является вещественным и имеет четную длину, периодограмма получается односторонней и содержит 512/2 + 1 точку.

import EngeePhased.Functions: periodogram
import EngeeDSP.Functions: hamming

n = collect(0:319)
x = cos.(pi ./ 4 .* n) .+ randn(size(n));
periodogram(x,hamming(length(x)))

periodogram 2

Длина ДПФ совпадает с длиной сигнала

Details

Получим периодограмму входного сигнала, представляющего собой дискретную синусоиду с угловой частотой рад/отсчет, к которой добавлен белый шум с распределением . Используем длину ДПФ, равную длине сигнала.

Сформируем входной сигнал. Длина сигнала составляет 320 отсчетов. Вычислим периодограмму, используя прямоугольное окно (по умолчанию) и длину ДПФ, равную длине сигнала. Поскольку сигнал является вещественным, по умолчанию возвращается односторонняя периодограмма длиной 320/2 + 1.

import EngeePhased.Functions: periodogram

n = collect(0:319);
x = cos.(pi ./4 .* n) .+ randn(size(n));
nfft = length(x);
periodogram(x,nfft=nfft)

periodogram 3

Периодограмма на наборе нормализованных частот

Details

Получим периодограмму входного сигнала, состоящего из двух дискретных синусоид с угловыми частотами и рад/отсчет на фоне аддитивного белого шума с распределением . Получим двусторонние оценки периодограммы на частотах и рад/отсчет.

import EngeePhased.Functions: periodogram

n = collect(0:319);
x = cos.(pi ./ 4 .* n) .+ 0.5 .* sin.(pi ./ 2 .* n) .+ randn(size(n));

pxx,w = periodogram(x,freq=[pi/4, pi/2],out=:data)
pxx
2×1 Matrix{Float64}:
 10.087037418886014
  1.7539619680699419

Сравним полученный результат с односторонней периодограммой. Значения двусторонней периодограммы составляют половину значений односторонней периодограммы. При вычислении периодограммы для определенного набора частот результатом является двусторонняя оценка.

periodogram 3

pxx1,w1 = periodogram(x,out=:data)

plot(w1 ./ π, pxx1; label="pxx1")
plot!(w ./ π, 2 .* pxx; label="2*pxx", seriestype=:scatter)
xlabel!("Normalized Frequency (× π rad/sample)")
title!("Periodogram Power Spectral Density Estimate")

periodogram 4

Периодограмма на наборе циклических частот

Details

Сформируем сигнал, состоящий из двух синусоид с частотами 100 и 200 Гц на фоне аддитивного белого шума с распределением . Частота дискретизации составляет 1 кГц. Вычислим двустороннюю периодограмму для частот 100 и 200 Гц.

import EngeePhased.Functions: periodogram

fs = 1000;
t = collect(0:0.001:1-0.001);
x = cos.(2 .* pi .* 100 .* t) .+ sin.(2 .* pi .* 200 .* t) .+ randn(size(t));

freq = [100, 200];
pxx,w = periodogram(x,freq=freq,fs=fs, out=:data)
pxx
2×1 Matrix{Float64}:
 0.28289790057149083
 0.26824883574824315

Центрированная периодограмма

Details

Получим периодограмму синусоидального сигнала с частотой 100 Гц на фоне аддитивного белого шума с распределением . Частота дискретизации составляет 1 кГц. Используем опцию центрирования для получения периодограммы, центрированной относительно нулевой частоты.

import EngeePhased.Functions: periodogram

fs = 1000;
t = collect(0:0.001:1-0.001);
x = cos.(2 .* pi .* 100 .* t) .+ randn(size(t));
periodogram(x,nfft=length(x), fs=fs, freqrange="onesided", centered=true)

periodogram 5

Доверительные интервалы 95%

Details

Сформируем сигнал, состоящий из двух синусоид с частотами 100 и 150 Гц на фоне аддитивного белого шума с распределением . Амплитуда обоих синусоидальных сигналов равна 1. Частота дискретизации составляет 1 кГц. Получим оценку СПМ методом периодограммы с 95%-ными доверительными интервалами.

import EngeePhased.Functions: periodogram

fs = 1000;
t = collect(0:1/fs:1-1/fs);
x = cos.(2 .* pi .* 100 .* t) .+ sin.(2 .* pi .* 150 .* t) + randn(size(t));
periodogram(x, nfft=length(x), fs=fs, ConfidenceLevel = 0.95)

periodogram 6

Оценка мощности синусоиды

Details

Оценим мощность синусоидального сигнала на заданной частоте, используя опцию "power".

Сформируем синусоидальный сигнал с частотой 100 Гц и длительностью 1 секунда при частоте дискретизации 1 кГц. Амплитуда синусоиды составляет 1.8, что соответствует мощности 1.82/2 = 1.62.

import EngeePhased.Functions: periodogram
import EngeeDSP.Functions: hamming

fs = 1000;
t = collect(0:1/fs:1-1/fs);
x = 1.8 .* cos.(2 .* pi .* 100 .* t);
pxx,f = periodogram(x, hamming(length(x)), length(x), fs=fs, spectrumtype="power", out = :data);
idx = argmax(pxx)
pwrest = pxx[idx]
println("The maximum power occurs at $(f[idx]) Hz")
println("The power estimate is $pwrest")
The maximum power occurs at 100.0 Hz
The power estimate is 1.6200000517531614

Дополнительно

Периодограмма

Details

Периодограмма представляет собой непараметрическую оценку спектральной плотности мощности (СПМ) стационарного в широком смысле случайного процесса. Периодограмма — это преобразование Фурье смещенной оценки последовательности автокорреляции. Для сигнала , дискретизированного с частотой fs отсчетов в единицу времени, периодограмма определяется следующим образом:

где — интервал дискретизации. В случае односторонней периодограммы значения на всех частотах, кроме 0 и частоты Найквиста , умножаются на 2 для сохранения общей мощности.

Если частоты выражены в радианах на отсчет, периодограмма определяется как

Диапазон частот в приведенных выше уравнениях может варьироваться в зависимости от значения аргумента freqrange.

Интеграл от истинной СПМ по одному периоду (равному для циклической частоты и для нормализованной частоты) равен дисперсии стационарного в широком смысле случайного процесса:

Для нормализованных частот пределы интегрирования следует изменить соответствующим образом.

Модифицированная периодограмма

Details

При вычислении модифицированной периодограммы входной временной ряд умножается на оконную функцию. Подходящая оконная функция является неотрицательной и убывает до нуля в начальной и конечной точках. Умножение временного ряда на оконную функцию обеспечивает плавное нарастание и спад значений данных, что помогает уменьшить эффект растекания в периодограмме.

Если — оконная функция, то модифицированная периодограмма определяется следующим образом:

где — интервал дискретизации.

Если частоты выражены в радианах на отсчет, модифицированная периодограмма определяется как

Диапазон частот в приведенных выше уравнениях может варьироваться в зависимости от значения аргумента freqrange.

Литература

  1. Auger, François, and Patrick Flandrin. Improving the Readability of Time-Frequency and Time-Scale Representations by the Reassignment Method. IEEE® Transactions on Signal Processing. Vol. 43, May 1995, pp. 1068–1089.

  2. Fulop, Sean A., and Kelly Fitz. Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications. Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.

Примеры