基于MATLAB的谱估计实现 或简写为
RX(m)?E{X[n]X[n?m]} (2-1)
当RX(m)满足
m??????RX(m)??的条件时，我们定义X[n]的功率谱密度为RX(m)的
离散傅里叶变换，并记为GX(?)
GX(?)?m????R??X(m)e?jm?T
式中，T是随机序列相邻各值的时间间隔。 2.2.2 随机过程功率谱密度的性质
1，功率谱密度为非负的。 2，功率谱密度是?的实函数。
3，对于实随机过程来说，功率谱密度是?的偶函数。 4，功率谱密度可积。
2.2.3 功率谱密度与自相关函数的关系
可以证明，平稳随机过程的自相关函数与功率谱密度之间构成傅里叶变换对，即
GX(?)??以及
RX(?)?12?????RX(?)e?j??d?
?????GX(?)e?j??d?
这一关系就是著名的维纳——辛钦定理。它给出了平稳随机过程的时域特性与频域特性之间的联系，它是分析随机信号的一个最重要，最基本的公式。
2.3 估计量的评价标准
2.3.1 无偏性
设X1,X2,…,Xn是总体X的一个样本。???是包含在总体X的分布中的待估参数，这里?是?的取值范围。
无偏性的描述如下
4
基于MATLAB的谱估计实现 ????(X,X,…,X)的数学期望E[??]存在，且对于任意???有 若估计量?12n?]?? E[?则称??是?的无偏估计量。
可以证明，样本均值（也就是总体的一阶矩）是总体均值的无偏估计量，样本方差是总体方差的无偏估计量。 2.3.2 有效性
????(X,X,…,X)与?????(X,X,…,X)都是?的无偏估计量，设?若对于任意1112n2212n???，有
?]?D[??] D[?12?有效。 且至少对于某一个???上式中的不等号成立，则称??1较?22.3.3 相合性
无偏性和有效性都是在样本容量n固定的前提下提出的。我们自然希望随着样本容量的增大，一个估计量的值稳定于待估参数的真值。这样，对估计量又有下述相合性的要求。
????(X,X,…,X)为参数?的估计量，若对于任意???，当n??时设?12n?(X,X,…,X)依概率收敛于?，则称??是?的相合估计量。 ?12n即，若对于任意???都满足：对于任意??0，有
?????}?1 limP{?n??则称??是?的相合估计量。
可以证明，样本k(k?1)阶矩是总体X的k阶矩的相合估计量。
相合性是对一个估计量的基本要求，若估计量不具有相合性，那么不论将样本容量n取得多么大，都不能将?估计得足够准确，这样的估计量是不可取的。
3 经典谱估计及其仿真
经典谱估计方法分为两种：BT法和周期图法。下面逐一阐述。
3.1 周期图
3.1.1 周期图的定义
对于一个离散随机过程，由公式(2-1)，以及样本的一阶矩就是对总体期望的无
5
基于MATLAB的谱估计实现 偏估计可得：
?[m]?1 RXNN?m?1?n?0X[n]X[n?m] (3-1)
?[m]是对此离散时间随机过程自相关函数的估计量。 其中N是信号的采样点数，RX定义矩形窗函数如下
?1,0?n?N?1 d[n]??0,其他?则式(3-1)可写成
?[m]?1RXNn????d[n]X[n]d[n?m]X[n?m]
??令d[n]X[n]?Y[n]，上式可简写为
?[m]?1 RXNn????Y[n]Y[n?m]???1Y[m]?Y[?m] N?[m]可看成是Y[m]与Y[?m]的卷积。而Y[m]的傅里叶变换 即RXm????Y[m]e???j?m?m????d[m]X[m]e???j?m??X[m]e?j?m?XN(ej?)
m?0N?1XN(ej?)是有限长序列X[n]的傅里叶变换，显然XN(ej?)是周期性的。当X[n]为实序
?列时，Y[n]也为实序列，根据傅里叶变换的性质，Y[?m]的傅里叶变换为XN(ej?)。
根据功率谱密度与自相关函数的关系，功率谱估计为自相关函数估计的傅里叶变换，即
?(?)?1X(ej?)X?(ej?)?1X(ej?)2 GXNNNNN直接将XN(ej?)模的平方除以N求得的功率谱估计的方法称为周期图法，其结果用IN(?)表示，即
1j?2?XN(e) (3-2) GX(?)?IN(?)?N周期图法是利用数据的傅里叶变换直接求得的，而不再计算自相关函数，所以又称直接法。由于序列的傅里叶变化可利用FFT计算而提高运算效率，这是周期图法
6
基于MATLAB的谱估计实现 的主要优点。 3.1.2 周期图的性能
为了了解周期图法的谱估计效果，我们来讨论它的估计均值和方差。
E[IN(?)]?E[211?XN(ej?)]?E[XN(ej?)XN(ej?)] NNN?11N?1j?n?E[?X(n)e?X(k)e?j?k] Nn?0k?0????1j?n ?E[?d[n]X[n]e?d[k]X[k]e?j?k]
Nn???k???????1 ?E[??d[n]d[k]X[n]X[k]e?j?(k?n)]
Nn???k???令m?k?n代入上式得
1????E[IN(?)]????d[n]d[n?m]E[X[n]X[n?m]]e?j?m?
m???N?n???????式中
m????b[m]RN???j?m[m]e （3-3） X1bN[m]?Nn????d[n]d[n?m]
???m,m?N?1?1? ?? N?0,其他?由于bN[m]为两个矩形窗函数的卷积，因此它是一个三角窗函数，其傅里叶变换为
?N?2sin??1?1??j?n?2B(?)?d[n]e? N???N?N?sin??n?????2
7
??? ??2基于MATLAB的谱估计实现 式(3-3)为两个函数乘积的傅里叶变换，等于他们各自傅里叶变换的卷积，即 E[IN(?)]?BN(?)?GX(?)
由上式可知，除非BN(?)是一个冲激函数，否则E[IN(?)]将不等于GX(?)，故周期图作为功率谱的估计是有偏差的。
当N??时
?m? limbN[m]?lim?1???1
N??N??N??BN?(?)??( )故 limN??EI[)G]X? ( (3-4) 此时 limN?(?N??因此，周期图作为功率谱估计，当N??时是无偏的，即渐近无偏。 下面讨论谱估计的方差。
2为分析简单起见，通常假设X[n]零均值，方差为?X的实高斯白噪声序列，即功2率谱密度为常数?X。按方差的定义应有
2 var[IN(?)]?E[IN(?)]?E2[IN(?)]
2由于假设X[n]零均值，方差为?X的实高斯白噪声序列，根据式(3-4)可知，当2N??时，E[IN(?)]?GX(?)??X。
2为了求E[IN先求IN(?)在两个频率?1和?2处的协方差，最后令???1??2。 (?)]，
根据式(3-2)可得
22??1E[IN(?1)IN(?2)]?E?2XN(ej?1)XN(ej?2)?
?N? ?1R[n]RN[k]RN[p]RN[q] 2????NNnkpq ?E[X[n]X[k]X[p]X[q]]e?j?1(n?k)e?j?2(p?q) (3-5)
其中RN[?]为矩形窗函数。利用正态白噪声、多元正态随机变量的多阶矩公式，有
E[X[n]X[k]X[p]X[q]]?E[X[n]X[k]]E[X[p]X[q]]?E[X[n]X[p]]E[X[k]X[q]]
8