pwelch

Оценка спектральной плотности мощности Уэлча.

Библиотека

EngeeDSP

Синтаксис

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

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

pxx,f,pxxc = pwelch(x,win,nOverlap,freqSpec,out=:data) — возвращает оценку СПМ pxx, частоты f и границы доверительного интервала pxxc, где
- win — разбиение сигнала на сегменты для оконной обработки;
- nOverlap — длина перекрытия соседних сегментов;
- freqSpec — количество точек дискретного преобразования Фурье (ДПФ) или частоты, в которых оно вычисляется.
Чтобы использовать значение по умолчанию для любого из этих аргументов, укажите их пустыми, например Float64[].

pxx,f,pxxc = pwelch(x,win,nOverlap,freqSpec,Fs,out=:data) — задает частоту дискретизации Fs и возвращает циклические частоты f в Гц.

pxx,f,pxxc = pwelch(___,freqRange,spectrumType,trace,out=:data) — дополнительно задает диапазон частот freqRange, тип спектра spectrumType и режим трассировки trace. Можно указать любой из предыдущих вариантов входного синтаксиса.

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

pwelch(___) — строит график оценки СПМ или спектра мощности Уэлча.

Аргументы

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

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

Details

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

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

# win — окно
пустой массив (по умолчанию) | целое число | вектор

Details

Окно, заданное как пустой массив Float32[], Float64[] или Int64[], положительное целое число или вектор.

Функция pwelch делит входной сигнал x на сегменты и применяет окно для поэлементного умножения отсчетов каждого сегмента на значения окна. Значения окна зависят от способа задания этого аргумента:

Float32[], Float64[] или Int64[] — функция применяет окно Хэмминга такой длины, чтобы разделить x на 8 сегментов с перекрытием в nOverlap отсчетов;
положительное целое число — функция делит x на сегменты длиной win и применяет окно Хэмминга той же длины;
вектор — функция делит x на сегменты длиной, равной длине вектора win, и применяет окно, заданное в аргументе win.

Если длина x не позволяет получить целое число сегментов с перекрытием nOverlap отсчетов, функция pwelch соответствующим образом усекает сигнал x.

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

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

Типы данных

Float32, Float64, Int64

# nOverlap — количество перекрывающихся отсчетов
пустой массив (по умолчанию) | целое число

Details

Количество перекрывающихся отсчетов, заданное как пустой массив Float32[], Float64[] или Int64[] или неотрицательное целое число.

Float32[], Float64[] или Int64[] — функция выбирает значение, обеспечивающее перекрытие сегментов на 50%.

Если длина сегмента не задана, функция устанавливает для этого аргумента значение , где — длина входного сигнала, а обозначает функцию округления вниз. Это действие эквивалентно делению сигнала на максимально длинные сегменты, чтобы получить как можно больше, но не более восьми сегментов с перекрытием 50%.
Неотрицательное целое число — смежные сегменты перекрываются на заданное в этом аргументе количество отсчетов.
- Если win — скаляр, значение nOverlap должно быть меньше win.
- Если win — вектор, значение nOverlap должно быть меньше длины win.

# freqSpec — частотная характеристика
пустой массив (по умолчанию) | целое число | вектор

Details

Частотная характеристика, заданная одним из следующих способов:

Float32[], Float64[] или Int64[] — функция вычисляет спектральные оценки, используя количество точек ДПФ, равное большему из значений 256 и ближайшей степени двойки, равной или большей длине окна;
положительное целое число — функция вычисляет спектральные оценки, используя количество точек ДПФ, заданное в этом аргументе:
- если freqSpec больше длины сегмента, функция pwelch дополняет данные сегмента нулями;
- если freqSpec меньше длины сегмента, функция разбивает сегмент по модулю freqSpec и суммирует полученные фрагменты;
вектор, содержащий не менее двух элементов, — функция вычисляет спектральные оценки на частотах, заданных в этом аргументе:
- если задана частота дискретизации Fs, функция рассматривает значения freqSpec как циклические частоты в тех же единицах, что и Fs;
- если частота дискретизации не задана, функция предполагает, что значения freqSpec — это нормализованные частоты в рад/отсчет.

Типы данных

Float32, Float64, Int64

# Fs — частота дискретизации
скаляр

Details

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

Если аргумент Fs задан как Float32[], Float64[] или Int64[], функция pwelch принимает частоту дискретизации входного сигнала x равной 1 Гц.
Если аргумент Fs не задан, функция pwelch принимает частоту дискретизации равной 2π рад/отсчет.

# freqRange — диапазон частот для оценки СПМ или спектра мощности
"onesided" | "twosided" | "centered"

Details

Диапазон частот оценки СПМ или спектра мощности, заданный как "onesided", "twosided" или "centered".

Функция pwelch возвращает pxx с количеством строк и частотным интервалом, зависящим от значения freqRange, четности количества точек ДПФ в freqSpec и от того, задан ли аргумент Fs.

freqRange

freqSpec

Количество строк в pxx

Частотный интервал для pxx

Fs не задан

Fs задан

"onesided" (по умолчанию, если x — вещественное число)

Четное

freqSpec / 2 + 1

[0, π] рад/отсчет

[0, Fs/2] циклов в единицу времени

Нечетное

(freqSpec + 1) / 2

[0, π) рад/отсчет

[0, Fs/2) циклов в единицу времени

"twosided" (по умолчанию, если x — комплексное число)

Четное или нечетное

freqSpec

[0, 2π) рад/отсчет

[0, Fs) циклов в единицу времени

"centered"

Четное

freqSpec

(−π, π] рад/отсчет

(−Fs/2, Fs/2] циклов в единицу времени

Нечетное

(−π, π) рад/отсчет

(−Fs/2, Fs/2) циклов в единицу времени

Этот аргумент не поддерживается, если freqSpec задан как вектор циклических или нормализованных частот.
Для значения "onesided" функция pwelch умножает мощность на 2 на всех частотах, кроме нулевой частоты и частоты Найквиста, чтобы сохранить полную мощность.
Значение "onesided" не поддерживается при комплексном значении x.

Типы данных

String

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

Details

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

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

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

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

Частота дискретизации	Соотношение масштабирования
`Fs` задан	`psd = pwelch(x, win, nOverlap, freqSpec, Fs, "psd") pow = pwelch(x, win, nOverlap, freqSpec, Fs, "power")` Тогда `pow` эквивалентно `psd * enbw(win, Fs)`, где `win` — это вектор окна.
`Fs` не задан	`psd = pwelch(x, win, nOverlap, freqSpec, "psd") pow = pwelch(x, win, nOverlap, freqSpec, "power")` Тогда `pow` эквивалентно `psd * enbw(win, 2 * π)`.

Fs задан

psd = pwelch(x, win, nOverlap, freqSpec, Fs, "psd")
pow = pwelch(x, win, nOverlap, freqSpec, Fs, "power")

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

Fs не задан

psd = pwelch(x, win, nOverlap, freqSpec, "psd")
pow = pwelch(x, win, nOverlap, freqSpec, "power")

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

# trace — режим трассировки
"mean" (по умолчанию) | "maxhold" | "minhold"

Details

Режим трассировки, заданный одним из следующих значений:

"mean" — функция pwelch возвращает оценку спектра Уэлча для каждого входного канала путем усреднения оценок спектра мощности всех сегментов в каждом частотном диапазоне;
"maxhold" — функция возвращает спектр максимального удержания для каждого входного канала, выбирая максимальное значение спектра мощности по сегментам в каждом частотном диапазоне;
"minhold" — функция возвращает спектр минимального удержания для каждого входного канала, выбирая минимальное значение спектра мощности по сегментам в каждом частотном диапазоне.

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

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

Используйте запятые для разделения имени и значения, а Name заключите в кавычки, либо используйте знак равенства для разделения имени и значения, а Name укажите без кавычек.

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

Details

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

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

# Parent — целевой родительский контейнер
пустой массив (по умолчанию) | объект Axes | объект Panel

Details

Целевой родительский контейнер, заданный как объект Axes или Panel.

Если задано Parent = h (или "Parent", h), функция pwelch строит график оценки СПМ или спектра мощности Уэлча на указанном целевом родительском контейнере независимо от наличия выходных аргументов.

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

Details

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

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

Для этого аргумента имя и значение разделяются знаком равенства (=).

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

# 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.

Типы данных

Float32, Float64, Int64

Примеры

Оценка Уэлча с использованием входных данных по умолчанию

Details

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

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

using Random
Random.seed!(0)

n = 0:319
x = cos.(π/4 * n) + randn(length(n))

Получим оценку СПМ Уэлча, используя окно Хэмминга и длину ДПФ по умолчанию. Длина сегмента по умолчанию составляет 71 отсчет, а длина ДПФ содержит 256 точек, что дает разрешение по частоте рад/отсчет. Поскольку сигнал является вещественным, оценка СПМ является односторонней, и имеется 256 / 2 + 1 точек.

import EngeeDSP.Functions: pwelch

pxx, f, pxxc = pwelch(x, out = :data)

([0.1863630724000894; 0.33918260379618265; … ; 0.24304012258269603; 0.11708917663794487;;], [0.0; 0.02454369260617026; … ; 3.117048960983623; 3.141592653589793;;], [0.08146885100079476 0.7719778085293696; 0.1818050422612173 0.8436289786035511; … ; 0.130271774739574 0.6044994291550672; 0.05118568052389489 0.48502229985464423])

Построим график оценки СПМ Уэлча.

pwelch(x)

pwelch 1

Повторим вычисления.

Разделим сигнал на участки длиной . Это действие эквивалентно разделению сигнала на максимально длинные сегменты, чтобы получить как можно больше, но не более 8 сегментов с перекрытием 50%.
Используем окно Хэмминга для сегментов.
Укажем перекрытие 50% между смежными сегментами.
Для вычисления БПФ используем точек, где .

Убедимся, что оба подхода дают идентичные результаты.

Nx = length(x)
nsc = floor(Int, Nx / 4.5)
nov = floor(Int, nsc / 2)
nff = max(256, 2^nextpow(2, nsc))

import EngeeDSP.Functions: hamming

pxxt = pwelch(x, hamming(nsc), nov, nff, out = :data)

maxerr = maximum(abs.(abs.((pxxt[1])[:]) - abs.(pxx[:])))

0.0

Разделим сигнал на 8 равных сегментов с перекрытием 50% между ними. Укажем ту же длину БПФ, что и на предыдущем шаге. Вычислим оценку СПМ Уэлча и убедимся, что результат совпадает с предыдущими оценками.

ns = 8
ov = 0.5
lsc = Int(floor(Nx / (ns - (ns - 1) * ov)))

pxxt8 = pwelch(x, lsc, floor(ov * lsc), nff, out = :data)

maxerr8 = maximum(abs.(abs.((pxxt8[1])[:]) - abs.(pxx[:])))

0.0

Оценка Уэлча с заданной длиной сегмента

Details

Создадим синусоиду с угловой частотой рад/отсчет и аддитивным белым шумом . Сбросим состояние генератора случайных чисел для получения воспроизводимых результатов. Сигнал содержит 512 отсчетов.

using Random
Random.seed!(0)

n = 0:511
x = cos.(π/3 * n) + randn(length(n))

Получим оценку СПМ Уэлча, разделив сигнал на сегменты длиной 132 отсчета. Сегменты сигнала умножаются на окно Хэмминга длиной 132 отсчета. Количество перекрывающихся отсчетов не указано, поэтому оно установлено равным 132 / 2 = 66. Длина ДПФ составляет 256 точек, что дает разрешение по частоте 2π / 256 рад/отсчет. Поскольку сигнал является вещественным, оценка СПМ является односторонней, и имеется 256 / 2 + 1 = 129 точек. Построим график СПМ как функцию нормализованной частоты.

import EngeeDSP.Functions: pwelch

segmentLength = 132
pxx, w, pxxc = pwelch(x, segmentLength, out = :data)

plot(w/π, 10 * log10.(pxx),
     xlabel = "Normalized Frequency (× π rad/sample)",
     title = "Welch Power Spectral Density Estimate")

pwelch 2

Оценка Уэлча с заданным перекрытием сегментов

Details

using Random
Random.seed!(0)

n = 0:319
x = cos.(π/4 * n) + randn(length(n))

Получим оценку СПМ Уэлча, разделив сигнал на сегменты длиной 100 отсчетов. Сегменты сигнала умножаются на окно Хэмминга длиной 100 отсчетов. Количество перекрывающихся отсчетов равно 25. Длина ДПФ составляет 256 точек, что дает разрешение по частоте 2π / 256 рад/отсчет. Поскольку сигнал является вещественным, оценка СПМ является односторонней, и имеется 256 / 2 + 1 точек.

import EngeeDSP.Functions: pwelch

segmentLength = 100
noverlap = 25
pxx, f, pxxc = pwelch(x, segmentLength, noverlap, out = :data)

plot(10 * log10.(pxx),
     xlabel = "Frequency Samples",
     title = "Welch Power Spectral Density Estimate")

pwelch 3

Оценка Уэлча с заданной длиной ДПФ

Details

using Random
Random.seed!(0)

n = 0:319
x = cos.(π/4 * n) + randn(length(n))

Получим оценку СПМ Уэлча, разделив сигнал на сегменты длиной 100 отсчетов с перекрытием по умолчанию 50%. Укажем длину ДПФ равной 640 точкам, чтобы частота π/4 рад/отсчет соответствовала 81-му бину ДПФ. Поскольку сигнал является вещественным, оценка СПМ является односторонней, и имеется 640 / 2 + 1 точек.

import EngeeDSP.Functions: pwelch

segmentLength = 100
nfft = 640
pxx, f, pxxc = pwelch(x, segmentLength, [], nfft, out = :data)

plot(10 * log10.(pxx),
     xlabel = "Frequency Points",
     title = "Welch Power Spectral Density Estimate")

pwelch 4

Оценка СПМ сигнала по методу Уэлча с частотой в Гц

Details

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

using Random
Random.seed!(0)

Fs = 1000
t = 0:1/Fs:5-1/Fs
x = cos.(2π * 100 * t) + randn(length(t))

Получим оценку СПМ Уэлча с усреднением по перекрывающимся сегментам для предыдущего сигнала. Используем длину сегмента 500 отсчетов с перекрытием 300 отсчетов. Используем 500 точек ДПФ, чтобы частота 100 Гц попадала непосредственно в бин ДПФ. Введем частоту дискретизации, чтобы получить вектор частот в Гц. Построим график результата.

import EngeeDSP.Functions: pwelch

pxx, f, pxxc = pwelch(x, 500, 300, 500, Fs, out = :data)

plot(f, 10 * log10.(pxx),
     xlabel = "Frequency (Hz)",
     ylabel = "PSD (dB/Hz)")

pwelch 5

Верхняя и нижняя границы доверительного интервала

Details

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

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

using Random
Random.seed!(0)

Fs = 1000
t = 0:1/Fs:1-1/Fs
x = cos.(2π * 100 * t) + sin.(2π * 150 * t) + randn(length(t))

Получим оценку с уровнем доверия 95%. Установим длину сегмента равной 200 отсчетам, а перекрытие — 50% (то есть 100 отсчетов).

import EngeeDSP.Functions: hamming, pwelch

L = 200
win = hamming(L)
nOverlap = 100

pxx, f, pxxc = pwelch(x, win, nOverlap, 200, Fs, ConfidenceLevel = 0.95, out = :data)

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

plot(f, 10 * log10.(pxx),
     xlim = (25, 250),
     xlabel = "Frequency (Hz)",
     ylabel = "PSD (dB/Hz)",
     title = "Welch Estimate with 95%-Confidence Bounds")

plot!(f, 10 * log10.(pxxc),
      linestyle = :dashdot,
      color = "#DD5400")

pwelch 6

Спектр мощности с центром в нулевой частоте

Details

using Random
Random.seed!(0)

Fs = 1000
t = 0:1/Fs:5-1/Fs

noisevar = 1/4
x = cos.(2π * 100 * t) + sqrt(noisevar) * randn(length(t))

Получим центрированный спектр мощности методом Уэлча, используя сегменты длиной 500 отсчетов, перекрытие 300 отсчетов и ДПФ длиной 500 точек.

import EngeeDSP.Functions: pwelch

pxx, f, pxxc = pwelch(x, 500, 300, 500, Fs, "centered", "power", out = :data)

Построим график результата. Мощность на частотах −100 и 100 Гц близка к ожидаемой мощности 1/4 для вещественной синусоиды с амплитудой 1. Отклонение от 1/4 обусловлено влиянием аддитивного шума.

plot(f, 10 * log10.(pxx),
     xlabel = "Frequency (Hz)",
     ylabel = "Magnitude (dB)")

pwelch 7

Оценка СПМ многоканального сигнала

Details

Сгенерируем 1024 отсчета многоканального сигнала, состоящего из трех синусоид с частотами , и рад/отсчет в аддитивном белом гауссовом шуме . Оценим СПМ Уэлча и построим график.

N = 1024
n = 0:N-1

w = π ./ [2, 3, 4]

using Random
Random.seed!(0)

x = cos.(w * collect(n)')' + randn(length(n), 3)

import EngeeDSP.Functions: pwelch

pwelch(x)

pwelch 8

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

Оценка СПМ методом усреднения перекрывающихся сегментов Уэлча

Details

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

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

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

Литература

Hayes, Monson H. Statistical Digital Signal Processing and Modeling. New York: John Wiley & Sons, 1996.
Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.