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

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

rng default

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

Получите валлийскую оценку PSD с помощью Окна Хэмминга по умолчанию и длины ДПФ. Длина сегмента по умолчанию является 71 выборкой, и длина ДПФ является 256 точками, приводящими к разрешению частоты $$2\pi /\mathrm{}\mathrm{}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%-м перекрытием между разделами. Задайте ту же длину БПФ как на предыдущем шаге. Вычислите кратковременное преобразование Фурье и проверьте, что оно дает тот же результат как предыдущие две процедуры.

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, делящую сигнал на сегменты 132 выборки в длине. Сегменты сигнала умножаются на Окно Хэмминга 132 выборки в длине. Количество перекрытых выборок не задано, таким образом, оно установлено в 132/2 = 66. Длина ДПФ является 256 точками, приводя к разрешению частоты $$2\pi /\mathrm{}\mathrm{}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, делящую сигнал на сегменты 100 выборок в длине. Сегменты сигнала умножаются на Окно Хэмминга 100 выборок в длине. Количество перекрытых выборок равняется 25. Длина ДПФ является 256 точками, приводящими к разрешению частоты $$2\pi /\mathrm{}\mathrm{}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, делящую сигнал на сегменты 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 предыдущего сигнала. Используйте продолжительность сегмента 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)$$ белый Гауссов шум. Частоты синусоид $$\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)