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