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

Создайте синусоиду с угловой частотой $$\pi /4$$ рад/выборка с дополнением $$N(0,1)$$ белый шум. Сигнал является 320 выборками в длине. Получите периодограмму с помощью прямоугольного окна по умолчанию и длины ДПФ. Длина ДПФ является следующей степенью двойки, больше, чем длина сигнала или 512 точек. Поскольку сигнал с действительным знаком и имеет даже длину, периодограмма является односторонней и существуют точки 512/2+1.

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
[pxx,w] = periodogram(x);
plot(w,10*log10(pxx))

Повторите график с помощью periodogram без выходных параметров.

periodogram(x)

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

Создайте синусоиду с угловой частотой $$\pi /4$$ радианы/выборка с дополнением $$N(0,1)$$ белый шум. Сигнал является 320 выборками в длине. Получите измененную периодограмму с помощью Окна Хэмминга и длины ДПФ по умолчанию. Длина ДПФ является следующей степенью двойки, больше, чем длина сигнала или 512 точек. Поскольку сигнал с действительным знаком и имеет даже длину, периодограмма является односторонней и существуют точки 512/2+1.

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

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

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

n = 0:319;
x = cos(pi/4*n)+randn(size(n));
nfft = length(x);
periodogram(x,[],nfft)

Получите периодограмму Волка (относительное солнечное пятно) данные о номере, выбираемые ежегодно между 1 700 и 1987.

Загрузите относительные данные о номере солнечного пятна. Получите периодограмму с помощью прямоугольного окна по умолчанию и количества точек ДПФ (512 в этом примере). Частота дискретизации для этих данных является 1 выборкой/год. Постройте периодограмму.

load sunspot.dat
relNums=sunspot(:,2);

[pxx,f] = periodogram(relNums,[],[],1);

plot(f,10*log10(pxx))
xlabel('Cycles/Year')
ylabel('dB / (Cycles/Year)')
title('Periodogram of Relative Sunspot Number Data')

Вы видите в предыдущей фигуре, что существует пик в периодограмме приблизительно в 0,1 циклах/год, которая указывает на период приблизительно 10 лет.

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

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

[pxx,w] = periodogram(x,[],[pi/4 pi/2]);
pxx

pxx = 1×2

   14.0589    2.8872

[pxx1,w1] = periodogram(x);
plot(w1/pi,pxx1,w/pi,2*pxx,'o')
legend('pxx1','2 * pxx')
xlabel('\omega / \pi')

Полученные значения периодограммы являются 1/2 значения в односторонней периодограмме. Когда вы оцениваете периодограмму в определенном наборе частот, вывод является двухсторонней оценкой.

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

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

freq = [100 200];
pxx = periodogram(x,[],freq,fs)

pxx = 1×2

    0.2647    0.2313

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

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

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

Получите периодограмму оценка PSD с 95%-доверительными-границами. Постройте периодограмму наряду с доверительным интервалом и увеличьте масштаб видимой области частоты около 100 и 150 Гц.

[pxx,f,pxxc] = periodogram(x,rectwin(length(x)),length(x),fs,...
    'ConfidenceLevel',0.95);

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

xlim([85 175])
xlabel('Hz')
ylabel('dB/Hz')
title('Periodogram with 95%-Confidence Bounds')

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

Получите периодограмму синусоиды на 100 Гц в дополнении $$N(0,1)$$ шум. Данные выбираются на уровне 1 кГц. Используйте опцию 'centered', чтобы получить сосредоточенную DC периодограмму и построить результат.

fs = 1000;
t = 0:0.001:1-0.001;
x = cos(2*pi*100*t)+randn(size(t));
periodogram(x,[],length(x),fs,'centered')

Сгенерируйте сигнал, который состоит из синусоиды на 200 Гц, встроенной в белый Гауссов шум. Сигнал выбирается на уровне 1 кГц в течение 1 секунды. Шум имеет отклонение 0,01 ². Сбросьте генератор случайных чисел для восстанавливаемых результатов.

rng('default')

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
N = length(t);
x = sin(2*pi*t*200)+0.01*randn(size(t));

Используйте БПФ, чтобы вычислить спектр мощности сигнала, нормированного длиной сигнала. Синусоида в интервале, таким образом, вся степень сконцентрирована в одной выборке частоты. Постройте односторонний спектр. Увеличьте масштаб к близости пика.

q = fft(x,N);
ff = 0:Fs/N:Fs-Fs/N;

ffts = q*q'/N^2

ffts = 0.4997

ff = ff(1:floor(N/2)+1);
q = q(1:floor(N/2)+1);

stem(ff,abs(q)/N,'*')
axis([190 210 0 0.55])

Используйте periodogram, чтобы вычислить спектр мощности сигнала. Задайте окно Hann и длину БПФ 1 024. Найдите процентную разницу между предполагаемой степенью на уровне 200 Гц и фактическим значением.

wind = hann(N);

[pun,fr] = periodogram(x,wind,1024,Fs,'power');

hold on
stem(fr,pun)

periodogErr = abs(max(pun)-ffts)/ffts*100

periodogErr = 4.7349

Повторно вычислите спектр мощности, но на этот раз используйте переназначение. Постройте новую оценку и сравните ее максимум со значением БПФ.

[pre,ft,pxx,fx] = periodogram(x,wind,1024,Fs,'power','reassigned');

stem(fx,pre)
hold off
legend('Original','Periodogram','Reassigned')

reassignErr = abs(max(pre)-ffts)/ffts*100

reassignErr = 0.0779

Оцените степень синусоиды на определенной частоте с помощью опции 'power'.

Создайте синусоиду на 100 Гц одна секунда в длительности, выбранной на уровне 1 кГц. Амплитуда синусоиды 1.8, который приравнивается к степени 1,8 ²/2 = 1.62. Оцените степень с помощью опции 'power'.

fs = 1000;
t = 0:1/fs:1-1/fs;
x = 1.8*cos(2*pi*100*t);
[pxx,f] = periodogram(x,hamming(length(x)),length(x),fs,'power');
[pwrest,idx] = max(pxx);
fprintf('The maximum power occurs at %3.1f Hz\n',f(idx))

The maximum power occurs at 100.0 Hz

fprintf('The power estimate is %2.2f\n',pwrest)

The power estimate is 1.62

Сгенерируйте 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);

periodogram(x)

Создайте функциональный periodogram_data.m, который возвращает измененную периодограмму оценка PSD с помощью окна, window. Длина окна и числа точек в ДПФ совпадает с длиной inputData.

type periodogram_data

function [pxx,f] = periodogram_data(inputData,window)
%#codegen
nfft = length(inputData);
[pxx,f] = periodogram(inputData,window,nfft);
end

Директива %#codegen указывает, что код MATLAB предназначается для генерации кода. Используйте codegen, чтобы сгенерировать файл MEX periodogram_mex. Опция -args задает входные спецификации для MEX - файл. input_data 1024x1 случайный вектор. Задайте Окно Хэмминга длины 1024. Чтобы просмотреть отчет генерации кода, можно добавить флаг -report в конце следующей строки кода.

codegen periodogram_data -args {randn(1024,1),hamming(1024)} -o periodogram_mex

Постройте периодограмму шумной синусоиды с помощью и periodogram_mex и periodogram_data.m.

t = linspace(0,1,1024);
x = 2*cos(2*pi*50*t')+randn(size(t'));

[pxMex,fMex] = periodogram_mex(x,hann(length(x)));

[pxMat,fMat] = periodogram_data(x,hamming(length(x)));
plot(fMex,10*log10(pxMex))
hold on
plot(fMat,10*log10(pxMat),'--')
hold off
grid on

title('Periodogram Power Spectral Density Estimate')
xlabel('Normalized Frequency (\times\pi rad/sample)')
ylabel('Power/frequeny (dB/(rad/sample))')

legend('Periodogram with MEX function','Periodogram with MATLAB function')