Оценка спектра с использованием комплексных данных - тест Марпла

Этот пример показывает, как выполнить спектральную оценку на данных временных рядов. Мы используем тест Марпла (Комплексные данные в L. Marple: S.L. Marple, Jr, Digital Spectral Analysis with Applications, Prentice Hall, Englewood Cliffs, NJ 1987.)

Тестовые данные

Начнем с загрузки тестовых данных:

load marple

Большинство стандартных программ в System Identification Toolbox™ поддерживают комплексные данные. Для графического изображения мы исследуем реальную и мнимую части данных отдельно, однако.

Сначала рассмотрим данные:

subplot(211),plot(real(marple)),title('Real part of data.')
subplot(212),plot(imag(marple)),title('Imaginary part of data.')

Figure contains 2 axes. Axes 1 with title Real part of data. contains an object of type line. Axes 2 with title Imaginary part of data. contains an object of type line.

В качестве предварительного шага анализа давайте проверим периодограмму данных:

per = etfe(marple);
w = per.Frequency;
clf
h = spectrumplot(per,w);
opt = getoptions(h);
opt.FreqScale = 'linear';
opt.FreqUnits = 'Hz';
setoptions(h,opt)

Figure contains an axes. The axes with title From: e@y1 To: y1 contains an object of type line. This object represents per.

Поскольку запись данных находится только 64 отсчетах, а периодограмма вычисляется для 128 частот, мы четко видим колебания из узкого частотного окна. Поэтому мы применяем некоторое сглаживание к периодограмме (соответствующее частотному разрешению 1/32 Гц):

sp = etfe(marple,32);
spectrumplot(per,sp,w);

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent per, sp.

Давайте теперь попробуем Blackman-Tukey подход к оценке спектра:

ssm = spa(marple); % Function spa performs spectral estimation
spectrumplot(sp,'b',ssm,'g',w,opt);    
legend({'Smoothed periodogram','Blackman-Tukey estimate'});

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Smoothed periodogram, Blackman-Tukey estimate.

Длина окна по умолчанию дает очень узкое окно задержки для этого небольшого объема данных. Мы можем выбрать большую задержку окна:

ss20 = spa(marple,20);
spectrumplot(sp,'b',ss20,'g',w,opt);
legend({'Smoothed periodogram','Blackman-Tukey estimate'});

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Smoothed periodogram, Blackman-Tukey estimate.

Оценка авторегрессивной (AR) модели

Параметрическая AR-модель 5 порядка вычисляется:

t5 = ar(marple,5);

Сравните с оценкой периодограммы:

spectrumplot(sp,'b',t5,'g',w,opt); 
legend({'Smoothed periodogram','5th order AR estimate'});

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Smoothed periodogram, 5th order AR estimate.

AR-команда фактически охватывает 20 различных методов оценки спектра. Вышеуказанное было то, что известно как «измененная ковариационная оценка» в книге Марпла.

Некоторые другие хорошо известные из них получают с:

tb5 = ar(marple,5,'burg');      % Burg's method
ty5 = ar(marple,5,'yw');        % The Yule-Walker method
spectrumplot(t5,tb5,ty5,w,opt);
legend({'Modified covariance','Burg','Yule-Walker'})

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 3 objects of type line. These objects represent Modified covariance, Burg, Yule-Walker.

Оценка модели AR с помощью инструментального переменного подхода

AR-моделирование также может быть сделано с помощью подхода Instrumental Variable. Для этого используем функцию ivar:

ti = ivar(marple,4); 
spectrumplot(t5,ti,w,opt);
legend({'Modified covariance','Instrumental Variable'})

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Modified covariance, Instrumental Variable.

Авторегрессионно-скользящее среднее значение (ARMA) модель спектров

Кроме того, System Identification Toolbox охватывает ARMA-моделирование спектров:

ta44 = armax(marple,[4 4]); % 4 AR-parameters and 4 MA-parameters
spectrumplot(t5,ta44,w,opt);
legend({'Modified covariance','ARMA'})

Figure contains an axes. The axes with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Modified covariance, ARMA.