经典功率谱估计的原理及MATLAB仿真（自相关函数BT法、周期图法、bartlett法、welch法）

---

前言

经典功率谱估计方法包括BT法（对自相关函数求傅里叶变换求功率谱）、周期图法、Bartlett法（分段求平均）、welch法（有重合分段求平均）。本文在总结各种方法的原理后将在MATLAB平台上进行仿真，完成对功率谱密度的估计。

---

提示：以下是本篇文章正文内容，转载请附上链接！

一、BT法

设$u_N(0),u_N(1),\cdot \cdot \cdot ,u_N(N-1)$ 为广义平稳随机过程$u(n)$ 的$N$ 个观测值（为描述方便，本章观测样本的序号为从0到$N-1$的整数)，且设$u_N(n)$其他时刻的值为零，则$u_N(n)$可表示为

$$u_N(n)=\{\begin{array}{ll}u(n),& 0⩽n⩽N-1\\ & \\ 0,& \text{其他}\end{array}\text{\hspace{1em}(1)}$$

$u(n)$ 的自相关函数$r(m)$ 可用时间平均进行估计，即

$$\hat{r}(m)=\frac{1}{N}\sum _{n=0}^{N-1}{u}_N(n)u_N^{\ast }(n-m),\mid m\mid ⩽N-1\text{\hspace{1em}(2)}$$
根据维纳-辛钦定理，对由式(2) 估计得到的自相关函数$\hat{r}(m)$求傅里叶变换，可得功率谱的估计为

$${\hat{S}}_{\mathrm{BT}}\left(\omega \right)=\sum _{m=-\infty }^{\infty }\hat{r}\left(m\right){\mathrm{e}}^{-\mathrm{j}\omega m}=\sum _{m=-N+1}^{N-1}\hat{r}\left(m\right){\mathrm{e}}^{-\mathrm{j}\omega m}\text{\hspace{1em}(3)}$$

考虑到自相关函数$\hat{r}(m)$ 在$\mid m\mid >N-1$ 时为零，且$\hat{r}(m)$ 在$m$ 接近$N-1$ 时性能较差， 式(3)经常表示为

$${\hat{S}}_{\mathrm{BT}}\left(\omega \right)=\sum _{m=-M}^M\hat{r}\left(m\right){\mathrm{e}}^{-\mathrm{j}\omega m},0⩽M⩽N-1\text{\hspace{1em}(4)}$$

以此结果作为对理论功率谱$S(\omega )$ 的估计，因为这种方法估计出的功率谱是通过自相关函数间接得到的，所以此方法又称为间接法。

当时延m≥0时，r(m)的均值可以表示为

E ⁡ { r ^ ( m ) } = E { 1 N ∑ n = 0 N − 1 u N ( n ) u N ∗ ( n − m ) } = 1 N ∑ n = m N − 1 E { u N ( n ) u N ∗ ( n − m ) } = 1 N ∑ n = m N − 1 r ( m ) = N − m N r ( m ) (5) \begin{aligned} \operatorname{E}\{\hat{r}(m)\}& =\text{ E}\left\{\frac1N\sum_{n=0}^{N-1}u_N\left(n\right)u_N^*\left(n-m\right)\right\} \\ &=\frac1N\sum_{n=m}^{N-1}E\{u_N(n)u_N^*(n-m)\} \\ &=\frac1N\sum_{n=m}^{N-1}r\left(m\right) \\ &=\frac{N-m}Nr\left(m\right)\tag{5} \end{aligned} E{r^(m)}​= E{N1​n=0∑N−1​uN​(n)uN∗​(n−m)}=N1​n=m∑N−1​E{uN​(n)uN∗​(n−m)}=N1​n=m∑N−1​r(m)=NN−m​r(m)​(5)

考虑到$m$ 的取值可正可负，所以有

E ⁡ { r ^ ( m ) } = N − ∣ m ∣ N r ( m ) (6) \operatorname{E}\{\hat{r}(m)\}=\frac{N-\mid m\mid}{N}r(m) \tag{6} E{r^(m)}=NN−∣m∣​r(m)(6)

BT法结论：
$\hat{r}(m)$ 是$r(m)$ 的渐近无偏估计。
$\hat{r}(m)$ 是$r(m)$ 的渐近一致估计。

二、周期图法

由于随机过程$u(n)$ 的$N$个观测值$u_N(n)$ 是确定信号，对其进行傅里叶变换，得

$$U_N(\omega )=\sum _{n=0}^{N-1}{u}_N(n){\mathrm{e}}^{-\mathrm{j}\omega n}\text{\hspace{1em}(7)}$$

根据帕塞瓦尔定理，上式的模的平方是确定信号$u_N(n)$的能量谱，对能量谱除以持续时间$N$, 其结果应是$u_N(n)$ 的功率谱估计，将其作为随机过程$u(n)$的功率谱的估计，表示为

$${\hat{S}}_{\mathrm{PER}}(\omega )=\frac{1}{N}\mid U_N(\omega ){\mid }^2\text{\hspace{1em}(8)}$$

该方法称为功率谱估计的周期图法（periodogram)。因为这种功率谱估计方法是直接通过观测数据的傅里叶变换求得的，称之为直接法。

周期图法结论：
$\hat{r}(m)$ 是$r(m)$ 的渐近无偏估计。
$\hat{r}(m)$ 不是$r(m)$ 的一致估计，估计方差趋近于一个有限值。

三、Bartlett法

Bartlett法的基本步骤是，将$N$ 点的观测数据$u_N(n)$ 分为$L$ 段，每段数据的长度为$M$, 整数$L$ 和$M$ 满足

$$L=\frac{N}{M}\text{\hspace{1em}(9)}$$

第$i$ 段数据加矩形窗后，有

$$u_N^i(n)=u(n+(i-1)M)w_M^{(R)}(n),0⩽n⩽M-1,1⩽i⩽L\text{\hspace{1em}(10)}$$
其中$,w_M^{(R)}(n)$ 是长度为$M$ 的矩形窗。

对于每段数据$u_{\mathbb{N}}^i\left(n\right)$, 先利用周期图法求得其功率谱${\hat{S}}_{\mathrm{PER}}^i(\omega )$, 限

$${\hat{S}}_{\mathrm{PER}}^i(\omega )=\frac{1}{M}|\sum _{n=0}^{M-1}{u}_N^i(n){\mathrm{e}}^{-\mathrm{j}\omega n}{|}^2,1⩽i⩽L\text{\hspace{1em}(11)}$$

然后计算各段功率谱估计的平均，得到平均周期图${\overline{\mathcal{S}}}_{\mathrm{PER}}\left(\omega \right)$, 即

$${\overline{S}}_{\mathrm{PER}}\left(\omega \right)=\frac{1}{L}\sum _{i=1}^L{\hat{S}}_{\mathrm{PER}}^i\left(\omega \right)=\frac{1}{LM}\sum _{i=1}^L|\sum _{n=0}^{M-1}{u}_N^i\left(n\right){\mathrm{e}}^{-\mathrm{j}\omega n}{|}^2\text{\hspace{1em}(12)}$$

Bartlett法结论：
Bartlett功率谱估计使得数据点数从$N$ 减小为$M=N/L$, 于是窗函数的频谱宽度增大为周期图法的$L$ 倍，因此频率分辨率下降为原来的$1/L$ 。
Bartlett功率谱估计的方差减小为周期图法的$1/L$, 因此Bartlett功率谱估计较周期图法的结果更为平滑。

四、welch法

Welch法也对$N$点的信号$u_N(n)$ 进行分段，但允许分段时每段信号样本重叠。若每段的样本长度为$M$, 信号被分为$L$ 段，取相邻两段的信号样本重叠$50\%$ 则$L$满足

$$L=\frac{N-M/2}{M/2}\text{\hspace{1em}(13)}$$

将每段信号$u_N^i(n)$ 和窗函数$w(n)$（可以采用矩形窗、三角窗、汉宁窗或海明窗等）相乘，然后按式(11) 得到每段信号的功率谱估计为

$${\hat{S}}_{\mathrm{PER}}^i(\omega )=\frac{1}{M}|\sum _{n=0}^{M-1}{u}_N^i(n)w(n){\mathrm{e}}^{-\mathrm{j}\omega n}{|}^2\text{\hspace{1em}(14)}$$

由此得到修正的周期图${\overline{\mathcal{S}}}_{\mathrm{PER}}\left(\omega \right)$ 为
$${\overline{\mathcal{S}}}_{\mathrm{PER}}\left(\omega \right)=\frac{1}{L}\sum _{i=1}^L{\hat{S}}_{\mathrm{PER}}^i(\omega )=\frac{1}{LM}\sum _{i=1}^L|\sum _{n=0}^{M-1}{u}_N^i(n)w(n){\mathrm{e}}^{-\mathrm{j}\omega n}{|}^2\text{\hspace{1em}(15)}$$

Welch法结论：
Welch方法允许分段数据样本的重叠，于是可以得到更多的周期图估计，从而进一步减小估计的功率谱密度的方差。而与使用矩形窗实现数据分段的Bartlett方法相比，Welch方法可以使用多种窗函数。通过窗函数加权，可以减小每个分段起始和结束位置附近的样本在估计中的作用，从而减小了相邻样本段之间的相关性。所以，Welch 方法可以更好地控制功率谱密度估计的方差特性。

五、MATLAB仿真

设随机过程$u_N(n)$由3个实正弦信号加噪声构成，为
$$u_N(n)=\sum _{k=1}^3{s}_k(n)+v(n)\text{\hspace{1em}(16)}$$
其中，$s_k(n)=A_k\cos (2\pi f_{k}n+{\phi }_k),k=1,2,3$,其归一化频率分别是$f_1=0.1,f_2=0.25$和$f_3=0.27,{\phi }_k$是相互独立并在$\left[\begin{array}{c}0,2\pi \end{array}\right]$上服从均匀分布的随机相位；$v(n)$是均值为0、方差${\sigma }^2=1$的实高斯白噪声序列。3个信号的信噪比分别为${\mathrm{SNR}}_1=30$dB$,{\mathrm{SNR}}_2=30$dB$,{\mathrm{SNR}}_3=27$dB。
参数设置如下：













通过以上结果可以看出点数对信号分辨率的影响，点数加大时，可以提高分辨率，我这里就不贴图了，感兴趣的可以自己下载代码验证。

六、MATLAB详细代码

经典功率谱估计超详细的MATLAB代码（BT法和周期图法和bartlett法和welch法）

总结

例如：以上就是今天要讲的内容，本文介绍了BT法、周期图法、Bartlett法、welch法四种经典功率谱估计方法，并在MATLAB上面完成了功率谱估计仿真。