Main Content

Spectrum Estimation Using Complex Data - Marple's Test Case

This example shows how to perform spectral estimation on time series data. We use Marple's test case (The complex data in L. Marple: S.L. Marple, Jr, Digital Spectral Analysis with Applications, Prentice-Hall, Englewood Cliffs, NJ 1987.)

Test Data

Let us begin by loading the test data:

load marple

Most of the routines in System Identification Toolbox™ support complex data. For plotting we examine the real and imaginary parts of the data separately, however.

First, take a look at the data:

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

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

As a preliminary analysis step, let us check the periodogram of the data:

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 object. The axes object with title From: e@y1 To: y1 contains an object of type line. This object represents per.

Since the data record is only 64 samples, and the periodogram is computed for 128 frequencies, we clearly see the oscillations from the narrow frequency window. We therefore apply some smoothing to the periodogram (corresponding to a frequency resolution of 1/32 Hz):

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

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

Let us now try the Blackman-Tukey approach to spectrum estimation:

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 object. The axes object with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Smoothed periodogram, Blackman-Tukey estimate.

The default window length gives a very narrow lag window for this small amount of data. We can choose a larger lag window by:

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

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

Estimating an Autoregressive (AR) Model

A parametric 5-order AR-model is computed by:

t5 = ar(marple,5);

Compare with the periodogram estimate:

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

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

The AR-command in fact covers 20 different methods for spectrum estimation. The above one was what is known as 'the modified covariance estimate' in Marple's book.

Some other well known ones are obtained with:

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 object. The axes object with title From: e@y1 To: y1 contains 3 objects of type line. These objects represent Modified covariance, Burg, Yule-Walker.

Estimating AR Model using Instrumental Variable Approach

AR-modeling can also be done using the Instrumental Variable approach. For this, we use the function ivar:

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

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

Autoregressive-Moving Average (ARMA) Model of the Spectra

Furthermore, System Identification Toolbox covers ARMA-modeling of spectra:

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 object. The axes object with title From: e@y1 To: y1 contains 2 objects of type line. These objects represent Modified covariance, ARMA.