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

Создайте синусоидальную волну с угловой частотой $$\mathrm{\delta /4}$$ рад/образец с аддитивным $$N\mathrm{(}\mathrm{}0,1$$) белым шумом. Сигнал имеет длину 320 выборок. Получите периодограмму, используя прямоугольное окно по умолчанию и длину DFT. Длина DFT - это следующая степень, которая на два больше, чем длина сигнала, или 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)

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

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

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

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

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

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

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

Загрузите данные относительного числа солнечных пятен. Получите периодограмму, используя прямоугольное окно по умолчанию и количество точек DFT (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 лет.

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

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));

Получить оценку ИПУ периодограммы с границами достоверности 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\mathrm{(}\mathrm{}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\mathrm{(}\mathrm{}0,1$$) белом гауссовом шуме. Частоты синусоид составляют $$\mathrm{}$$δ/2, $$\mathrm{}$$λ/3 и $$\mathrm{}$$λ/4 рад/образец. Оцените ИПУ сигнала по периодограмме и постройте его график.

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) для генерации функции MEX.

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

Code generation successful.

Вычислите оценку PSD для 1024-образной шумной синусоиды, используя функцию периодограммы и созданную функцию MEX. Задайте нормализованную частоту синусоиды $\mathrm{}\mathrm{2\delta /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')