pwelch

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

Описание

пример

pxx = pwelch(x) возвращает оценку степени спектральной плотности (PSD), pxx, из входного сигнала, x, найденное использование перекрытого средства оценки усреднения сегмента валлийцев. Когда x вектор, он обработан как один канал. Когда x матрица, PSD вычисляется независимо для каждого столбца и хранится в соответствующем столбце pxx. Если x с действительным знаком, pxx односторонняя оценка PSD. Если x с комплексным знаком, pxx двухсторонняя оценка PSD. По умолчанию, x разделен на самые длинные сегменты, чтобы получить как близко к, но не превысить 8 сегментов с 50%-м перекрытием. Каждый сегмент является оконным с Окном Хэмминга. Модифицированные периодограммы усреднены, чтобы получить оценку PSD. Если вы не можете разделить длину x точно в целое число сегментов с 50%-м перекрытием, x является усеченным соответственно.

пример

pxx = pwelch(x,window) использует входной вектор или целое число, window, разделить сигнал на сегменты. Если window вектор, pwelch делит сигнал на сегменты, равные в длине к длине window. Модифицированные периодограммы вычисляются с помощью сегментов сигнала, умноженных на вектор, window. Если window целое число, сигнал разделен на сегменты длины window. Модифицированные периодограммы вычисляются с помощью Окна Хэмминга длины window.

пример

pxx = pwelch(x,window,noverlap) использование noverlap выборки перекрытия от сегмента до сегмента. noverlap должно быть положительное целое число, меньшее, чем window если window целое число. noverlap должно быть положительное целое число меньше, чем длина window если window вектор. Если вы не задаете noverlap, или задайте noverlap как пустое, количество по умолчанию перекрытых выборок составляет 50% длины окна.

пример

pxx = pwelch(x,window,noverlap,nfft) указывает, что количество дискретного преобразования Фурье (DFT) указывает, чтобы использовать в оценке PSD. nfft по умолчанию большие из 256 или следующая степень 2 больших, чем длина сегментов.

[pxx,w] = pwelch(___) возвращает нормированный вектор частоты, w. Если pxx односторонняя оценка PSD, w охватывает интервал [0, π] если nfft является четным и [0, π), если nfft является нечетным. Если pxx двухсторонняя оценка PSD, w охватывает интервал [0,2π).

пример

[pxx,f] = pwelch(___,fs) возвращает вектор частоты, f, в циклах в единицу времени. Частота дискретизации, fs, количество выборок в единицу времени. Если модуль времени является секундами, то f находится в циклах/секунда (Гц). Для сигналов с действительным знаком, f охватывает интервал [0, fs/2], когда nfft является четным и [0, fs/2) когда nfft является нечетным. Для сигналов с комплексным знаком, f охватывает интервал [0, fs). fs должен быть пятый вход к pwelch. Чтобы ввести частоту дискретизации и все еще использовать значения по умолчанию предыдущих дополнительных аргументов, задайте эти аргументы как пустые, [].

[pxx,w] = pwelch(x,window,noverlap,w) возвращает двухсторонние валлийские оценки PSD на нормированных частотах, заданных в векторе, w. Векторный w должен содержать по крайней мере два элемента, потому что в противном случае функция интерпретирует его как nfft.

[pxx,f] = pwelch(x,window,noverlap,f,fs) возвращает двухсторонние валлийские оценки PSD на частотах, заданных в векторе, f. Векторный f должен содержать по крайней мере два элемента, потому что в противном случае функция интерпретирует его как nfft. Частоты в f находятся в циклах в единицу времени. Частота дискретизации, fs, количество выборок в единицу времени. Если модуль времени является секундами, то f находится в циклах/секунда (Гц).

пример

[___] = pwelch(x,window,___,freqrange) возвращает валлийскую оценку PSD по частотному диапазону, заданному freqrange. Допустимые опции для freqrange : 'onesided', 'twosided', или 'centered'.

пример

[___] = pwelch(x,window,___,trace) возвращает оценку спектра хранения максимум если trace задан как 'maxhold' и возвращает оценку спектра хранения минимум если trace задан как 'minhold'.

пример

[___,pxxc] = pwelch(___,'ConfidenceLevel',probability) возвращает probability × 100% доверительных интервалов для PSD оценивают в pxxc.

пример

[___] = pwelch(___,spectrumtype) возвращает оценку PSD если spectrumtype задан как 'psd' и возвращает спектр мощности если spectrumtype задан как 'power'.

пример

pwelch(___) без выходных аргументов строит валлийскую оценку PSD в окне текущей фигуры.

Примеры

свернуть все

Получите валлийскую оценку PSD входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум.

Создайте синусоиду с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум. Сбросьте генератор случайных чисел для восстанавливаемых результатов. Сигнал имеет длину Nx=320 выборки.

rng default

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

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

pxx = pwelch(x);

pwelch(x)

Повторите расчет.

  • Разделите сигнал на разделы длины nsc=Nx/4.5. Это действие эквивалентно делению сигнала в самые длинные сегменты, чтобы получить как близко к, но не превысить 8 сегментов с 50%-м перекрытием.

  • Окно разделы с помощью Окна Хэмминга.

  • Задайте 50%-е перекрытие между непрерывными разделами

  • Чтобы вычислить БПФ, использовать max(256,2p) точки, где p=log2nsc.

Проверьте, что два подхода дают идентичные результаты.

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

t = pwelch(x,hamming(nsc),nov,nff);

maxerr = max(abs(abs(t(:))-abs(pxx(:))))
maxerr = 0

Разделите сигнал на 8 разделов равной длины с 50%-м перекрытием между разделами. Задайте ту же длину БПФ как на предыдущем шаге. Вычислите валлийскую оценку PSD и проверьте, что она дает тот же результат как предыдущие две процедуры.

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

t = pwelch(x,lsc,floor(ov*lsc),nff);

maxerr = max(abs(abs(t(:))-abs(pxx(:))))
maxerr = 0

Получите валлийскую оценку PSD входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой π/3 рад/выборка с дополнением N(0,1) белый шум.

Создайте синусоиду с угловой частотой π/3 рад/выборка с дополнением N(0,1) белый шум. Сбросьте генератор случайных чисел для восстанавливаемых результатов. Сигнал имеет 512 выборок.

rng default

n = 0:511;
x = cos(pi/3*n)+randn(size(n));

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

segmentLength = 132;
[pxx,w] = pwelch(x,segmentLength);

plot(w/pi,10*log10(pxx))
xlabel('\omega / \pi')

Получите валлийскую оценку PSD входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум.

Создайте синусоиду с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум. Сбросьте генератор случайных чисел для восстанавливаемых результатов. Сигнал является 320 выборками в длине.

rng default

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

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

segmentLength = 100;
noverlap = 25;
pxx = pwelch(x,segmentLength,noverlap);

plot(10*log10(pxx))

Получите валлийскую оценку PSD входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум.

Создайте синусоиду с угловой частотой π/4 рад/выборка с дополнением N(0,1) белый шум. Сбросьте генератор случайных чисел для восстанавливаемых результатов. Сигнал является 320 выборками в длине.

rng default

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

Получите валлийскую оценку PSD, делящую сигнал на сегменты 100 выборок в длине. Используйте перекрытие по умолчанию 50%. Задайте длину ДПФ, чтобы быть 640 точками так, чтобы частота π/4 рад/выборка соответствует интервалу ДПФ (интервал 81). Поскольку сигнал с действительным знаком, оценка PSD является односторонней и существуют точки 640/2+1.

segmentLength = 100;
nfft = 640;
pxx = pwelch(x,segmentLength,[],nfft);

plot(10*log10(pxx))
xlabel('rad/sample')
ylabel('dB / (rad/sample)')

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

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));

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

[pxx,f] = pwelch(x,500,300,500,fs);

plot(f,10*log10(pxx))

xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

Создайте сигнал, состоящий из трех шумных синусоид и щебета, произведенного на уровне 200 кГц в течение 0,1 секунд. Частоты синусоид составляют 1 кГц, 10 кГц и 20 кГц. Синусоиды имеют различные амплитуды и уровень шума. Бесшумный щебет имеет частоту, которая запускается на уровне 20 кГц и увеличивается линейно до 30 кГц во время выборки.

Fs = 200e3; 
Fc = [1 10 20]'*1e3; 
Ns = 0.1*Fs;

t = (0:Ns-1)/Fs;
x = [1 1/10 10]*sin(2*pi*Fc*t)+[1/200 1/2000 1/20]*randn(3,Ns);
x = x+chirp(t,20e3,t(end),30e3);

Вычислите валлийскую оценку PSD и хранение максимум и содержите минимум спектры сигнала. Постройте график результатов.

[pxx,f] = pwelch(x,[],[],[],Fs);
pmax = pwelch(x,[],[],[],Fs,'maxhold');
pmin = pwelch(x,[],[],[],Fs,'minhold');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
legend('pwelch','maxhold','minhold')

Повторите процедуру, на этот раз вычислив сосредоточенные оценки спектра мощности.

[pxx,f] = pwelch(x,[],[],[],Fs,'centered','power');
pmax = pwelch(x,[],[],[],Fs,'maxhold','centered','power');
pmin = pwelch(x,[],[],[],Fs,'minhold','centered','power');

plot(f,pow2db(pxx))
hold on
plot(f,pow2db([pmax pmin]),':')
hold off
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
legend('pwelch','maxhold','minhold')

Этот пример иллюстрирует использование доверительных границ с Перекрытым усреднением сегмента валлийцев (WOSA) оценка PSD. В то время как не необходимое условие для статистического значения, частоты в оценке валлийцев, где более низкая доверительная граница превышает верхнюю доверительную границу для окружения оценок PSD ясно, указывают на значительные колебания во временных рядах.

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

rng default
fs = 1000;
t = 0:1/fs:1-1/fs;
x = cos(2*pi*100*t)+sin(2*pi*150*t)+randn(size(t));

Получите оценку WOSA с 95%-доверительными-границами. Установите длину сегмента, равную 200, и перекройте сегменты на 50% (100 выборок). Постройте WOSA оценка PSD наряду с доверительным интервалом и увеличьте масштаб необходимой области частоты около 100 и 150 Гц.

L = 200;
noverlap = 100;
[pxx,f,pxxc] = pwelch(x,hamming(L),noverlap,200,fs,...
    'ConfidenceLevel',0.95);

plot(f,10*log10(pxx))
hold on
plot(f,10*log10(pxxc),'-.')
hold off

xlim([25 250])
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
title('Welch Estimate with 95%-Confidence Bounds')

Более низкая доверительная граница в мгновенной близости 100 и 150 Гц значительно выше верхней доверительной границы вне близости 100 и 150 Гц.

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

rng default

fs = 1000;
t = 0:1/fs:5-1/fs;

noisevar = 1/4;
x = cos(2*pi*100*t)+sqrt(noisevar)*randn(size(t));

Получите сосредоточенный DC спектр мощности с помощью метода валлийцев. Используйте продолжительность сегмента 500 выборок с 300 перекрытыми выборками и длину ДПФ 500 точек. Постройте результат.

[pxx,f] = pwelch(x,500,300,500,fs,'centered','power');

plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
grid

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

Сгенерируйте 1 024 выборки многоканального сигнала, состоящего из трех синусоид в дополнении N(0,1) белый Гауссов шум. Частоты синусоид π/2, π/3, и π/4 рад/выборка. Оцените PSD сигнала с помощью метода валлийцев и постройте его.

N = 1024;
n = 0:N-1;

w = pi./[2;3;4];
x = cos(w*n)' + randn(length(n),3);

pwelch(x)

Входные параметры

свернуть все

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

Пример: cos(pi/4*(0:159))+randn(1,160) одноканальный сигнал вектора-строки.

Пример: cos(pi./[4;2]*(0:159))'+randn(160,2) двухканальный сигнал.

Типы данных: single | double
Поддержка комплексного числа: Да

Окно в виде строки или вектор-столбца или целого числа. Если window вектор, pwelch делит x в перекрывающиеся сегменты длины равняются длине window, и затем умножает каждый сегмент сигнала с вектором, заданным в window. Если window целое число, pwelch разделен на сегменты длины, равной целочисленному значению, и Окно Хэмминга равной длины используется. Если длина x не может быть разделен точно в целое число сегментов с noverlap количество перекрывающихся выборок, x является усеченным соответственно. Если вы задаете window как пустое, Окно Хэмминга по умолчанию используется, чтобы получить восемь сегментов x с noverlap перекрывающиеся выборки.

Типы данных: single | double

Количество перекрытых выборок в виде положительного целого числа, меньшего, чем длина window. Если вы не используете noverlap или задайте noverlap как пустое, значение используется, чтобы получить 50%-е перекрытие между сегментами.

Количество ДПФ указывает в виде положительного целого числа. Для входного сигнала с действительным знаком, x, оценка PSD, pxx имеет длину (nfft/2 + 1), если nfft является четным, и (nfft + 1)/2, если nfft является нечетным. Для входного сигнала с комплексным знаком, x, оценка PSD всегда имеет длину nfft. Если nfft задан как пустой, nfft по умолчанию используется.

Если nfft больше длины сегмента, данные дополнены нулем. Если nfft меньше длины сегмента, сегмент перенесен с помощью datawrap заставить длину равняться nfft.

Типы данных: single | double

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

Нормированные частоты в виде строки или вектор-столбца по крайней мере с двумя элементами. Нормированные частоты находятся в рад/выборке.

Пример: w = [pi/4 pi/2]

Типы данных: double

Частоты в виде строки или вектор-столбца по крайней мере с двумя элементами. Частоты находятся в циклах в единицу времени. Единица времени задана частотой дискретизации, fs. Если fs имеет модули выборок/секунда, затем f имеет модули Гц.

Пример: fs = 1000; f = [100 200]

Типы данных: double

Частотный диапазон для PSD оценивает в виде того 'onesided', 'twosided', или 'centered'. Значением по умолчанию является 'onesided' для сигналов с действительным знаком и 'twosided' для сигналов с комплексным знаком. Частотные диапазоны, соответствующие каждой опции,

  • 'onesided' — возвращает одностороннюю оценку PSD входного сигнала с действительным знаком, x. Если nfft является четным, pxx имеет длину nfft/2 + 1 и вычисляется на интервале [0, π] рад/выборка. Если nfft является нечетным, длина pxx (nfft + 1)/2 и интервал [0, π), рад/выборка. Когда fs опционально задан, соответствующие интервалы [0, fs/2] циклы/единица времени и [0, fs/2) циклы/единица времени для четной и нечетной длины nfft соответственно.

  • 'twosided' — возвращает двухстороннюю оценку PSD или для входа с комплексным знаком или для с действительным знаком, x. В этом случае, pxx имеет длину nfft и вычисляется на интервале [0,2π), рад/выборка. Когда fs опционально задан, интервал [0, fs) циклы/единица времени.

  • 'centered' — возвращает двухстороннюю оценку PSD в центре или для входа с комплексным знаком или для с действительным знаком, x. В этом случае, pxx имеет длину nfft и вычисляется на интервале (–π, π] рад/выборка для даже длины nfft и (–π, π) рад/выборка для нечетной длины nfft. Когда fs опционально задан, соответствующие интервалы (–fs/2, fs/2] циклы/единица времени и (–fs/2, fs/2) циклы/единица времени для четной и нечетной длины nfft соответственно.

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

Проследите режим в виде одного из 'mean', 'maxhold', или 'minhold'. Значением по умолчанию является 'mean'.

  • 'mean' — возвращает валлийскую оценку спектра каждого входного канала. pwelch вычисляет валлийскую оценку спектра в каждом интервале частоты путем усреднения оценок спектра мощности всех сегментов.

  • 'maxhold' — возвращает спектр хранения максимум каждого входного канала. pwelch вычисляет спектр хранения максимум в каждом интервале частоты путем хранения максимального значения среди оценок спектра мощности всех сегментов.

  • 'minhold' — возвращает спектр хранения минимум каждого входного канала. pwelch вычисляет спектр хранения минимум в каждом интервале частоты путем хранения минимального значения среди оценок спектра мощности всех сегментов.

Вероятность покрытия для истинного PSD в виде скаляра в области значений (0,1). Выход, pxxc, содержит нижние и верхние границы probability × 100%-я оценка интервала для истинного PSD.

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

свернуть все

Оценка PSD, возвращенная как неотрицательный вектор-столбец с действительным знаком или матрица. Каждый столбец pxx оценка PSD соответствующего столбца x. Модули оценки PSD находятся в единицах величины в квадрате данных временных рядов на модульную частоту. Например, если входные данные находятся в вольтах, оценка PSD находится в модулях вольт в квадрате на модульную частоту. Какое-то время ряд в вольтах, если вы принимаете сопротивление 1 Ω и задаете частоту дискретизации в герц, оценка PSD, находится в ваттах на герц.

Типы данных: single | double

Нормированные частоты, возвращенные как вектор-столбец с действительным знаком. Если pxx односторонняя оценка PSD, w охватывает интервал [0, π] если nfft является четным и [0, π), если nfft является нечетным. Если pxx двухсторонняя оценка PSD, w охватывает интервал [0,2π). Для сосредоточенной DC оценки PSD, w охватывает интервал (–π, π] для даже nfft и (–π, π) для нечетного nfft.

Типы данных: double

Циклические частоты, возвращенные как вектор-столбец с действительным знаком. Для односторонней оценки PSD, f охватывает интервал [0, fs/2], когда nfft является четным и [0, fs/2) когда nfft является нечетным. Для двухсторонней оценки PSD, f охватывает интервал [0, fs). Для сосредоточенной DC оценки PSD, f охватывает интервал (–fs/2, fs/2] циклы/единица времени для даже длины nfft и (–fs/2, fs/2) циклы/единица времени для нечетной длины nfft.

Типы данных: double | single

Доверительные границы, возвращенные как матрица с элементами с действительным знаком. Размер строки матрицы равен длине оценки PSD, pxx. pxxc имеет вдвое больше столбцов как pxx. Нечетные столбцы содержат нижние границы доверительных интервалов, и четные столбцы содержат верхние границы. Таким образом, pxxc(m,2*n-1) более низкая доверительная граница и pxxc(m,2*n) верхняя доверительная граница, соответствующая оценке pxx(m,n). Вероятность покрытия доверительных интервалов определяется значением probability входной параметр.

Типы данных: single | double

Больше о

свернуть все

Перекрытый сегмент валлийцев, составляющий в среднем спектральную оценку

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

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

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

Ссылки

[1] Hayes, Монсон Х. Статистическая цифровая обработка сигналов и моделирование. Нью-Йорк: John Wiley & Sons, 1996.

[2] Stoica, Петр и Рэндольф Моисей. Спектральный анализ сигналов. Верхний Сэддл-Ривер, NJ: Prentice Hall, 2005.

Расширенные возможности

Генерация кода C/C++
Генерация кода C и C++ с помощью MATLAB® Coder™.

Представлено до R2006a