Смещение и изменчивость в периодограмме

Этот пример показывает, как уменьшить смещение и изменчивость в периодограмме. Использование окна может уменьшить смещение в периодограмме, а использование окон с усреднением может уменьшить изменчивость.

Используйте широкополосные стационарные авторегрессивные (AR) процессы, чтобы показать эффекты смещения и изменчивости в периодограмме. Процессы AR представляют удобную модель, потому что их PSD имеют выражение закрытой формы. Создайте AR (2) модель следующей формы:

y(n)-0.75y(n-1)+0.5y(n-2)=ε(n),

где ε(n) является нулевой средней последовательностью белого шума с некоторой заданным отклонением. В этом примере примите отклонение и период дискретизации равными 1. Чтобы симулировать предыдущий процесс AR (2), создайте полнополюсный (БИХ) фильтр. Просмотрите величину ответ фильтра.

B2 = 1;
A2 = [1 -0.75 0.5];
fvtool(B2,A2)

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains an object of type line.

Этот процесс является полосно-пропускающим. Динамическая область значений PSD составляет приблизительно 14,5 дБ, что можно определить с помощью следующего кода.

[H2,W2] = freqz(B2,A2,1e3,1);
dr2 = max(20*log10(abs(H2)))-min(20*log10(abs(H2)))
dr2 = 14.4984

Исследуя размещение полюсов, вы видите, что этот процесс AR (2) является стабильным. Два полюса находятся внутри модуля круга.

fvtool(B2,A2,'Analysis','polezero')

Figure Filter Visualization Tool - Pole-Zero Plot contains an axes and other objects of type uitoolbar, uimenu. The axes with title Pole-Zero Plot contains 4 objects of type line, text.

Затем создайте процесс AR (4), описанный следующим уравнением:

y(n)-2.7607y(n-1)+3.8106y(n-2)-2.6535y(n-3)+0.9238y(n-4)=ε(n).

Используйте следующий код для просмотра величины ответа этой БИХ-системы.

B4 = 1;
A4 = [1 -2.7607 3.8106 -2.6535 0.9238];
fvtool(B4,A4)

Figure Filter Visualization Tool - Magnitude Response (dB) contains an axes and other objects of type uitoolbar, uimenu. The axes with title Magnitude Response (dB) contains an object of type line.

Исследуя размещение полюсов, можно увидеть, что этот процесс AR (4) также стабильен. Четыре полюса находятся внутри модуля круга.

fvtool(B4,A4,'Analysis','polezero')

Figure Filter Visualization Tool - Pole-Zero Plot contains an axes and other objects of type uitoolbar, uimenu. The axes with title Pole-Zero Plot contains 4 objects of type line, text.

Динамическая область значений этого PSD приблизительно 65 дБ, намного больше, чем модель AR (2).

[H4,W4] = freqz(B4,A4,1e3,1);
dr4 = max(20*log10(abs(H4)))-min(20*log10(abs(H4)))
dr4 = 64.6213

Чтобы симулировать реализации из этих процессов AR (p), используйте randn и filter. Установите генератор случайных чисел в настройки по умолчанию, чтобы получить повторяемые результаты. Постройте график реализаций.

rng default
x = randn(1e3,1);
y2 = filter(B2,A2,x);
y4 = filter(B4,A4,x);

subplot(2,1,1)
plot(y2)
title('AR(2) Process')
xlabel('Time')

subplot(2,1,2)
plot(y4)
title('AR(4) Process')
xlabel('Time')

Figure contains 2 axes. Axes 1 with title AR(2) Process contains an object of type line. Axes 2 with title AR(4) Process contains an object of type line.

Вычислите и постройте график периодограмм реализаций AR (2) и AR (4). Сравните результаты с истинным PSD. Обратите внимание, что periodogram преобразует частоты в миллигерц для графического изображения.

Fs = 1;
NFFT = length(y2);

subplot(2,1,1)
periodogram(y2,rectwin(NFFT),NFFT,Fs)
hold on
plot(1000*W2,20*log10(abs(H2)),'r','linewidth',2)
title('AR(2) PSD and Periodogram')

subplot(2,1,2)
periodogram(y4,rectwin(NFFT),NFFT,Fs)
hold on
plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2)
title('AR(4) PSD and Periodogram')
text(350,20,'\downarrow Bias')

Figure contains 2 axes. Axes 1 with title AR(2) PSD and Periodogram contains 2 objects of type line. Axes 2 with title AR(4) PSD and Periodogram contains 3 objects of type line, text.

В случае процесса AR (2) оценка периодограммы соответствует форме истинного PSD, но показывает значительную изменчивость. Это связано с низкими степенями свободы. Выраженные отрицательные отклонения (в дБ) в периодограмме объясняются взятием журнала случайной переменной хи-квадрат с двумя степенями свободы.

В случае процесса AR (4) периодограмма повторяет форму истинного PSD на низких частотах, но отклоняется от PSD на высоких частотах. Это эффект свертки с ядром Фехера. Большая динамическая область значений процесса AR (4) по сравнению с процессом AR (2) делает смещение более выраженным.

Уменьшите смещение, продемонстрированное в процессе AR (4), используя конусность, или окно. В этом примере используйте окно Хэмминга, чтобы сужать реализацию AR (4) перед получением периодограммы.

figure
periodogram(y4,hamming(length(y4)),NFFT,Fs)
hold on
plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2)
title('AR(4) PSD and Periodogram with Hamming Window')
legend('Periodogram','AR(4) PSD')

Figure contains an axes. The axes with title AR(4) PSD and Periodogram with Hamming Window contains 2 objects of type line. These objects represent Periodogram, AR(4) PSD.

Обратите внимание, что оценка периодограммы теперь следует истинному PSD AR (4) во всей частотной области значений Найквиста. Оценки пародограммы все еще имеют только две степени свободы, поэтому использование окна не уменьшает изменчивость пародограммы, но решает проблему смещения.

При непараметрической спектральной оценке двумя методами для увеличения степеней свободы и уменьшения изменчивости периодограммы являются среднее с перекрывающимся сегментом Уэлча и спектральная оценка с множеством указателей.

Получите оценку многозначности AR (4) временных рядов с использованием временного продукта половины полосы пропускания 3,5. Постройте график результата.

NW = 3.5;

figure
pmtm(y4,NW,NFFT,Fs)
hold on
plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2)
legend('Multitaper Estimate','AR(4) PSD')

Figure contains an axes. The axes with title Thomson Multitaper Power Spectral Density Estimate contains 2 objects of type line. These objects represent Multitaper Estimate, AR(4) PSD.

Метод многозначности дает оценку PSD со значительно меньшей изменчивостью, чем периодограмма. Поскольку метод multitaper также использует окна, вы видите, что смещение периодограммы также решено.

См. также

|

Для просмотра документации необходимо авторизоваться на сайте