【DSP】12 平稳随机信号以及谱估计
1.平稳随机信号
1.1. 引言——关于随机变量
由概率论可知，我们可以用一个随机变量 \(X\) 来描述自然界中的随机事件，若\(X\)的取值是连续的，则\(X\)是连续型随机变量。若 \(X\)的取值是离散的，则\(X\)是离散型随机变量，如服从二项式分布、泊松分布的随机变量。对随机变量\(X\)，一般用它的分布函数、概率密度及数字特征来描述。
随机信号的特征描述
信号随时间变化不具有明确的规律性，不能准确预测，不能用明确的数学关系来描述，但是可以描述其统计特征，如分布函数、概率密度函数、以及数字特征。
分布函数: $$P(x)=P{X\leq x}$$
数字特征
数字特征用来从统计的角度描述随机变量的特性，首先介绍矩的概念：
表示X的m阶原点矩。
表示X的m阶中心矩。
常用的数字特征如下：
均值:\(\mu_x=\mathbb{E}\{X\}\)
方差:\(\sigma_x^2=\mathbb{E}\{|X-\mu_x|^2\}\)
均方:\(D_X^2=\psi_x^2=\mathbb{E}\{|X|^2\}\)
协方差：\(cov[X,Y]=\mathbb{E}\{(X-\mu_x)^*(Y-\mu_y)\}\)
协方差矩阵：\(\mathbf{\Sigma}=\mathbb{E}\{(\mathbf{X-\mu_x})^*(\mathbf{X-\mu_x})^T\}\)
随机信号概念分析
- 对于一个特定的时刻，如\(t=t_1\)，如果\(x(t)\)看作随机变量\(X(t)\)，显然\(x_1(t_1),x_2(t_1)\cdots x_N(t_1)\)是一个随机变量（在同一时刻做了无数次测量）
- 当t取不同值时，可以得到不同的一组随机变量，时间采样趋近无穷时，得到的就是随机过程or随机信号
- 样本无穷多，持续时间无穷长，所以，随机信号是功率信号；
- 对任一时刻\(t_j\)，的集合构成一个随机变量。 随着tj的变化，我们会得到无穷多个随机变量。 所以，随机信号是依赖于时间 (t or n)的随机变量。
在实际应用中更希望处理的是离散化的信号。若是对上述信号离散化，得到离散随机信号\(X(n)=X(nTs)\)
🐱🏍平稳随机信号
若
\(p_X(x_1,\cdots,x_N;n_1,\cdots,n_N)=p_X(x_1,\cdots,x_N;n_{1+k},\cdots,n_{N+k}) \quad \forall k\in\mathbb{Z}\)
则称是N阶平 稳 的。如果上式对\(N\in(1,\infty)\)都成立， 则称\(X(n)\)是 严 平 稳( strict-sense stationary) ， 或狭义平稳的随机信号。
:::tips
严平稳的随机信号可以说基本上不存在， 而且其定义也无法应用于实际。因此， 人们研究和应用最多的是宽平稳 （wide-sense stationary， WSS)信号， 又称广义平稳信号。对随机信号\(X(n)\)， 若其均值为常数， 即
\(\mu_x(n)=\mathbb{E}\{X(n)\}=\mu_x\)
且方差也为常数
\(\sigma^2_X(n)=\mathbb{E}\{|X(n)-\mu_X|^2\}=\sigma_X^2\)
:::
- 宽平稳信号方差也为常数
- 自相关函数与起点无关。只和两点的时间差有关
- 互协方差函数也和时间的起点无关。
自相关函数的性质
- 自相关函数在m=0处取得最大值，并且是非负的
- 若\(X(n)\)是实信号，则\(r_X(m)=r_X(-m)\),实偶函数；若\(X(n)\)为复信号，则\(r_X(m)=r_X^*(-m)\)
- 相关矩阵是非负定的
功率谱密度
一般的方法，不是对信号直接进行傅立叶变换，而是对信号的自相关函数作傅立叶变换，这时得到的不再是频谱，而是功率谱（Power Spectrum Density, PSD)。
🐙两种定义
- \(P_x(e^{j\omega})=\lim_{M\to \infty} \frac{1}{2M+1}|\sum_{n=-M}^{M}x(n)e^{-j\omega n}|^2\)
- \(P_X(e^{j\omega})=\sum_{m=-\infty}^{\infty}r_X(m)e^{j\omega n}\)
Markov过程
🐾平稳随机信号通过线性系统
两边不能直接取傅立叶变换和 Z 变换，因为输入、输出都是功率信号。
如下4个关系成立：
\(\begin{align}
r_y(m)&=r_x(m)*h(m)*h(-m)\\
P_y(e^{j\omega})&=P_x(e^{j\omega})|H(e^{j\omega})|^2\\
r_{xy}(m)&=r_x(m)*h(m)\\
P_{xy}(e^{j\omega})&=P_x(e^{j\omega})H(e^{j\omega})
\end{align}\)
🎶各态遍历性
在具备一定的条件下，观察时间足够长的平稳过程的一个样本函数\(x(n,i)\)的“时间平均（Time Average）” 等于其集总平均。 于是，可以用其任一个样本来得到其数字特征。此性质称为“各态遍历性（Ergodic)”。即：
\(\mu_x=\lim_{M\to\infty}\frac{1}{2M+1}\sum_{n=-M}^{M}x(n)\)
平稳信号不一定是各态遍历的。
最小平方估计
一般来说，希尔伯特空间中线性逼近问题的求解方法称为最小二乘法，通常有三种不同的表现形式：投影法、求导法和匹配法。
:::info
投影法好用啊!一般都用投影法了,简单易理解😝
:::
s
2.经典功率谱估计
概述
平稳随机信号功率谱的两个定义
从自相关函数
\(P_X(e^{j\omega})=\sum_{m=-\infty}^{\infty}r_X(m)e^{-j\omega n}\)
tips:\(r_X(m)\)是对所有样本取了集总平均
这里利用了性质为功率谱密度和\(r_X(m)\)互为傅里叶变换对
直接从定义
对单一无限长的样本求求功率谱:
\(P_x(e^{j\omega})=\lim_{M\to \infty}
\frac{1}{2M+1}|\sum_{n=-M}^{M}x(n)e^{-j\omega n}|^2\)
对于\(P_x(e^{j\omega})\)求集总平均:\(P_X(e^{j\omega})=\mathbb{E}\{P_x(e^{j\omega})\}\)
单一样本不能收敛到所有样本的功率谱, 因此必须有求平均运算.
\(\begin{align}
P_X(e^{j\omega})&=\lim_{M\to \infty}\mathbb{E}\{
\frac{1}{2M+1}|\sum_{n=-M}^{M}x(n)e^{-j\omega n}|^2
\}\\
&=\lim_{M\to \infty}\mathbb{E}\{ \frac{|X(e^{j\omega})|^2}{2M+1}\}
\end{align}\)
取的是随机信号的单个样本\(x(n)\),取极限和求均值计算后变成集总平均
实际情况
实际情况下不能求得所有样本,时间也是有限长.即只能对\(x_M(n)\)进行数据处理, 之间的关系为:
截短后的信号\(x_M(n)\)可以看作功率信号做傅里叶变换, 得到其功率谱
\(P_M(e^{j\omega})=\frac{1}{2M+1}|\sum_{n=-M}^{M}x_M(n)e^{-j\omega n}|\)
用\(x_M(n)\)的功率谱对单个样本或所有样本的功率谱做估计, 效果如何?
🐙两个定义是否等效
\(\begin{aligned}
P_{X}\left(e^{j \omega}\right) &=\lim _{M \rightarrow \infty} E\left\{\frac{1}{2 M+1}\left|\sum_{n=-M}^{M} x_{M}(n) e^{-j \omega n}\right|^2\right\} \\
&=\lim _{M \rightarrow \infty} E\left\{\frac{1}{2 M+1} \sum_{n=-M}^{M} \sum_{m=-M}^{M} x_{M}^{*}(n) x_{M}(m) e^{-j \omega(m-n)}\right\} \\
&=\lim _{M \rightarrow \infty} \frac{1}{2 M+1} \sum_{n=-M}^{M} \sum_{m=-M}^{M} r_{X}(m-n) e^{-j \omega(m-n)}
\end{aligned}\)
先将模平方拆成共轭相乘的形式, 之后将双求和变成单求和
求和变成单求和
\(\sum_{m=-M}^{M} \sum_{m=-M}^{M}g(m-n)=\sum_{k=-2M}^{2M}(2M+1-|k|)g(k)\)
所以实际是一种类似三角窗的情况
| m-n(-2M~2M) | m(-M~M) | n(-M~M) |
|---|---|---|
| -2M | m=-M 1种情况 | n=M 1种情况 |
| -2M+1 | m=-M&m=-M+1 2种 | n=-M+1&n=-M 2种 |
| -2M+2 | 3种情况 | |
| …… | ||
| 0 | M+1种情况 | |
| …… | ||
| 2M | 1种情况 |
所以
\(P_X(e^{j\omega})=\lim_{M\to\infty}\sum_{k=-2M}^{2M}(1-\frac{|k|}{2M+1})r_X(k)e^{-j\omega k}\\
\to^{M\to\infty}\to\sum_{k=-\infty}^{\infty}r_X(k)e^{-j\omega k}\)
所以两种定义相同。
自相关函数的估计
定义
\(r_X(m)=E\{ X^*(n)X(n+m)\}\)集总自相关
\(r_x(m)=\lim_{N\to\infty}\frac{1}{2N+1}\sum_{n=-N}^{N}x^*(n)x(n+m)\)时间自相关（\(x(n)\)为一个无限长的样本）
\(r(m)=\lim_{n\to\infty}\frac{1}{N}\sum_{n=0}^{N-1}x(n)x(n+m)\) \(x(n)\)为实因果信号
实际上能得到的只有一个样本的N个观测值
如何做出估计？
两种方法
- 利用时间自相关的定义直接计算；
2. 先计算出\(r_X(m)\)的能量谱，再对该能量谱作反变换。
🦄直接估计法
\(\hat{r_x}(m)=\frac{1}{N}\sum_{n=0}^{N-1-|m|}x(n)x(n+|m|)\)
从估计方法上看，实际上是把随机信号“视为”单样本有限长的确定性信号。
质量评估
\(bia\{\hat{r}(m)\}=E\{\hat{r}(m)\}-r(m)\)
\(E\{\hat{r}(m)\}=E\{\frac{1}{N}\sum_{n=0}^{N-1-|m|}x(n)x(n+m)\}=\frac{1}{N}\sum_{n=0}^{N-1-|m|}E\{x(n)x(n+m)\}\)
因为有\(r_X(m)=E\{ X^*(n)X(n+m)\}\)，所以\(E\{\hat{r}(m)\}=\frac{1}{N}\sum_{n=0}^{N-1-|m|}r(m)\)（直接把里面的x的求均值看作自相关？（集总平均？））所以\(E\{\hat{r}(m)\}=\frac{N-|m|}{N}r(m)\)
即\(bia\{\hat{r}(m)\}=\frac{-|m|}{N}r(m)\)，偏差是对于真实值的加权
- \(|m|\)固定，\(N\to\infty\),\(bia\{\hat{r}(m)\}=0\)：渐进无偏估计
- \(N\)固定，\(|m|\ll N\),\(\hat{r}(m)\)接近\(r(m)\)，显然取点越多估计精读越高
- \(E\{\hat{r}(m)\}=w(m)r(m),w(m)=\frac{N-|m|}{N}\)（三角窗，来源于有限长的数据可以看作是无限长的数据被加上了矩形窗，两个矩形窗的卷积（自相关）为三角窗）。注意矩形窗加在数据上，三角窗加在相关 函数上，体现在估计的自相关函数的均值上
方差
\(var[\hat{r}(m)]=E\{[\hat{r}(m)-E\{\hat{r}(m)\}]^2\}\)
最后有四阶统计量，太复杂不写了
结论是\(N\to\infty,var[\hat{r}(m)]=0\),渐进一致估计
🦄自相关函数的计算
假设已知单个样本的\(N\)点数据，估计\(\hat{r}_x(m)\)
1.直接按定义
\(\hat{r_x}(m)=\frac{1}{N}\sum_{n=0}^{N-1-|m|}x(n)x(n+|m|)\)
利用性质\(\hat{r}(m)=\hat{r}(-m)\)
✨快速算法
2.利用FFT
1.将\(x_N(n)\)补N个零，得到\(x_{2N}(n)\)
2.对\(x_{2N}(n)\)做FFT得到\(X_{2N}(k)\)
3.对\(X_{2N}(k)\)求幅平方, 的\(|X_{2N}(k)|^2\)
4.得\(\frac{1}{N}|X_{2N}(k)|^2\),做IFFT,得到\(\hat{r}_0(m)\),对其做移位后得到\(\hat{r}_x(m)\)
经典谱估计的基本方法
如何用单一样本的有限长数据去估计原随机信号真实的功率谱
🐱🏍直接法(周期图 Periodogram 法)
\(X_N(\omega)=\sum_{n=0}^{N-1}x(n)e^{-j\omega n},\hat{P}_{PER}(\omega)=\frac{1}{N}|X_N(\omega)|^2\)
\(X_N(k)=\sum_{n=0}^{N-1}x(n)W_N^{nk},\hat{P}_{PER}(\omega)=\frac{1}{N}|X_N(k)|^ 2\)
- 对\(x(n)\)做DTFT（DFT），得到频谱；
- 对该频谱求幅平方，再除以N，即得到“周期图”功率谱，
以此作为对真谱的估计
🐱🏍间接法(Blackman and Tukey, BT法)
\(\hat{r_x}(m)=\frac{1}{N}\sum_{n=0}^{N-1-|m|}x(n)x(n+|m|)\)
\(\hat{P}_{BT}(\omega)=\sum_{m=-M}^{M}\hat{r}_x(m)e^{-j\omega n}, M\leq N-1\)
因为先要估计自相关函数，所以又称间接法
直接法与间接法的关系
关于点数
:::info
\(x_N(n)\)——N点DFT——\(X_N(e^{j\omega})\)————\(\frac{|X_N(e^{j\omega})|}{N}=\hat{P}_{PER}^N(\omega)\)
——补0——\(x_{2N}(n)\)——2N点DFT——\(\frac{|X_{2N}(e^{j\omega})|}{N}=\hat{P}_{PER}^{2N}(\omega)\)
\(x_N(n)\)——xcorr——\(\hat{r}(m)\)——\(\hat r_0(m)\)——2N+1点DFT——\(\hat{P}_{BT}^{N}(\omega)\)
——截短——\(r_M(m)\)——2(M+1)点DFT——\(\hat{P}_{BT}^{2(M+1)}(\omega)\)
:::
经典谱估计的质量
分两种情况讨论：
1.\(M=N-1\)
2.\(M\leq N-1\)
M=N-1
周期图和自相关法是等效的，可统一考虑。此时，周期图和自相关法都是渐近无偏估计；
经典功率谱估计 不是一致估计；
原因：功率谱的定义中既要求极限，又要求均值；而实际的估计方法，仅靠单次实现的有限长，无极限、又无均值运算，因此产生上述问题
\(r(m)\)的傅立叶变换\(P(\omega)\)并不具有各态遍历性。
增大数据长度，效果如何
N 增大,\(D_0(\omega)\)的主瓣B将变窄，因此，引起不相关的区域进一步增多，从而引起谱曲线的更加起伏，实际上是方差变大。通常，增加 N ，会提高谱的分辨率，对经典谱估计来说，增加 N 固然会有利于提高分辨率，但谱曲线的起伏令使用者难以接受，这是经典谱估计的一个固有矛盾。
🦄经典谱估计的改进
主要是改进方差的性能
平滑（Smoothing)
用\(v(m)\)对\(\hat{r}(m)\)加窗实现：\(\hat{P}_{BT}(\omega)=\hat{P}_{PER}(\omega)*V(\omega)\)
平均（Average）
理论依据:** L个独立同分布随机变量和的分布，方差减小L倍**
将一个较长的信号分成若干段，对每一段求功率谱，每一 段的功率谱都是随机变量，然后平均之。类似相干平均，用以 弥补经典谱估计中缺少的求均值运算。
注意：信号应是平稳的，且每一段的统计特性基本一样。
经典谱估计算法比较
3.参数模型功率谱估计
平稳随机信号的参数模型
why使用参数模型？
经典功率谱估计存在先天不足：
- 分辨率低（受窗函数长度的限制）;
- 方差性能不好；
- 方差和分辨率之间的矛盾
参数模型的基本思路
-
假定所研究的平稳过程\(x(n)\)是由一白噪声序列\(u(n)\)激励一线性系统所产生的输出
-
由\(x(n)\)的先验知识，如\(r_x(m)\) ，估计\(H(z)\)的参数；
:::tips -
由\(H(z)\)的参数来估计\(x(n)\)的功率谱\(P(e^{j\omega})=P_u(e^{j\omega})|H(e^{j\omega})|^2=\sigma_u^2\frac{|B(e^{j\omega})|^2}{|A(e^{j\omega})|^2}\)
:::
🦄数学描述
对于LST系统，有两种描述方式：
差分方程：
\(x(n)=-\sum_{k=1}^{p}a_kx(n-k)+\sum_{k=0}^{q}b_ku(n-k)\)
卷积描述：
\(x(n)=\sum_{k=0}^{\infty}h(k)u(n-k)\)
转移函数的两种表示形式：
\(H(z)=\frac{B(z)}{A(z)},B(z)=\sum_{k=0}^{q}b_kz^{-k},A(z)=1+\sum_{k=1}^{p}a_kz^{-k}\)
\(H(z)=\sum_{k=0}^{\infty}h(k)z^{-k}\)
功率谱：
\(P_x(z)=\sigma_u^2H(z)H(z^{-1})\)
AR（Auto-Regressive，自回归）模型
假设\(b_1=b_2=\cdots=b_q=0,b_0=1\)
🦄全极点模型：
现在的输出是现在的输入和过去p个输出的加权和。
易反映谱中的峰值；线性，用的最多，被研究的也最多，性能很好
MA（Moving-Average,移动平均）模型
假设\(a_1=a_2=\cdots=a_p=0,a_0=1\)
易反映谱中的谷值；全零模型，看起来简单；但是非线性；
ARMA（Auto-Regressive Moving-Average，自回归－移动平均）模型
🎶AR模型的正则方程与参数计算
🐵正则方程
目的：找到已知参数和未知参数的关系，以便求解未知参数
未知参数：\(a_1,a_2\cdots a_p,\sigma_u^2\)，共\(p+1\)个；
已知参数：\(r_x(m),m=0,1,\cdots p\)
方程建立：通过差分方程
\(x(n)=-\sum_{k=1}^{p}a_kx(n-k)+u(n)\)
两边同时乘上\(x(n+m)\)，之后求均值：
\(\begin{align}r_x(m)&=\mathbb{E}\{x(n)x(n+m)\}=\mathbb{E}\{-[\sum_{k=1}^{p}a_kx(n+m-k)+u(n+m)]x(n)\}\\
&=-\sum_{k=1}^{p}a_kr_x(m-k)+r_{xu}(m)
\end{align}\)
\(x(n)\)和\(u(n)\)的互相关:
\(r_{xu}(m)=\sigma^2h(-m)=\left\{
\begin{array}[ll]
&0&m\geq1\\
\sigma^2h(0)&m=0
\end{array}
\right.\)
又因为\(h(0)=\lim_{z\to\infty}H(z)=1\)
✨所以最终推出正则方程：
\(r_x(m)=\left\{
\begin{array}[l]
&-\sum_{k=1}^{p}a_kr_x(m-k)&m\geq1\\
-\sum_{k=1}^{p}a_kr_x(k)+\sigma^2&m=0
\end{array}
\right.\)
矩阵表示形式：
$\left[\begin{array}[llllll]
&r_x(0)&r_x(1)&r_x(2)\cdots r_x(p)\
r_x(1)&r_x(0)&r_x(1)\cdots r_x(p-1)\
r_x(2)&r_x(1)&r_x(0)\cdots r_x(p-2)\
\vdots\
r_x(p)&r_x(p-1)&r_x(p-2)\cdots r_x(0)
\end{array}\right]\left[\begin{array}[c]
&1\
a_1\
a_2\
\vdots\
a_p
\end{array}\right]=\left[\begin{array}[c]
&\sigma^2\
0\
0\
\vdots\
0
\end{array}\right]\(<br />🙊(这里利用了性质\)r(m)=r(-m)$）
🙈又称Yule-Walker 方程
🐱👓线性预测问题
问题描述:
设\(x(n)\)在n时刻之前的p个数据：\(x(n-p),x(n-p+1),\cdots,x(n-1)\) 已知，现希望用它们预测\(x(n)\)。
可以利用之前AR模型求出的参数(即假设\(x(n)\)是又某个全极系统产生的?)
\(\hat{x}(n)=-\sum_{k=1}^{p}\alpha_kx(n-k)\\
e(n)=x(n)-\hat{x}(n)\)
均方误差:
\(\rho=\mathbb{E}\{e^2(n)\}\)
所以要求一组最佳的系数\(\alpha_1,\alpha_2\cdots,\alpha_p\),及\(\rho_{min}\)
利用正交原理:\(\mathbb{E}\{\boldsymbol{x}e(n)\}=0,\boldsymbol{x}=[x(n-1),x(n-2),\cdots,x(n-p)]\)
🐳现代数字信号处理第二部分
所以可以得到方程:
:::tips
\(\rho_{min}=\mathbb{E}\{x(n)[x(n)-\hat{x}(n)]\}=r_x(0)+\sum_{k=1}^{p}\alpha_kr_x(k)\\
r_x(m)=-\sum_{k=1}^{p}\alpha_kr_x(m-k),m=1,2,\cdots,p\)
:::
🙈又称Wiener-Hopf 方程
对同一信号\(x(n)\)，都使用其\(r_x(m)\)，得到了两组方程:
Yule-Walk 方程
\(r_x(m)=\left\{
\begin{array}[l]
&-\sum_{k=1}^{p}a_kr_x(m-k)&m\geq1\\
-\sum_{k=1}^{p}a_kr_x(k)+\sigma^2&m=0
\end{array}
\right.\)
Wiener-Hopf 方程
\(r_x(m)=\left\{
\begin{array}[l]
&-\sum_{k=1}^{p}\alpha_kr_x(m-k)&m\geq1\\
-\sum_{k=1}^{p}\alpha_kr_x(k)+\rho_{min}&m=0
\end{array}
\right.\)
对同一信号，二者是相同的,即一个p阶AR模型的系数可以用来构成一个p阶的线性预测器
此外\(e(n)==u(n)\)
分析
- \(\rho_{min}\) 应等于AR模型激励白噪声的功率 \(\sigma^2\) ；
- 线性预测器的误差序列等效于激励AR模型的白噪声序列；
- 由LP的含意，因此AR模型也可以看作是在最小平方意义上对数的拟合；
- 由于LP包含了对数据的外推，因此，对应的谱估计所用数据的范围比实际的应有扩展，因此可以提高分辨率
Levinson-Durbin快速算法
对于AR模型，实际情况下一开始不能立刻把阶数确定下来，（阶数过高可能会将噪声引起的误差也表现在功率谱里）因此希望能够给出一种方法来从低阶到高阶迭代的计算。
利用Toeplitz 矩阵特点，从低阶到高阶进行计算：
需要计算的参数
\(a_m(k);m=1,2,\cdots,p; k=1,2,\cdots,m\)
定义反射系数:\(k_m=a_m(m)\)
\(m=1, a_1(1),\rho_1\)
\(m=2,a_2(1),a_2(2),\rho_2\)
...
\(m=p,a_p(1),a_p(2),\cdots,a_p(p),\rho_\)
第一阶:
\(\left\{\begin{array}{l}
a_{1}(1)=-r_{x}(1) / r_{x}(0) \triangleq k_{1} \\
\rho_{1}=r_{x}(0)-r_{x}^{2}(1) / r_{x}(0)=r_{x}(0)\left[1-a_{1}^{2}(1)\right] \\
\rho_{1}=\rho_{0}\left[1-k_{1}^{2}\right]
\end{array}\right.\)
🐱🐉递推公式
\(\left\{\begin{array}{l} k_{m}=-\left[\sum_{k=1}^{m-1} a_{m-1}(k) r_{x}(m-k)+r_{x}(m)\right] / \rho_{m-1} \\ a_{m}(k)=a_{m-1}(k)+k_{m} a_{m-1}(m-k), k=1,2, \cdots, m-1 \\ \rho_{m}=\rho_{m-1}\left[1-k_{m}^{2}\right] \end{array}\right.\)
:::tips
递推过程中,要始终保持\(|k_m|<1\),即有\(0<\rho_{p}<\rho_{p-1}<\cdots<\rho_{1}<\rho_{0}\)
:::
summary
p阶AR模型（LP）有三组参数：
- p+1个自相关函数:\(r_x(0),\cdots,r_x(p)\)
- p+1个AR模型参数:\(a_p(1),\cdots,a_p(p),\rho_p\)
- p+1个反射系数:\(r_x(0),k_1,\cdots,k_p\)
三组参数可互相导出
AR模型谱估计的性质与阶次选择
AR模型谱估计的性质
:::info
平滑性: AR模型是一有理分式，估计出的谱比经典法平滑。不需要像周期图那样再做平滑或平均，因此，不需要为此去牺牲分辨率
:::
🐱💻外推
AR模型包含了对\(x(n)\)的“预测”或“外推”。实际上，这包含着自相关函数的“外推”
\(P_{AR}(e^{j\omega})=\frac{\rho_p}{|1+\sum_{k=1}^{p}a_ke^{-j\omega k}|^2}=\sum_{\color{red}{m=-\infty}}^{\color{red}\infty}r_a(m)e^{-j\omega m}\)
可以证明
\(r_{a}(m)=\left\{\begin{array}{ll}
r_{x}(m), & |m| \leq p \\
-\sum_{k=1}^{p} a_{k} r_{a}(m-k), & |m|>p
\end{array}\right.\)
而经典谱估计中无外推
\(\hat{P}_{B T}\left(e^{j \omega}\right)=\sum_{\color{red}{m=-p}}^{\color{red}p} \hat{r}_{x}(m) e^{-j \omega m}\)
:::tips
因此AR谱的分辨率要高于经典谱估计
:::
最大熵估计
Burg 最大熵谱估计的思路是:
已知某随机信号自相关函数 \(r_x(m)\)的 p+1个值\(r_x(0),\cdots,r_x(p)\) ，现希望以这 p+1个值对 的自相关函数予以外推。 外推的方法很多，Burg的准则是：外推后的自相关函数对应的时 间序列具有最大的熵，即是最随机的。
:::tips
最大熵法对自相关函数进行外推和AR(p)谱估计的外推性质相同。 自相关函数约束条件下，高斯随机过程的最大熵谱与AR(p)谱完全等效。
:::
AR谱的匹配性质
:::tips
若用AR谱去匹配信号的谱，则误差序列的谱应由常数谱来匹配， 体现 \(A(z)\)的白化性质。
:::
\(\rho=E\left\{e^{2}(n)\right\}=r_{e}(0)=\frac{1}{2 \pi} \int_{-\pi}^{\pi} P_{e}\left(e^{j \omega}\right) d \omega=\frac{1}{2 \pi} \int_{-\pi}^{\pi} P_{x}\left(e^{j \omega}\right)\left|A\left(e^{j \omega}\right)\right|^{2} d \omega\)
且
\(P_{A R}\left(e^{j \omega}\right)=\frac{\rho_{\min }}{\left|A\left(e^{j \omega}\right)\right|^{2}}\)
于是有:
\(\rho=\frac{\rho_{\min }}{2 \pi} \int_{-\pi}^{\pi} \frac{P_{x}\left(e^{j \omega}\right)}{P_{A R}\left(e^{j \omega}\right)} d \omega\)
要令\(\rho\)相对AR模型中的参数\(a_1,\cdots,a_\)最小,给出了yule-walk方程在频域的又一解释
:::info
给定平稳信号 \(x(n)\)的功率谱，希望用一模型的谱\(P_{AR}(e^{j\omeg})\)来匹配它， 匹配的原则是使二者比值的积分最小。
:::
所以说,增加 p，等效地扩大了\(r_a(m),r_x(m)\)相等的部分
阶次选择
理论上阶数越高越接近真实的功率谱,\(p\to\infty\)\(P_{AR}=P\)
:::info
实际上，对自相关函数的估计不可避免地存在误差，大多数AR谱估计方法的标准方程的系数矩阵是满秩的，虽然AR过程实际的阶可能很低，但是解标准方程得到的解，其高阶的值可能是非零的，这些额外的非零参数将产生额外的极点(超过p个)，若额外极点靠近单位圆，就会形成明显的虚假谱峰。 虚假谱峰产生与模型阶选得过高有关， 即相对于观测数据记录长度N来说过高。而过低的阶有可能降低谱估计的分辨率，因此应该正确选取模型的阶
:::
🐱👤AR模型的稳定性
问题:
能否保证\(H(z)=\frac{1}{A(z)}\)是稳定的?
\(A(z)\)是由\(\boldsymbol{R} \boldsymbol{a}=\left[\begin{array}{l} \sigma^{2} \\ \boldsymbol{O}_{p} \end{array}\right]\)解出参数\(a_1,\cdots,a_\),式中自相关函数是估计出的，所以是否稳定取决于自相关函数的性质
:::tips
之前已经证明了自相关函数矩阵式非负定的,即\(\det{(R)}\geq0;\)
:::
:::info
定理1：若\(\mathbf{R}_{p+1}\)正定，则求出的\(a_1\sim a_p\)保证 A(z) 的根都在单位圆内，且唯一。 (AR模型的最小相位性质)
:::
:::info
定理2：若\(x(n)\) 由 p个复正弦组成，即：\(x(n)=\sum_{k=1}^{p}A_ke^{j(\omega_kn+\varphi_n)}\)，则\(\det{(\mathbf{R}_{p+1})}=0;\to\)矩阵奇异
且\(\det{(\mathbf{R}_k)}>0,k=1,\cdots,p\)
:::
:::info
定理3：若\(x(n)\)由 p个正弦组成，\(x(n)\)又称纯谐波过程，则\(x(n)\)是完全可预测的，即可以做到：\(\rho_{min}=0\)
:::
任意一个过程都可以分解成纯谐波过程和其他过程的叠加？
关于信号建模本质的讨论
通过一个白噪声\(u(n)\)激励一个线性系统\(H(z)\)，得到的\(\hat{x}(n)\)并不一定等于\(x(n)\)，即\(\hat{x}(n)\neq x(n)\),也就是说并不能达到时域上的准确建模。实际上我们只做到了自相关函数的匹配和功率谱的匹配：
\(r_{\hat{x}}(m)=r_x(m)\\ P_{\hat{x}}(\omega)=P_x(\omega)\)
统计意义上准确建模
因此，我们讨论过的信号建模是在二阶统计意义上的建模，要求的是自相关函数和功率谱这些二阶统计量的匹配。
:::info
定义：若平稳过程\(x(n)\)存在\(\gamma\)阶模型，使得模型的输出\(\hat{x}(n)\)和\(x(n)\)在\(\gamma\)阶统计意义上一致，则称\(x(n)\)可在\(\gamma\)阶统计意义上准确建模。
:::
阶次大于2的统计分析，称为“高阶谱分析（High-Order Spectral Analysis)”。
👾Wold 分解定理
:::info
任一平稳过程\(x(n)\)都可以分解为一个规则过程\(x_1(n)\)和一个纯正弦过程\(x_2(n)\)的和。二者是不相关的。即\(x(n)=x_1(n)+x_2(n)\);\(x_1(n)\)可表示为一个无穷阶的MA过程; 即\(x_1(n)=1+\sum_{k=1}^{\infty}b_ku(n-k)\)
:::
分解定理也可等效为：任一宽平稳过程的功率谱都可表示为一连续谱和一线谱的和，即
\(P_x(e^{j\omega})=P_{x_1}(e^{j\omega})+P_{x_2}(e^{j\omega})=\sigma^2_u|B(e^{j\omega})|+\sum_{k=1}A^_k\delta(\omega-\omega_k)\)
:::tips
\(x_1(n)\)也可表示为一个有限阶的AR过程： 和MA过程有着相同的单位抽样响应
:::
:::info
引申：一个有限阶的AR模型可近似为一个高阶的MA模型，反之亦然；一个有限阶的ARMA模型也可以用一个阶次足够高的 AR或MA模型来近似。
:::
🚴♂️关于线性预测的进一步讨论
上一节使用的AR模型等效于一个p阶的线性预测器。即 Yule-Walker方程等效于Wiener-Hopf 方程。但估计的功率谱的分辨率不理想，其原因是仅用了前向预测，
👺双向预测
\(\hat{x}^f(n)=-\sum_{k=1}^{p}a^f_kx(n-k)前向预测\\
\hat{x}^b(n)=-\sum_{k=1}^{p}a^b_kx(n+k)后向预测\\\)
对同一组数据:\(\hat{x}^b(n-p)=-\sum_{k=1}^{p}a_k^bx(n-p+k)\)
\(e^b(n-p)=x(n-p)-\hat{x}^b(n-p) \buildrel \Delta \over =e^b(n)\)
后向预测误差功率:
\(\rho^b=\mathbb{E}\{e^b(n)\}\)
使用正交原理可以得:
:::tips
后向预测的Wiener-Hopf Eq.
\(\left\{\begin{align}
\rho_{\min }^{b}&=r_{x}(0)+\sum_{k=1}^{p} a^{b}(k) r_{x}(k) \\
r_{x}(m)&=-\sum_{k=1}^{p} a^{b}(k) r_{x}(m-k), m=1,2, \ldots, p
\end{align}\right.\)
:::
所以有:
\(\begin{aligned}
\rho_{\min }^{b} &=\rho_{\min }^{f} \\
a^{b}(k) &=a^{f}(k)
\end{aligned}\)
上述结果表明，使用已知的 p 个数据，我们可以实现前向预测，也可以实现后向预测，两种情况下可各自得到对等的**WienerHopf **方程。
:::info
将它们单独使用，所得分辨率都不理想。可以设想，如将二者结合起来，即同时使前向、后向预测误差功率为最小，应能得到更好的分辨率。
:::
lattice结构*
前、后向预测误差序列有如下的关系：
:::info
\(\left.\begin{array}{rl}
e_{m}^{f}(n) & =e_{m-1}^{f}(n)+k_{m} e_{m-1}^{b}(n-1) \\
e_{m}^{b}(n) & =e_{m-1}^{b}(n-1)+k_{m}^{*} e_{m-1}^{f}(n)
\end{array}\right\} \quad m=1,2, \cdots, p\)
初始条件:
\(e^f_0(n)=e^b_0(n)=x(n)\)
反射系数
\(k_{m}=\frac{-<e_{m-1}^{b}(n-1), e_{m-1}^{f}(n)>}{\left\|e_{m-1}^{f}(n-1)\right\|^{*}\left\|e_{m-1}^{b}(n-1)\right\|}=-\frac{\operatorname{cov}\left(e_{m-1}^{b}(n-1), e_{m-1}^{f}(n)\right)}{\sqrt{\operatorname{var}\left(e_{m-1}^{f}(n)\right)} \sqrt{\operatorname{var}\left(e_{m-1}^{b}(n-1)\right)}}\)
:::
证明略
:::tips
上述关系引出了线性预测中的Lattice结构。这一结构在现代谱 估计、语音信号处理中有着重要的应用
:::
👅关于实际数据加窗
之前所求的关系都是集总平均。对实际的信号：单个样本有限长，求均值要简化,即:
\(\rho^f=\sum_{n}|e^f(n)|^2\buildrel 替代\over\to \mathbb{E}\{e^f(n)\}\\
\rho^b=\sum_{n}|e^b(n)|^2\buildrel 替代\over\to \mathbb{E}\{e^b(n)\}\)
问题在于: n的取值范围?
\(e^{f}(n)=x(n)-\hat{x}^{f}(n)=x(n)+\sum_{k=1}^{p} a^{f}(k) x(n-k)\)
做前向预测时,n的取值可以时\(0\sim N-1+p\),于是有:
记为\(\mathbf{X}_{0 (N+p)\times (p+1)}\)
此外,记:
\(\left.\begin{array}{l}
\boldsymbol{X}_{1}=\left[\begin{array}{ccc}
x(p) & \cdots & x(0) \\
\vdots & \vdots & \vdots \\
x(N-1) & \cdots & x(N-1-p)
\end{array}\right]_{(N-p) \times(p+1)} \\
\boldsymbol{X}_{2}=\left[\begin{array}{ccc}
x(0) & & \\
\vdots & \ddots & \\
x(p) & \cdots & x(0) \\
\vdots & \vdots & \vdots \\
x(N-1) & \cdots & x(N-1-p)
\end{array}\right]_{N \times(p+1)} \\
\boldsymbol{X}_{3}=\left[\begin{array}{ccc}
x(p) & \cdots & x(0) \\
\vdots & \vdots & \vdots \\
x(N-1) & \cdots & x(N-1-p) \\
& \ddots & \vdots \\
& & x(N-1)
\end{array}\right]_{N \times(P+1)}
\end{array}\right\}\)
:::info
使用\(X_0 、X_2 、X_3\)的算法，用来估计预测误差功率所使用的 预测误差值涉及到了n<0和n>N-1范围内的观测数据，实际上引入并假定这些值为0，即对随机过程的无限长取样序列加了窗。 三角区域对应的预测误差使用了很多假定为0的数据，因此有很大的误差，所以虽然最大限度地使用了预测误差值，反而效果不好。因此后面提到的算法中，使用了\(X_1\)的协方差法或者 Marple算法的谱估计的分辨率比使用了\(X_0\)的自相关法高。
:::
AR模型参数求解算法
AR模型系数求解算法很多，人们目前仍在探讨新的求解算法。 目前，常用的算法是：
-
自相关法
-
Burg算法
-
改进的协方差算法（modified covariance)
:::info
各算法之间的主要区别： -
\(e^f(n)\),n的取值范围
-
仅用前向预测，还是前后向都预测,即\(\rho^f\)最小还是\(\rho^f+\rho^\)
-
递推算法: 是由\(x(n)\)求\(r_x(m)\),然后由\(r_x(m)\)递推还是直接由\(x(n)\)递推?
:::
自相关法
- 只用前向预测\(e^f(n)\)
- 使用\(X_0\),前后加窗，分辨率不好。
\(\boldsymbol{e}^f_p=[e_p^f(0),e_p^f(1),\cdots,e_p^f(N-1+p)]\\
a_p^f=[a^f(1),a^f(2),\cdots,a^f(p)]\)
目标:
\(\min: \rho^f=\sum_{n=0}^{N-1+p}|e_p^f(n)|^2=(\boldsymbol{e}^f_p)^H\boldsymbol{e}^f_p\)
还是最小误差均方估计,实际上等效于Yule-Walker方程,用Levinson递推
Burg算法
- \(\rho^{fb}=\frac{1}{2}(\rho^f+\rho^b)\)
- 前后都不加窗(取\(X_1\))
令\(\frac{\rho^{fb}}{\partial k_m}=0\),求\(k_\)
- 用递推公式直接对数据递推
:::info
递推公式:
- 先求\(\hat{k}_{m}=\frac{-2 \sum_{n=m}^{N-1} e_{m-1}^{f}(n) e_{m-1}^{b^{*}}(n-1)}{\sum_{n=m}^{N-1}\left|e_{m-1}^{f}(n)\right|^{2}+\sum_{n=m}^{N-1}\left|e_{m-1}^{b}(n-1)\right|^{2}}=a_m(m)\)
- 首先令\(e_0^f(n)=e_0^b(n)=x(n)\),求\(k_1\)
- 然后求\(\rho_1,a_1(1)\):\(\left\{\begin{array}{l} \hat{a}_{m}(k)=\hat{a}_{m-1}(k)+\hat{k}_{m} \hat{a}_{m-1}^{*}(m-k) \\ \hat{a}_{m}(m)=\hat{k}_{m}, k=1,2, \cdots, m-1 \\ \hat{\rho}_{m}=\left(1-\left|\hat{k}_{m}\right|^{2}\right) \hat{\rho}_{m-1} \end{array}\right.\)
- 根据递推关系求\(e_1^f(n),e_1^b(n)\):\(\left\{\begin{array}{l} e_{m}^{f}(n)=e_{m-1}^{f}(n)+k_{m} e_{m-1}^{b}(n-1) \\ e_{m}^{b}(n)=e_{m-1}^{b}(n-1)+k_{m}^{*} e_{m-1}^{f}(n), \quad m=1,2, \cdots, p \\ e_{0}^{f}(n)=e_{0}^{b}(n)=x(n) \end{array}\right.\)
- 求\(k_2\)
- 类似step3,求系数;
- 重复上述过程;
:::
:::tips
summary:
Burg法步骤一般就是先求前后项误差,然后根据误差求反射系数\(k_m\),然后根据反射系数算本阶的参数系数;
:::
改进的协方差法(Maple法)
和burg法基本一样,唯一区别在于
令\(\frac{\partial\rho^{fb}}{\partial a_m(i)}=0;(i=1\sim m; m=1\sim p)\)
上述最小化的结果是得到一个协方差方程;该矩阵不是Toeplitz矩阵，因此不能用Levinson算法求解
\(\left[\begin{array}{lllll}
c_{x}(1,1) &c_{x}(1,2)& \cdots & c_{x}(1, p) \\
c_{x}(2,1) &c_{x}(2,2) &\cdots & c_{x}(2, p) \\
\vdots \\
c_{x}(p, 1)& c_{x}(p, 2) &\cdots &c_{x}(p, p)
\end{array}\right]\left[\begin{array}{l}
\widehat{a}(1) \\
\hat{a}(2) \\
\vdots \\
\hat{a}(p)
\end{array}\right]=-\left[\begin{array}{l}
c_{x}(1,0) \\
c_{x}(2,0) \\
\vdots \\
c_{x}(p, 0)
\end{array}\right]\)
MA 模型及其正则方程
MA模型是一个全零点模型
\(\begin{align}
x(n)&=u(n)+\sum_{k=1}^{q}b_ku(n-k)\\
H(z)&=B(z)=1+\sum_{k=1}^{q}b_kz^{-k}\\
P(e^{j\omega})&=\sigma_u^2|B(e^{j\omega})|^2
\end{align}\)
类似于AR模型, 求相关函数:
$\begin{align}
r_x(m)&=\mathbb{E}{x(n)x(n+m)}\
&=\sum_{k=0}^{q}b_kr_{xu}(m-k)\
&=\left{
\begin{array}
&\sigma2\sum_{k=0}b(k)b(k+m)&m=0,1,2,\cdots,q\
0&m>q
\end{array}
\right.
\end{align}$
正则方程是一个非线性方程组
:::info
由于MA模型的正则方程是非线性方程，所以人们 提出了很多的求解算法，如谱分解、基于迭代的方法、基于高阶AR模型近似的方法。后者最好用，基础是Wold分解定理。
:::
用高阶的AR模型来近似MA模型
:::tips
Wold分解定理：任一平稳过程可分解为一规则过程和一纯正弦过程（连续谱和线谱的和）
:::
\(H_{\infty}=\frac{1}{A(z)}=\frac{1}{1+\sum_{k=1}^{\infty}a(k)z^{-k}}=B(z)=1+\sum_{k=1}^{q}b(k)z^{-k}\)
于是有\(A_{\infty}(z)B(z)=1\)
\(\begin{align}
a(k)*b(k)&=\delta(k)\\
a(m)+\sum_{k=1}^{q}b(k)a(m-k)&=\delta(m)
\end{align}\)
观察上式，发现形式上和线性预测的误差公式一样！
:::info
线性预测:
\(\hat{x}(n)=-\sum_{k=1}^{p}\alpha_kx(n-k)\)
\(e(n)=x(n)+\sum_{k=1}^{p}\alpha_kx(n-k)\)
:::
如果把\(a_p(k)\)看作用来预测的\(x(n)\),\(b(k)\)看作AR模型的参数,那么可以用AR模型的参数来在求一组AR模型的参数, 来对MA模型做近似, 这组参数即为MA模型的参数\(b(k)\)
即
\(\hat{a}(m)+\sum_{k=1}^{q}b(k)\hat{a}(m-k)=e(m)\)
:::info
与之前推导的式子的区别只有等式右侧为\(e(m)\)而非\(\delta(m)\),(只要令\(e(0)=1,e(m)=0;m\neq0)\)即可,实际上\(e(m)\)由有限长数据计算出,不可能为0;因此只需要令\(\hat{\rho}_{MA}=\sum_{m}|e(m)|^2\)最小
:::
可以看出这即为AR模型.
SUMMARY
:::tips
- 由\(x(n),n=0,1,\cdots,N-1,\)建立p阶AR模型(\(p>>q\)),得到\(\hat{a}_p(k)\)
- 对\(\hat{a}_p(k)\)建立q阶AR模型求得系数\(b(k)\)
:::
浙公网安备 33010602011771号