Получить оценку Welch PSD входного сигнала, состоящего из дискретно-временной синусоиды с угловой частотой $$\mathrm{\delta /4}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом.

Создайте синусоидальную волну с угловой частотой $$\mathrm{\delta /4}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом. Сбросьте генератор случайных чисел для воспроизводимых результатов. Сигнал имеет длину $${}_{}Nx\mathrm{}\mathrm{=}\mathrm{}$$320 выборок.

rng default

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

Получите оценку Welch PSD, используя окно Хэмминга по умолчанию и длину DFT. Длина сегмента по умолчанию равна 71 выборке, а длина ДПФ равна 256 точкам, что дает разрешение по частоте $$\mathrm{}\mathrm{}\mathrm{}\mathrm{2\delta /256}$$ рад/выборка. Поскольку сигнал имеет действительное значение, периодограмма является односторонней и имеются точки 256/2 + 1. Постройте график оценки PSD Welch.

pxx = pwelch(x);

pwelch(x)

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

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

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% перекрытием между секциями. Укажите ту же длину БПФ, что и на предыдущем шаге. Вычислите оценку Welch 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

Получить оценку Welch PSD входного сигнала, состоящего из дискретно-временной синусоиды с угловой частотой $$\mathrm{\delta /3}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом.

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

rng default

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

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

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

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

Получить оценку Welch PSD входного сигнала, состоящего из дискретно-временной синусоиды с угловой частотой $$\mathrm{\delta /4}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом.

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

rng default

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

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

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

plot(10*log10(pxx))

Получить оценку Welch PSD входного сигнала, состоящего из дискретно-временной синусоиды с угловой частотой $$\mathrm{\delta /4}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом.

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

rng default

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

Получить оценку Welch PSD, делящую сигнал на сегменты длиной 100 выборок. Используйте перекрытие по умолчанию 50%. Задайте длину DFT, равную 640 точкам, так чтобы частота $$\mathrm{\delta /4}$$ рад/выборка соответствовала ячейке DFT (ячейка 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 для перекрывающегося сегмента Welch предыдущего сигнала. Используйте длину сегмента 500 выборок с 300 перекрывающимися выборками. Используйте 500 точек DFT, чтобы частота 100 Гц падала непосредственно на ячейку DFT. Введите частоту дискретизации для вывода вектора частот в Гц. Постройте график результата.

[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 Welch и спектры максимального и минимального удержания сигнала. Постройте график результатов.

[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')

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

Создайте сигнал, состоящий из наложения 100 Гц и 150 Гц синусоидальных волн в аддитивном белом 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 образцов). Постройте график оценки PSD WOSA вместе с доверительным интервалом и увеличьте масштаб интересующей области частот около 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\mathrm{(}\mathrm{}\mathrm{}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));

Получение спектра мощности с центром постоянного тока методом Уэлча. Используйте длину сегмента 500 выборок с 300 перекрывающимися выборками и длиной DFT 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 обусловлено действием аддитивного шума.

Генерируют 1024 выборки многоканального сигнала, состоящего из трех синусоид в аддитивном $$N\mathrm{(}\mathrm{}0,1$$) белом гауссовом шуме. Частоты синусоид составляют $$\mathrm{}$$δ/2, $$\mathrm{}$$λ/3 и $$\mathrm{}$$λ/4 рад/образец. Оцените PSD сигнала с помощью метода Уэлча и постройте его график.

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

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

pwelch(x)