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

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

rng default

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

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

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

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 Уэлча входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой $$\pi /3$$ рад/образец с добавкой $$N(0,1)$$ белый шум.

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

rng default

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

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

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

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

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

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

rng default

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

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

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

plot(10*log10(pxx))

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

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

rng default

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

Получите оценку PSD Welch, разделяющую сигнал на сегменты 100 выборки в длину. Используйте перекрытие по умолчанию 50%. Задайте длину ДПФ, чтобы она составляла 640 точек, так чтобы частота $$\pi /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 перекрытого сегмента Welch предыдущего сигнала. Используйте длину сегмента 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 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 перекрытого сегмента (WOSA) Welch. Хотя не является необходимым условием для статистической значимости, частоты в оценке Уэлча, где нижняя доверительная граница превышает верхнюю доверительную границу для окружающих оценок 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(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 перекрываемыми выборками и длину ДПФ 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(0,1)$$ белый Гауссов шум. Частоты синусоидов $$\pi /2$$, $$\pi /3$$, и $$\pi /4$$ рад/образец. Оцените PSD сигнала с помощью метода Уэлча и постройте график.

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

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

pwelch(x)