Получите периодограмму входного сигнала, состоящего из синусоиды дискретного времени с угловой частотой $$\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)

Получите периодограмму данных числа Вольфа (относительное пятно солнца), дискретизированных ежегодно между 1700 и 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, четко указывают на значительные колебания во временных рядах.

Создайте сигнал, состоящий из суперпозиции 100 Гц и 150 Гц синусоид аддитивного белого 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' опция для получения периодограммы постоянного тока и построения графика результата.

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 для вычисления степени спектра сигнала. Задайте окно Ханна и длину БПФ 1024. Найдите процентное различие между расчетной степенью в 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

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

periodogram(x)

Создайте функцию periodogram_data.m который возвращает модифицированную оценку периодограммы степени спектральной плотности (PSD) входного сигнала с помощью окна. Функция задает количество дискретных точек преобразования Фурье, равное длине входного сигнала.

type periodogram_data

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

Использование codegen (MATLAB Coder), чтобы сгенерировать MEX-функцию.

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

Code generation successful.

Вычислите оценку PSD шумной синусоиды с 1024 образцом, используя функцию периодограммы и MEX-функцию, которую вы сгенерировали. Задайте синусоидальную нормированную частоту $2\pi /5$ рад/образец и окно Ханна. Постройте график двух оценок, чтобы убедиться, что они совпадают.

N = 1024;
x = 2*cos(2*pi/5*(0:N-1)') + randn(N,1);
periodogram(x,hann(N))
[pxMex,fMex] = periodogram_data(x,hann(N));
hold on
plot(fMex/pi,pow2db(pxMex),':','Color',[0 0.4 0])
hold off
grid on
legend('periodogram','MEX function')