利用复数数据进行频谱估计 - 马普尔的测试案例

打开实时脚本

本示例演示了如何对时间序列数据进行谱估计。我们采用马普尔的测试案例 (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(2,1,1)
plot(real(marple))
title('Real part of data.')
subplot(2,1,2)
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.

作为初步分析步骤，让我们检查一下数据的周期图：

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, ylabel Power (dB) contains an object of type line. This object represents per.

由于数据记录仅包含 64 个采样点，而频谱图是针对 128 个频率计算的，因此我们可以清楚地看到来自窄频带的振荡。因此，我们对周期图进行了一些平滑处理（对应的频率分辨率为 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, ylabel Power (dB) contains 2 objects of type line. These objects represent per, sp.

现在，让我们尝试使用布莱克曼-图基法进行频谱估计：

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

自回归 (AR) 模型的估计

通过以下方式计算出一个参数化五阶 AR 模型：

t5 = ar(marple,5);

与频谱图估计值进行比较：

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, ylabel Power (dB) contains 2 objects of type line. These objects represent Smoothed periodogram, 5th order AR estimate.

实际上，AR-command 涵盖了 20 种不同的频谱估计方法。上述内容就是马普尔 (Marple) 的书中所谓的“修正协方差估计”。

其他一些众所周知的示例可通过以下方式获得：

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, ylabel Power (dB) contains 3 objects of type line. These objects represent Modified covariance, Burg, Yule-Walker.

使用工具变量法估计 AR 模型

AR 建模也可以采用工具变量法来实现。为此，我们使用函数 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, ylabel Power (dB) 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 object. The axes object with title From: e@y1 To: y1, ylabel Power (dB) contains 2 objects of type line. These objects represent Modified covariance, ARMA.