第五章 功率谱估计 笔记

第五章 功率谱估计

5.1 引言

所谓的谱估计或功率谱估计，就是用已观测到的一定数量的样本数据估计一个平稳随机信号的功率谱。随机信号是一种无始无终，能量无限的信号，其傅里叶变换不不收敛，没有频谱。但是其功率大多有限，功率反映单位时间内的能量，故对随机信号而言我们常用功率谱描述其频率特性，且它在各个点上的值无法先验确定，因此只能用它的各种非统计量来表征。其中自相关函数是平稳随机过程一个重要的统计量，并特别指出其傅里叶变换正是该信号的功率谱密度。

谱估计：确定信号在特定频率处的功率含量。
$$x(n)⟶H(e^{j\omega })⟶y(n)$$

$$\begin{array}{rl}{P}_{yy}(\omega )& =P_{xx}(\omega )|H(e^{j\omega }){|}^2\\ E(y^2)& =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{P}_{yy}(\omega )d\omega \\ & =\frac{1}{2\pi }{\int }_{-\pi }^{\pi }{P}_{xx}(\omega )|H(e^{j\omega }){|}^{2}d\omega \\ & \approx 2\cdot \frac{1}{2\pi }{P}_{xx}({\omega }_0)\Delta \omega =\frac{\Delta \omega }{\pi }{P}_x({\omega }_0)\end{array}$$

即

$$P_{xx}({\omega }_0)\approx \frac{\pi }{\Delta \omega }E[y^2]$$

$\Delta \omega$ 越小，估计越准确。

功率谱估计分为两大类

1）参数方法：模型辨识法、capon最小方差、……

2）非参数法：BT法、周期图法、……

设 $\hat{\theta }$ 是对真实值 $\theta$ 的估计值，一个好的估计应满足

1）无偏估计：$E[\hat{\theta }]=\theta$ , $B=Bia[\hat{\theta }]=\theta -E[\hat{\theta }]$ 称为偏差；

2）最小方差估计：最小化 ${\sigma }^2=E[\hat{\theta }-E^2[\hat{\theta }]]$

均方误差：$E[e^2]=E[(\theta -\hat{\theta }{)}^2]=B^2+{\sigma }^2$

一致估计：（N是数据量）

$N\to \infty ,Bia(\hat{\theta })\to 0$

$N\to \infty ,{\sigma }^2\to 0$

5.2 经典估计方法

5.2.1 相关图法（自相关函数估计）

• 相关函数的估计

假设 $x(n)$ 是各态历经的，可以用时间平均代替集合平均，${\hat{\Phi }}_{xx}^{{}^{\prime }}(m)$ 长度为 $2N-1$ 。
$${\hat{\Phi }}_{xx}^{{}^{\prime }}(m)=\frac{1}{N-|m|}\sum _{n=0}^{N-1-|m|}x(n)x(n+m),|m|\le N-1,n+m\le N-1$$

$$E[{\hat{\Phi }}_{xx}^{{}^{\prime }}(m)]=\frac{1}{N-|m|}\sum _{n}E[x(n)x(n+m)]={\Phi }_{xx}(m)无偏估计$$

这里可以直接求 ${\hat{\Phi }}_{xx}^{{}^{\prime }}(m)$ 的离散傅里叶变换，得到 $x(n)$ 的功率谱的估计值（间接法（BT法））。
$${\hat{P}}_{xx}(\omega )=\sum _{m=-M}^M{\hat{\Phi }}_{xx}^{{}^{\prime }}(m)e^{-j\omega m},M\le N-1$$

当 $m$ 接近 $N$ 时，估计的方差很大。对此问题进行改进，使用三角窗加权。（改进）
$$\begin{array}{rl}{\hat{\Phi }}_{xx}(m)& =\frac{1}{N}\sum _{n=0}^{N-1-|m|}x(n)x(n+m)\\ & =\frac{N-|m|}{N}{\hat{\Phi }}_{xx}^{{}^{\prime }}(m)\\ & =W_B(m){\hat{\Phi }}_{xx}^{{}^{\prime }}(m)\\ W_B(m)=1-\frac{|m|}{N}三角窗& \end{array}$$
这时，有偏估计
$$E[{\hat{\Phi }}_{xx}(m)]=W_B(m)E[{\hat{\Phi }}_{xx}^{{}^{\prime }}(m)]=W_B(m){\Phi }_{xx}(m)$$
当 $N\to \infty$ 时
$$\begin{gathered} E[{\hat{\Phi }}_{xx}(m)]\to {\Phi }_{xx}渐进无偏估计 \\ Var[{\hat{\Phi }}_{xx}^{{}^{\prime }}(m)]=W_B^2(m)Var[{\Phi }_{xx}^{{}^{\prime }}(m)]<Var[{\hat{\Phi }}_{xx}(m)]方差变小 \end{gathered}$$
而且 ${\hat{\Phi }}_{xx}^{{}^{\prime }}(m)$ 的均方误差也小于 ${\hat{\Phi }}_{xx}(m)$ 的均方误差。所以一般用 ${\hat{\Phi }}_{xx}^{{}^{\prime }}(m)$ 作为自相关函数的估计，而较少使用 ${\hat{\Phi }}_{xx}(m)$ 。后来为了表示方便，我们仍用 ${\hat{\Phi }}_{xx}(m)$ 表示，即
$${\hat{\Phi }}_{xx}(m)=\frac{1}{N}\sum _{n=0}^{N-|m|-1}x(n)x(n+m)$$
无论是有偏估计还是无偏估计都可以通过增大 $N$ 减小误差。

5.2.2 周期图法

• 周期图法定义（直接法）

按上一节所述的自相关函数的定义
$${\hat{\Phi }}_{xx}(m)=\frac{1}{N}\sum _{n=0}^{N-|m|-1}x(n)x(n+m)$$
定义矩形窗函数
$$d(n)=\{\begin{array}{cc}1,& 0\le n\le N-1\\ 0,& 其它\end{array}$$
令 $d(n)x(n)=y(n)$ ，上式可以简写为
$$\begin{array}{rl}{\hat{\Phi }}_{xx}(m)& =\frac{1}{N}\sum _{n=-\infty }^{\infty }d(n)x(n)d(n+m)x(n+m)\\ & =\frac{1}{N}\sum _{n=-\infty }^{\infty }y(n)y(n+m)\\ & =y(n)\ast y(-m)\end{array}$$
$y(m)$ 的傅里叶变换
$$\begin{array}{rl}\sum _{m=-\infty }^{\infty }y(m)e^{-j\omega m}& =\sum _{m=-\infty }^{\infty }d(m)x(m)e^{-j\omega m}\\ & =\sum _{m=0}^{N-1}x(m)e^{-j\omega m}\\ & =X_N(e^{-j\omega })\end{array}$$
$y(-m)$ 的傅里叶变换为 $X_N^{\ast }(e^{-j\omega })$ 根据相关图法，功率谱估计为自相关函数估计的傅里叶变换，即
$${\hat{P}}_{xx}(\omega )=\frac{1}{N}{X}_N(e^{-j\omega })X_N^{\ast }(e^{-j\omega })=\frac{1}{N}|X_N(e^{-j\omega }){|}^2$$
直接将 $X_N(\omega )$ 模的平方除以 $N$ 求得的功率谱估计的方法称为周期图法，其结果用 $I(\omega )$ 表示
$$I(\omega )={\hat{P}}_{xx}(\omega )=\frac{1}{N}|X_N(e^{-j\omega }){|}^2$$
（优点不再计算自相关函数，直接利用数据的傅里叶变换求得。）

• 周期图法的谱估计性能

为了了解周期图法的估计效果，我们来讨论它的估计均值和方差。（均值）
$$\begin{array}{rl}E[I_N(\omega )]& =E[\frac{1}{N}|X_N(e^{j\omega }){|}^2]=\frac{1}{N}E[X_N^{\ast }(e^{j\omega })X_N(e^{j\omega })]\\ & =\frac{1}{N}E[\sum _{n=0}^{N-1}x(n)e^{j\omega n}\sum _{k=0}^{N-1}x(k)e^{-j\omega k}]\\ & =\frac{1}{N}E[\sum _{n=-\infty }^{\infty }d(n)x(n)e^{j\omega n}\sum _{k=-\infty }^{\infty }d(k)x(k)e^{-j\omega k}]\\ & =\frac{1}{N}E[\sum _{n=-\infty }^{\infty }\sum _{k=-\infty }^{\infty }d(n)d(k)x(n)x(k)e^{-j\omega (k-n)}]\end{array}$$
令 $m=k-n$ ，代入上式
$$\begin{array}{rl}E[I_N(\omega )]& =\sum _{m=-\infty }^{\infty }\frac{1}{N}[\sum _{n=-\infty }^{\infty }d(n)d(n+m)E[x(n)x(n+m)]e^{-j\omega m}]\\ & =\sum _{m=-\infty }^{\infty }{q}_N(m){\varphi }_{xx}(m)e^{-j\omega m}\\ 其中，& q_N(m)=\frac{1}{N}\sum _{n=-\infty }^{\infty }d(n)d(n+m)=\{\begin{array}{cc}1-\frac{|m|}{N},& 对于|m|\le N-1\\ 0,& 其它\end{array}\end{array}$$
由于 $q_N(m)$ 为两个矩形函数的卷积，为三角窗函数，故可对上式傅里叶变换，
$$\begin{array}{rl}E[I_N(\omega )]& =Q_N(\omega )\ast P_{xx}(\omega )\\ & =\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{Q}_N(\theta )P_{xx}(\omega -\theta )d\theta \end{array}$$
由上式可知，除非 $Q_N(\omega )$ 为 $\delta$ 函数，否则 $E[I_N(\omega )]$ 将不等于 $P_{xx}(\omega )$ ，故周期图作为功率谱估计是有偏估的。普分辨率取决于主瓣宽度， $N$ 越大（时间越长）谱分辨率越高。

当 $N\to \infty$ 时，
$$\begin{gathered} \underset{N\to \infty }{\lim }{q}_N(m)=\underset{N\to \infty }{\lim }(1-\frac{|m|}{N})=1 \\ {Q}_{N\to \infty }(\omega )=\delta (\omega ) \\ 故\underset{N\to \infty }{\lim }E[I_N(\omega )]=P_{xx}(\omega )渐进无偏估计 \end{gathered}$$

（方差）下面讨论方差，假设序列 $x(n)(0\le n\le N-1)$ 是均值为零、方差为 ${\sigma }_x^2$ 的实高斯白噪声序列，按方差定义
$$var[I_N(\omega )]=E[I_N^2(\omega )]-E^2[I_N(\omega )]\approx P_{xx}^2(\omega )[1+(\frac{sinN\omega }{Nsin\omega }\left)\right]$$
当 $N\to \infty$ 时
$$var[I_N(\omega )]\to P_{xx}^2(\omega )$$
明显 $N\to \infty$ 时，方差不趋向于零，周期图不是一致估计，方差大。

（改进）定理 如果 $x_1,x_2,\cdots ,x_L$ 是互不相关的随机变量，且每一个有均值 $\mu$ 及方差 ${\sigma }^2$ ，则它们的数学平均 $\overline{x}=(x_1+x_2+\cdots +x_L)/L$ 的均值等于 $\mu$ ，方差为 ${\sigma }^2/L$ 。
$$\begin{array}{rl}E[\overline{x}]& =\frac{1}{L}E(x_1+x_2+\cdots +x_L)=\frac{1}{L}\cdot L\mu =\mu \\ var[\overline{x}]& =E[(\overline{x}-E[\overline{x}]{)}^2]=E[{\overline{x}}^2]-E^2[\overline{x}]\\ & =\frac{1}{L^2}E[(x_1+x_2+\cdots +x_L{)}^2]-{\mu }^2\\ & =\frac{1}{L^2}\{E[x_1^2+x_2^2+\cdots +x_L^2]+\sum _{j=1}^L\sum _{i=1,i\ne j}^{L}E[x_i{x}_j]\}-{\mu }^2\end{array}$$
由于 $x_1,x_2,\cdots ,x_L$ 不相关
$$\begin{array}{rl}\sum _{j=1}^L\sum _{i=1,i\ne j}^{L}E[x_i{x}_j]& =\sum _{j=1}^{L}E[x_j]\sum _{i=1,i\ne j}^{L}E[x_i]\\ & =L\mu (L-1)\mu \\ & =L^2{\mu }^2-L{\mu }^2\end{array}$$
故
$$\begin{array}{rl}var[\overline{x}]& =\frac{1}{L^2}[\sum _{i=1}^{L}E[x_i^2]+L^2{\mu }^2-L{\mu }^2]-{\mu }^2\\ & =\frac{1}{L^2}[\sum _{i=1}^{L}E[x_i^2]-L{\mu }^2]\\ & =\frac{1}{L^2}\{[E(x_1^2)-E^2(x_1)]+[E(x_2^2)-E^2(x_2)]+\cdots +[E(x_L^2)-E^2(x_L)]\}\\ & =\frac{1}{L^2}L{\sigma }^2=\frac{{\sigma }^2}{L}\end{array}$$
具有 $L$ 个独立同分布随机变量平均的方差是每个随机变量方差的 $\frac{1}{L}$ ，当 $L\to \infty$ 时， $var[\overline{x}]\to 0$ ，可达到一致估计的目的。因而将若干个独立估计值进行平均可有效降低估计量的方差。

平均周期图法

将 $x(n)$ 分成 $L$ 段，每段有 $M$ 个样本，则 $N=LM$ ，第 $i$ 段序列的样本可写成
$$x^i(n)=x(n+iM-M)0\le n\le M-1,1\le i\le L$$
第 $i$ 段的周期图为
$$I_M^i(\omega )=\frac{1}{M}|\sum _{n=0}^{M-1}{x}^i(n)e^{-j\omega n}{|}^2$$
如果 $m>M$ ，${\varphi }_{xx}(m)$ 很小，则可假定各段的周期图 $I_M^i(\omega )$ 是互相独立的，总的谱估计可定义为 $L$ 段周期图的平均，即
$${\hat{P}}_{xx}(\omega )=\frac{1}{L}\sum _{i=1}^L{I}_M^i(\omega )$$
如各段数据互相独立，则估计方差将只有原来的 $1/L$ 所以当 $L$ 越大，方差越小，越接近一致估计，但是随着 $L$ 增大，$M$ 减小，分辨率下降，将使估计变成有偏。故实际中必须兼顾分辨率与方差的要求适当选择 $L$ 和 $M$ 的值。

具体方法（Welch法）

1）先将 $N$ 个数据分成 $L$ 段，每段 $M$ 个数据，$N=LM$ ；

2）选择适当的窗函数 $\omega (n)$ 依次对每段数据作相应的加权，然后确定每一段的周期图
$$I_M^i(\omega )=\frac{1}{MU}|\sum _{n=0}^{M-1}{x}^i(n)\omega (n)e^{-j\omega n}|,1\le i\le L$$
这里 $U=\frac{1}{M}{\sum }_{n=0}^{M-1}{\omega }^2(n)$ 为归一化因子；

3）对分段周期图进行平均得到功率谱估计
$${\hat{P}}_{xx}(\omega )=\frac{1}{L}\sum _{i=1}^L{I}_M^i(\omega )$$

优点：

• 对窗函数没有特殊要求，均可使谱估计非负；
• 分段时减少因分段数增加给分辨率带来的影响，可使各段之间有重叠，例如 50% 重叠。

5.3 谱估计的参数化模型方法

我们可以将随机过程看成由白噪声 $\omega (n)$ 经一线性系统形成的
$$\begin{gathered} \omega (n)⟶H(z)⟶x(x) \\ H(z)=\frac{B(z)}{A(z)}=\frac{{\sum }_{k=0}^q{b}_k{z}^{-k}}{{\sum }_{k=0}^p{a}_k{z}^{-k}} \end{gathered}$$
当输入白噪声的功率谱密度 $P_{xx}(\omega )={\sigma }_{\omega }^2$ 时，输出的功率谱密度为
$$P_{xx}(\omega )={\sigma }_{\omega }^2|H(e^{j\omega }){|}^2={\sigma }_{\omega }^2|\frac{B(e^{j\omega })}{A(e^{j\omega })}{|}^2$$
下面假设 $a_0=1,b_0=1$ 我们将线性系统分为三种情况。

• AR模型

除 $b_0=1$ 外，所有的 $b_k(1\le k\le q)$ 均为零，此时系统函数为
$$H(z)=\frac{1}{A(z)}=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$
其差分方程为
$$x(n)=-\sum _{k=1}^p{a}_{k}x(n-k)+\omega (n)$$
模型输出功率谱
$$P_{xx}(\omega )=\frac{{\sigma }_{\omega }^2}{|A(e^{j\omega }){|}^2}=\frac{{\sigma }_{\omega }^2}{|1+{\sum }_{k=1}^p{a}_k{e}^{-j\omega k}{|}^2}$$

• MA模型

除 $a_0=1$ 外，所有的 $a_k(1\le k\le p)$ 均为零，此时系统函数为
$$H(z)=B(z)=1+\sum _{k=1}^q{b}_k{z}^{-k}$$
其差分方程为
$$x(n)=\sum _{k=1}^q{a}_{k}x(n-k)$$
模型输出功率谱
$$P_{xx}(\omega )={\sigma }_{\omega }^2|B(e^{j\omega }){|}^2={\sigma }_{\omega }^2|1+\sum _{k=1}^q{b}_k{e}^{-j\omega k}{|}^2$$

• ARMA模型

设 $a_0=1,b_0=1$ 其余的 $a_k,b_k$ 不全为零，这是一个 “极点、零点” 模型，称为 $ARMA(p,q)$ 模型，或直接简称 ARMA 模型。

5.4 自回归（AR）模型方法

5.4.1 AR模型的 Yule-Walker 方法

有前面讨论可知，$p$ 阶 AR 模型的系统函数和差分方差分别为

除 $b_0=1$ 外，所有的 $b_k(1\le k\le q)$ 均为零，此时系统函数为
$$H(z)=\frac{1}{A(z)}=\frac{1}{1+{\sum }_{k=1}^p{a}_k{z}^{-k}}$$
其差分方程为
$$x(n)=-\sum _{k=1}^p{a}_{k}x(n-k)+\omega (n)$$
模型输出功率谱
$$P_{xx}(\omega )={\sigma }_{\omega }^2|B(e^{j\omega }){|}^2={\sigma }_{\omega }^2|1+\sum _{k=1}^q{b}_k{e}^{-j\omega k}{|}^2$$

若已知参数 $a_1,a_2,\cdots ,a_p$ 及 ${\sigma }_w^2$ ，就可以得到信号的功率谱估计。现在我饿们研究这些参数与自相关函数的关系。将 AR 模型的差分方程代入 $x(n)$ 的自相关函数表达式，得
$$\begin{array}{rl}{\varphi }_{xx}(m)& =E[x(n)x(n+m)]\\ & =E\{x(n)[-\sum _{k=1}^p{a}_{k}x(n+m-k)+\omega (n+m)\left]\right\}\\ & =-\sum _{k=1}^p{a}_k{\varphi }_{xx}(m-k)+E[x(n)\omega (n+m)]\end{array}$$
由差分方程可知，$x(n)$ 只与 $\omega (n)$ 相关而与 $\omega (n+m)$ 无关 $(m\ge 1)$ ，故有
$$\begin{gathered} E[x(n)\omega (n+m)]=E[\omega (n)\omega (n+m)]=\{\begin{array}{cc}{\sigma }_w^2,& m=0\\ 0,& m>0\end{array} \\ {\varphi }_{xx}(m)=\{\begin{array}{cc}-{\sum }_{k=1}^p{a}_k{\varphi }_{xx}(m-k),& m>0\\ -{\sum }_{k=1}^p{a}_k{\varphi }_{xx}(m-k)+{\sigma }_w^2,& m=0\end{array} \end{gathered}$$
写成矩阵形式，得
$$\left[\begin{array}{ccccc}{\varphi }_{xx}(0)& {\varphi }_{xx}(-1)& {\varphi }_{xx}(-2)& \cdots & {\varphi }_{xx}(-p)\\ {\varphi }_{xx}(1)& {\varphi }_{xx}(0)& {\varphi }_{xx}(-1)& \cdots & {\varphi }_{xx}[-(p-1)]\\ {\varphi }_{xx}(2)& {\varphi }_{xx}(1)& {\varphi }_{xx}(0)& \cdots & {\varphi }_{xx}[-(p-2)]\\ ⋮& ⋮& ⋮& & ⋮\\ {\varphi }_{xx}(p)& {\varphi }_{xx}(p-1)& {\varphi }_{xx}(p-2)& \cdots & {\varphi }_{xx}(0)\end{array}\right]\left[\begin{array}{c}1\\ a_1\\ a_2\\ ⋮\\ a_p\end{array}\right]=\left[\begin{array}{c}{\sigma }_w^2\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]$$
上式就是 AR 模型得 Yule-Walker 方程。

5.4.2 AR 谱估计与线性预测谱估计等效

线性预测器是在过去 $p$ 个样本 $x(n-1),x(n-2),\cdots ,x(n-p)$ 的基础上预测当前值 $x(n)$ 的，即
$$\hat{x}(n)=-\sum _{k=1}^p{a}_{pk}x(n-k)$$
由于两个 Yule-Walker 方程相同，当相应的自相关函数值相同时，它们的解也必然相同，即 $a_k=a_{pk}(k=1,2,\cdots ,p){\epsilon }_{min}={\sigma }_w^2$ 也就是说最佳线性预测器参数恰好等于 AR 模型参数，最小预测误差功率也就等于激励的白噪声功率。

用时刻 $n$ 以前的 $p$ 个取样样本数据预测 $x(n)$ ，可做出滤波解释。若将 $x(n)$ 作为输入，滤波器传输函数为
$$H^{-1}(z)=\sum _{k=0}^p{a}_{pk}{z}^{-1}$$
那么输出将是预测误差 $e(n)$ 即白噪声 $\omega (n)$ 。

​

5.4.3 最大熵谱估计及其与 AR 谱估计的等效性

最大熵谱估计则是基于一段已知的自相关序列进行外推，以得到未知的自相关值。这样因对自相关序列加窗而使谱估计特性变坏的弊端就被去掉了。

方法：外推后的自相关序列所对应的时间序列应当具有最大熵，通过这种外推得到的自相关序列求出的谱则称为最大熵谱估计（MESE）。

熵：熵代表一不确定度，最大熵为最大不确定度。最大熵法外推的自相关函数序列将是随机或是不可预测的，或者说它的谱将是平坦的或者是最白的

• 当 $x$ 为离散值时，熵 $H$ 定义为

$$H=-\sum _l{p}_{i}ln{p}_i$$

• 当 $x$ 为连续值时，熵被定义为

$$H=-\int p(x)lnp(x)dx=-E[ln(p(x))]$$

当 { x ( n ) } \{x(n) \} {x(n)} 为零均值、高斯分布的随机过程时， $N$ 维联合概率函数为
$$P(x_1,x_2,\cdots ,x_N)=(2\pi {)}^{-N/2}(det\Phi (N){)}^{\frac{1}{2}}exp(-\frac{1}{2}{X}^T[\Phi (N){]}^{-1}X)$$
将一维熵的结果推广到 $N$ 维，得到 $N$ 维高斯分布信号的熵为
$$H=ln[(2\pi e{)}^{N/2}(det\Phi (N){)}^{1/2}]$$
要使熵最大，就要求 $det(\Phi (N))$ 最大。如果已知 ${\varphi }_{xx}(0),{\varphi }_{xx}(1),\cdots ,{\varphi }_{xx}(N)$ ，欲求 ${\varphi }_{xx}(N+1)$ 时，应该在最大熵条件下，求取下式所列方程的解：
$$\frac{\partial H}{\partial {\varphi }_{xx}(N+1)}=0$$
由于最大熵与最大 $det[\Phi (N+1)]$ 等效，因而可求解
$$\frac{\partial det[\Phi (N+1)]}{\partial {\varphi }_{xx}(N+1)}=0$$
最后 AR 模型的结果与按照最大熵外推 ${\varphi }_{xx}(N+1)$ 得到的结果完全相同。证明了当 $x(n)$ 为高斯分布时的最大熵谱估计与 AR 模型法是等价的。

5.4.4 Levinson-Durbin 递推算法

一阶 AR 模型的 Yule-Walker 矩阵方程为
$$\left[\begin{array}{cc}{\varphi }_{xx}(0)& {\varphi }_{xx}(1)\\ {\varphi }_{xx}(1)& {\varphi }_{xx}(0)\end{array}\right]\left[\begin{array}{c}1\\ a_{11}\end{array}\right]=\left[\begin{array}{c}{\sigma }_1^2\\ 0\end{array}\right]$$
解方程中的未知参数 $a_{11}$ 和 ${\sigma }_1^2$ 为
$$\begin{gathered} a_{11}=-\frac{{\varphi }_{xx}(1)}{{\varphi }_{xx}(0)} \\ {\sigma }_1^2=(1-|a_{11}{|}^2){\varphi }_{xx}(0),|a_{11}|<1 \end{gathered}$$
二阶
$$\left[\begin{array}{ccc}{\varphi }_{xx}(0)& {\varphi }_{xx}(1)& {\varphi }_{xx}(2)\\ {\varphi }_{xx}(1)& {\varphi }_{xx}(0)& {\varphi }_{xx}(1)\\ {\varphi }_{xx}(2)& {\varphi }_{xx}(1)& {\varphi }_{xx}(0)\end{array}\right]\left[\begin{array}{c}1\\ a_{21}\\ a_{22}\end{array}\right]=\left[\begin{array}{c}{\sigma }_2^2\\ 0\\ 0\end{array}\right]$$
参数
$$\begin{array}{rl}{a}_{22}& =-\frac{{\varphi }_{xx}(0){\varphi }_{xx}(2)-{\varphi }_{xx}^2(1)}{{\varphi }_{xx}^2(0)-{\varphi }_{xx}^2(1)}=-\frac{{\varphi }_{xx}(2)+a_{11}{\varphi }_{xx}(1)}{{\sigma }_1^2}\\ a_{21}& =-\frac{{\varphi }_{xx}(0){\varphi }_{xx}(1)-{\varphi }_{xx}(1){\varphi }_{xx}(2)}{{\varphi }_{xx}^2(0)-{\varphi }_{xx}^2(1)}=a_{11}+a_{22}{a}_{11}\\ {\sigma }_2^2& =(1-|a_{22}{|}^2){\sigma }_1^2\end{array}$$
以此得到类推公式
$$\begin{array}{rl}{a}_{kk}& =-\frac{{\varphi }_{xx}(k)+{\sum }_{l=1}^{k-1}{a}_{k-1,l}{\varphi }_{xx}(k-l)}{{\sigma }_{k-1}^2}\\ a_{ki}& =a_{k-1,i}+a_{kk}{a}_{k-l,k-i}\\ {\sigma }_k^2& =(1-|a_{kk}{|}^2){\sigma }_{k-1}^2,{\sigma }_0^2=\varphi (0)\end{array}$$

• AR 模型谱估计的优点

1. AR 模型是一个有理分式，AR 模型谱比经典法的功率谱光滑。

2. 分辨率高，AR 模型隐含着数据和自相关的外推

$$\begin{gathered} \hat{x}(n)=-\sum _{k=1}^p{a}_{k}x(n-k) \\ {\varphi }_{xx}(m)=-\sum _{k=1}^p{a}_k{\varphi }_{xx}(m-k),m\ge 1 \end{gathered}$$

• AR 模型谱估计的缺点

1. 分辨率明显受信噪比影响，低信噪比时，分辨率下降明显，对噪声敏感。
2. 谱估计的质量受阶次影响，而阶次难以确定。

• 误差

前向预测
$$\begin{gathered} \hat{x}(n)=f_{1}x(n-1)+f_{2}x(n-2)f_1=-a_1,f_2=-a_2 \\ {e}_f(n)=x(n)-\hat{x}(n) \end{gathered}$$
后向预测
$$\begin{gathered} \hat{x}(n-2)=b_{1}x(n-1)+b_{2}x(n) \\ {e}_b(n)=x(n-2)-\hat{x}(n-2) \end{gathered}$$
故 $E[e_f^2(n)]=E[e_b^2(n)]$ ，又
$$\begin{gathered} \hat{x}(n-1)=f\cdot x(n-1) \\ E[e^2(n)]=E[x(n)-f\cdot x(n-1){]}^2={\varphi }_{xx}(0)+f^2{\varphi }_{xx}(0)-2f\cdot {\varphi }_{xx}(1)=0 \end{gathered}$$
矩阵形式
$$\left[\begin{array}{cc}{\varphi }_{xx}(0)& {\varphi }_{xx}(1)\\ {\varphi }_{xx}(1)& {\varphi }_{xx}(0)\end{array}\right]\left[\begin{array}{c}1\\ -f\end{array}\right]=\left[\begin{array}{c}E[e^2(n)]\\ 0\end{array}\right]$$
有如下条件
$$\begin{gathered} \frac{dE[e^2(n)]}{df}=2f{\varphi }_{xx}(0)-2{\varphi }_{xx}(1)=0 \\ f=\frac{{\varphi }_{xx}(1)}{{\varphi }_{xx}(0)}=\rho (1)相关系数 \end{gathered}$$

5.5 白噪声中正弦波频率的估计及谱估计的其它方法

5.5.1 最大似然法

5.5.2 Capon 最小方差法

假设 $w^{\ast }(i)i=1,2,\cdots ,N-1$ ，为 $N$ 阶 FIR 滤波器的冲激响应函数， $x(n)$ 为一个零均值的平稳随机过程，经过有限冲激响应滤波器的输出 $y(n)$ 为
$$y(n)=\sum _{i=0}^{N-1}{w}^{\ast }(i)x(n+i)=W^{H}X$$
式中 $W=[w(0),w(1),\cdots ,w(N-1){]}^T$ 为滤波器系数矢量； $X=[x(n),x(n+1),\cdots ,x(n+N-1){]}^T$ 为输入信号矢量。 $y(n)$ 的平均功率（即方差）为
$$E[|y(n){|}^2]=E[W^{H}X{X}^{H}W]=W^H{\Phi }_{x}W$$
式中 ${\Phi }_x=E[X{X}^H]$ 为输入信号的自相关矩阵。

复正弦可表示为
$$x(n)=Aexp(j{\omega }_{i}n)$$
所以输入信号矢量为
$$\begin{gathered} X=Aexp(j{\omega }_{i}n)[1,exp(j{\omega }_i),\cdots ,exp(j(N-1){\omega }_i){]}^T=x(n)e \\ e=[1,exp(j{\omega }_i),\cdots ,exp(j(N-1){\omega }_i){]}^T \end{gathered}$$
所以 $y(n)=W^{H}X=x(n)W^{H}e$ 为了使复正弦 ${\omega }_i$ 完全通过，必须有
$$\begin{gathered} |W^{H}e|=1 \\ {W}^{H}e{e}^{H}W=I \end{gathered}$$
由于实际信号总包含有噪声，求滤波器系数矢量，使滤波器输出平均功率 $E[|y(n){|}^2]$ 最小，这是一个条件极值问题，可以引入 Lagrange 乘子构成代价函数 $J$
$$J=W^H{\Phi }_{x}W-\lambda (W^{H}e{e}^{H}W-1)$$
令 $\frac{\partial J}{\partial W}=0$ 可得
$$W^H{\Phi }_x=\lambda W^{H}e{e}^H$$
两边同时取复共轭转置，而且注意到 ${\Phi }_x^H={\Phi }_x$
$${\Phi }_{x}W=\lambda e{e}^{H}W$$
又因为 $e^{H}W$ 是一个模为 1 的复数，可以表示为
$$\begin{gathered} e^{H}W=exp(j\phi ) \\ W=\lambda exp(j\phi ){\Phi }_x^{-1}e \end{gathered}$$
两边同时右乘 $W$ 得 $W^H\Phi W=\lambda W^{H}e{e}^{H}W=\lambda$ 变换得
$$\begin{gathered} \lambda =\frac{1}{e^H{\Phi }_x^{-1}e} \\ W=\frac{exp(j\omega ){\Phi }_x^{-1}e}{e^H{\Phi }_x^{-1}e} \end{gathered}$$
因为 $exp(j\phi )$ 只产生相移，并不影响功率谱得估计，所以可以略去
$$W=\frac{{\Phi }_x^{-1}e}{e^H{\Phi }_x^{-1}e}$$
此时输出功率最小，输出平均功率为
$$E[|y(n){|}^2{]}_{min}=W^H\Phi W=\frac{1}{e^H{\Phi }_x^{-1}e}$$
因此 Capon 方法的功率谱估计表达式为
$${\hat{P}}_{Capon}({\omega }_i)=\frac{1}{e^H{\hat{\Phi }}_x^{-1}e}$$
式中 ${\hat{\Phi }}_x^{-1}$ 为自相关矩阵 ${\Phi }_x$ 的估计。

用最优滤波器 $w_{opt}$ 的输出功率 $E[|y(n){|}^2]$ 作为 $x(n)$ 在 ${\omega }_0$ 处的功率谱密度 $P_{xx}({\omega }_0)$ 的估计，对于任意的 $\omega$ 有
$${\hat{P}}_{Capon}({\omega }_i)=\frac{1}{e^H{\hat{\Phi }}_x^{-1}e}$$

$$P_{Capon}(\omega )=\sum _{m=0}^p\frac{1}{P_{AR}^{(m)}(\omega )}m阶次调和平均$$

Capon 最小方差谱估计的分辨率低于同阶次的 AR 模型谱估计。