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

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

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

где $$\lambda (n$$) - нулевая средняя последовательность белого шума с некоторой заданной дисперсией. В этом примере предположим, что дисперсия и период выборки равны 1. Для моделирования предыдущего процесса AR (2) создайте всеполюсный (IIR) фильтр. Просмотрите отклик фильтра на величину.

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

Этот процесс является полосовым. Динамический диапазон 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')

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

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

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

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

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

Динамический диапазон этого 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')

Вычислите и постройте график периодограмм реализаций 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')

В случае процесса 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')

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

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

Получение многоабонентской оценки временного ряда 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')

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

См. также

periodogram | pmtm