PAM的功率谱 P s ( f ) = σ a 2 T s ∣ G T ( f ) ∣ 2 + m a 2 T s 2 ∑ k = − ∞ + ∞ ∣ G T ( k R s ) ∣ 2 δ ( f − k R s ) P_s(f) = \frac{\sigma_a^2}{T_s}|G_T(f)|^2 + \frac{m_a^2}{T_s^2}\sum_{k=-\infty}^{+\infty}|G_T(kR_s)|^2\delta(f-kR_s) Ps​(f)=Ts​σa2​​∣GT​(f)∣2+Ts2​ma2​​k=−∞∑+∞​∣GT​(kRs​)∣2δ(f−kRs​)

我们首先有这样一个常识，对于一个确定的信号，它唯一对应了一个确定的频谱，这是信号与系统课程中我们学到的看世界的两种角度——时域、频域。

现在我们面对的是随机的信号，不知道信号在下一时刻的取值是1还是0，但是我们知道它出现1的概率，以及它出现0的概率。也就是说，不同时间观察到的信号是不一样的，也就是意味着观察到的频谱是不一样的。现在我们试图在变化中抓住本质，抓住那个不变的东西。

从随机信号分析这门课程里，我们学习到了自相关函数。对于一个广义平稳随机信号而言，它的均值为一个常数，自相关函数只与观察时差有关，而与观察时刻无关。与观察时刻无关就令人十分满意了，因为这意味着就算我明天来观察这个信号，今天所计算出来的自相关函数仍然是有用的。

但是自相关函数仍然不好看，因为它终归是从时域角度观察的。由维纳-辛钦定理我们知道，自相关函数的傅里叶变换就是这个信号的功率谱密度函数，常常简称其为功率谱。因此，我们得出个结论，广义平稳随机信号的功率谱是不变的，是一个能体现随机信号统计特性的统计量，其物理意义是单位频率的平均功率（单位是W/Hz）。

写成数学公式： R X ( τ ) = E [ X ( t + τ ) X ( t ) ] S X ( j ω ) = ∫ − ∞ + ∞ R X ( τ ) e − j ω τ d τ R_X(\tau) = E[X(t+\tau)X(t)]\\ S_X(j\omega)=\int_{-\infty}^{+\infty}R_X(\tau)e^{-j\omega\tau}d\tau RX​(τ)=E[X(t+τ)X(t)]SX​(jω)=∫−∞+∞​RX​(τ)e−jωτdτ

从上述分析我们知道了功率谱是个好东西，但是现实是残酷的，因为我们不可能观察一个信号无穷时长，自然就得不到无穷宽的时差，也就得不到绝对精准的功率谱。

我们现在有的仅仅是有限时长的一段离散信号，这个时候就可以称其为序列了。但是我们试图在仅有的条件下，尽可能地反应现实世界。

在随机信号分析中，我们知道了各态历经这个概念，当观察时间长到足以把随机信号能经历的所有“状态”都经历一遍，那么咱就不观察了。我们就用此时得到的足够长的离散序列来求自相关序列，并通过快速傅里叶变换得到功率谱。

这部分主要涉及到自相关序列的无偏估计和有偏估计，关于其无偏性和有效性的数学推导较为复杂，但学过先导知识里面的课程之后，还是能够顺着参考书里面的步骤推导出来的。这里只列出结果。 u ( 0 ) , u ( 1 ) , . . . , u ( N − 1 ) 是观察到的样本序列。 r u ( m ) 是 自 相 关 序 列 r ^ u ( m ) 是 r u ( m ) 的 有 偏 估 计 ， r ^ ′ u ( m ) 是 r u ( m ) 的 无 偏 估 计 u(0),u(1),...,u(N-1)\text{是观察到的样本序列。} \\r_u(m)是自相关序列\\ \hat r_u(m)是r_u(m)的有偏估计，\hat r{'}_u(m)是r_u(m)的无偏估计 u(0),u(1),...,u(N−1)是观察到的样本序列。ru​(m)是自相关序列r^u​(m)是ru​(m)的有偏估计，r^′u​(m)是ru​(m)的无偏估计

r ^ u ( m ) = 1 N ∑ n = 0 N − 1 − ∣ m ∣ u ( n ) u ∗ ( n + m ) , ∣ m ∣ ≤ N − 1 E [ r ^ u ( m ) ] = N − ∣ m ∣ N r u ( m ) V a r [ r ^ u ( m ) ] = ( N − ∣ m ∣ N ) 2 { ( 1 N − ∣ m ∣ ) 2 ∑ r = − ( N − 1 − ∣ m ∣ ) N − 1 − ∣ m ∣ ( N − ∣ m ∣ − ∣ r ∣ ) [ r u 2 ( r ) + r u ( r − m ) r u ( r + m ) ] } \hat r_u(m) = \frac{1}{N}\sum_{n=0}^{N-1-|m|}u(n)u^*(n+m),|m|\le N-1 \\ E[\hat r_u(m)] = \frac{N-|m|}{N} r_u(m) \\ Var[\hat r_u(m)] = (\frac{N-|m|}{N})^2 \{ (\frac{1}{N-|m|})^2 \sum_{r = -(N-1-|m|)}^{N-1-|m|}(N-|m|-|r|) [r^2_u(r)+r_u(r-m)r_u(r+m)] \} r^u​(m)=N1​n=0∑N−1−∣m∣​u(n)u∗(n+m),∣m∣≤N−1E[r^u​(m)]=NN−∣m∣​ru​(m)Var[r^u​(m)]=(NN−∣m∣​)2{(N−∣m∣1​)2r=−(N−1−∣m∣)∑N−1−∣m∣​(N−∣m∣−∣r∣)[ru2​(r)+ru​(r−m)ru​(r+m)]}

无偏估计 r ^ u ′ ( m ) = 1 N − ∣ m ∣ ∑ N = 0 N − 1 − ∣ m ∣ u ( n ) u ∗ ( n + m ) , ∣ m ∣ ≤ N − 1 E [ r ^ u ′ ( m ) ] = r u ( m ) V a r [ r ^ u ′ ( m ) ] = ( 1 N − ∣ m ∣ ) 2 ∑ r = − ( N − 1 − ∣ m ∣ ) N − 1 − ∣ m ∣ ( N − ∣ m ∣ − ∣ r ∣ ) [ r u 2 ( r ) + r u ( r − m ) r u ( r + m ) ] \hat r'_u(m) = \frac{1}{N-|m|}\sum_{N=0}^{N-1-|m|}u(n)u^*(n+m),|m|\le N-1 \\ E[\hat r'_u(m)] = r_u(m) \\ Var[\hat r'_u(m)] =(\frac{1}{N-|m|})^2 \sum_{r = -(N-1-|m|)}^{N-1-|m|}(N-|m|-|r|) [r^2_u(r)+r_u(r-m)r_u(r+m)] r^u′​(m)=N−∣m∣1​N=0∑N−1−∣m∣​u(n)u∗(n+m),∣m∣≤N−1E[r^u′​(m)]=ru​(m)Var[r^u′​(m)]=(N−∣m∣1​)2r=−(N−1−∣m∣)∑N−1−∣m∣​(N−∣m∣−∣r∣)[ru2​(r)+ru​(r−m)ru​(r+m)]

结论

有偏估计的方差和均方误差比无偏估计的小，并且不会带来负的功率谱估计值，因此通常使用有偏估计式来估计自相关序列。

所谓经典方法，就是通过有偏估计式 r ^ u ( m ) \hat r_u(m) r^u​(m)进行傅里叶变换来计算出功率谱。

（1）间接法

所谓间接法就是先求出自相关序列的有偏估计式，然后对有偏估计式进行傅里叶变换，由维纳辛钦定理得到功率谱的估计，这种方法最早由Blackman和Turkey两个人提出的，又叫BT法。

这种方法通常需要对估计式加窗进行修正（这个时候其实也可以认为是在平滑周期图），减小自相关序列截断的影响（从无限到有限的截断）。加窗的本质就是在减少频谱泄露。

同时注意窗型和窗长：窗长和窗型同时影响分辨率，窗型主要影响谱估计的方差。

（2）直接法

直接法是指，先求出序列的离散傅里叶变换，然后求模平方，得到能量谱，再除以观察时间就得到功率谱，这种方法又称周期图法。为什么称为周期图（periodogram）法呢，这个和离散傅里叶变换的原理有关，离散傅里叶变换就是认为现在观察的时间已经足够长了，采样频率也足够高了，那么再去观察的话，不过是对现有观察序列的重复，那么也就是说现有序列是这个本该是无限长序列的一个周期。

显然这个假设在大多数时候是不符合实际的，而且谱分辨率也是低的，但是由于可以采用FFT算法，它的计算效率很高，在对谱分辨率要求不高的地方，这种方法还是很常用的。不过，后面还有对周期图法的改进方案。

（3）Bartlett法

这种方法的思路很简单，就是将原本的序列切割成两段，每段各自按照periodogram的方法求一个功率谱出来，然后把两个功率谱加起来除以二，从而将缩小原始的周期图法的方差。显然，可以切割成多段。但是这种切割是要付出代价的，代价就是谱分辨率降低。为什么会降低呢？因为原序列N个点，切割成两段，每段就是N/2个点，在采样率相同的情况下，N点FFT当然比N/2点FFT的分辨率高了（再深入的话就偏题了）。

（4）Welch法

Welch法是结合了Bartlett法和Blackman-Turkey法，并加以改进的一种方法。第一个改进的地方是切割成重叠的分段，用房地产行业的术语来说就是公摊面积。原本N个点分成两截，每截N/2，现在一人拿出N/4来公摊，那就是每截0.75N了，不过有0.5N是重复的，也就是帧长0.75N，帧移只有0.25N，重叠了0.5N。这样在一定程度上解决了谱分辨率和谱估计方差的矛盾，但有重叠，那么互相关性就会增强嘛，也就是说谱估计方差没有Bartlett法那么理想。

同时，也采用了BT法的思路，给每段序列加窗，改进的地方在于加了一个归一化因子——窗函数的能量，这样使得任何窗函数都可以得到非负的谱估计。

1、《通信原理（第二版）》李晓峰等著。

2、《谱估计与自适应信号处理教程》赵晓晖编著。

3、《MATLAB辅助现代工程数字信号处理》李益华主编。