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(2,1,1)
plot(real(marple))
title('Real part of data.')
subplot(2,1,2)
plot(imag(marple))
title('Imaginary part of data.')

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)

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);

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

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

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

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

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

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