복소수 데이터를 사용한 스펙트럼 추정 - Marple의 테스트 케이스

이 예제에서는 시계열 데이터에서 스펙트럼 추정을 수행하는 방법을 보여줍니다. Marple의 테스트 케이스(The complex data in 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.')

예비 분석 단계로, 데이터 주기도를 검사해 보겠습니다.

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

데이터 레코드는 64개 샘플뿐이고 주기도는 128개 주파수에 대해 계산되었기 때문에 좁은 주파수 윈도우에서 진동을 명확하게 확인할 수 있습니다. 따라서 주기도에 평활화를 적용합니다(1/32Hz의 주파수 분해능에 해당).

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

이제 블랙맨-튜키(Blackman-Tukey) 방법을 스펙트럼 추정에 사용해 보겠습니다.

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

디폴트 윈도우 길이는 이 소량의 데이터에 대해 매우 좁은 지연 윈도우를 제공합니다. 다음 방법으로 더 큰 지연 윈도우를 선택할 수 있습니다.

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

자기회귀(AR) 모델 추정하기

모수적 5차 AR 모델은 다음과 같이 계산됩니다.

t5 = ar(marple,5);

주기도 추정값과 비교합니다.

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

사실상 AR 명령은 스펙트럼 추정을 위한 20개의 서로 다른 방법을 포함합니다. 위에 나온 방법은 Marple 교재에서 '수정된 공분산 추정값(the modified covariance estimate)'으로 알려져 있습니다.

잘 알려진 몇 가지 다른 방법의 결과는 다음과 같이 얻습니다.

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'})

도구 변수 방법을 사용하여 AR 모델 추정하기

AR 모델링은 도구 변수 방법을 사용하여 수행할 수도 있습니다. 이를 위해 여기서는 ivar 함수를 사용합니다.

ti = ivar(marple,4); 
spectrumplot(t5,ti,w,opt);
legend({'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'})