Этот пример показывает, как уменьшить смещение и изменчивость в периодограмме. Использование окна может уменьшить смещение в периодограмме, а использование окон с усреднением может уменьшить изменчивость.
Используйте широкополосные стационарные авторегрессивные (AR) процессы, чтобы показать эффекты смещения и изменчивости в периодограмме. Процессы AR представляют удобную модель, потому что их PSD имеют выражение закрытой формы. Создайте AR (2) модель следующей формы:
где является нулевой средней последовательностью белого шума с некоторой заданным отклонением. В этом примере примите отклонение и период дискретизации равными 1. Чтобы симулировать предыдущий процесс AR (2), создайте полнополюсный (БИХ) фильтр. Просмотрите величину ответ фильтра.
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')
Обратите внимание, что оценка периодограммы теперь следует истинному 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')
Метод многозначности дает оценку PSD со значительно меньшей изменчивостью, чем периодограмма. Поскольку метод multitaper также использует окна, вы видите, что смещение периодограммы также решено.
