\[
\newcommand{\bs}{\boldsymbol}
\newcommand{\bsX}{\boldsymbol{X}}
\newcommand{\bf}{\mathbf}
\newcommand{\msc}{\mathscr}
\newcommand{\mca}{\mathcal}
\newcommand{\T}{\text{T}}
\newcommand{\rme}{\mathrm{e}}
\newcommand{\rmi}{\mathrm{i}}
\newcommand{\rmj}{\mathrm{j}}
\newcommand{\rmd}{\mathrm{d}}
\newcommand{\rmm}{\mathrm{m}}
\newcommand{\rmb}{\mathrm{b}}
\newcommand{\and}{\land}
\newcommand{\or}{\lor}
\newcommand{\exist}{\exists}
\newcommand{\sube}{\subseteq}
\newcommand{\lr}[3]{\left#1 #2 \right#3}
\newcommand{\intfy}{\int_{-\infty}^{+\infty}}
\newcommand{\sumfy}[1]{\sum_{#1=-\infty}^{+\infty}}
\newcommand{\vt}{\vartheta}
\newcommand{\ve}{\varepsilon}
\newcommand{\vp}{\varphi}
\newcommand{\Var}{\text{Var}}
\newcommand{\Cov}{\text{Cov}}
\newcommand{\edef}{\xlongequal{def}}
\newcommand{\prob}{\text{P}}
\newcommand{\Exp}{\text{E}}
\newcommand{\t}[1]{\text#1}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newcommand{\versionofnewcommand}{\text{260125}}
\]
Linear Prediction Coding (LPC) 线性预测编码
本节主要强调平稳性假设. 实际上, 在随机过程中, 平稳性是非常重要的一点. 因此, 我们也从随机过程开始.
假设离散的随机序列 \(X_0,\cdots, X_n,\cdots\) 具有平稳性, 即它的某种形式在统计意义上有时间不变性.
本节我们把注意力放在 "宽平稳, Wide Sense Stationary", 即随机序列的一阶距和二阶距具有时间不变性.
\[
E(X_k)=m_k\equiv m\\
R_{ij}=E(X_i X_j)=R_{|i-j|}
\]
其中, 第二点即是说, 随机序列中两项的相关值 (或者根据上一节, 称为内积) 只依赖于两项的距离, 而与它们在序列中的具体位置无关.
Toeplitz 矩阵
在生活中, 我们常常做的事是用随机序列中已经发生的一串信息, 估计未来下一时刻发生的信息.
\[
\begin{aligned}
&&\boldsymbol X&=(X_1,\cdots,X_{n-1})^\text T\to X_n\\
\Rightarrow& &X_n&=\boldsymbol \alpha^\text T X\\
\Rightarrow& &\boldsymbol\alpha_{opt}&=(E(\boldsymbol{XX}^\text{T}))^{-1}E(\boldsymbol{X}X_n)
\end{aligned}
\]
其中, 分母这个矩阵的逆是很难算的. 然而, 平稳性竟可以让我们不必算这个逆, 以更快的速度解出 \(\boldsymbol \alpha\)
\[
E(\boldsymbol{XX^\text T})=
{
\begin{bmatrix}
R_{0,0}&\cdots&R_{0,n-1}\\
\vdots &\ddots&\vdots \\
R_{n-1,0} &\cdots &R_{n-1,n-1}
\end{bmatrix}
}
\overset{stationary}=
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1}\\
R_1 &\ddots &\ddots &\vdots \\
\vdots &\ddots &\ddots &R_1 \\
R_{n-1} &\cdots &R_1 &R_{0}
\end{bmatrix}
}
\]
矩阵的自由度原本为 \(n^2\) , 现在神奇地竟变为 \(n\) 了. 我们称右边这样的矩阵为 Toeplitz 矩阵.
再看另一项,
\[
E(\boldsymbol X\cdot X_n)=
\begin{bmatrix}
R_{0,n}\\
\vdots\\
R_{n-1,n}
\end{bmatrix}\overset{stationary}=
\begin{bmatrix}
R_{n}\\
\vdots\\
R_{1}
\end{bmatrix}
\]
这样就有:
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1}\\
R_1 &\ddots &\ddots &\vdots \\
\vdots &\ddots &\ddots &R_1 \\
R_{n-1} &\cdots &R_1 &R_{0}
\end{bmatrix}
}
\cdot
\begin{bmatrix}
\alpha_{1}\\
\vdots\\
\alpha_{n}
\end{bmatrix}=
\begin{bmatrix}
R_{n}\\
\vdots\\
R_{1}
\end{bmatrix}
\]
Levinson-Durbin 迭代
想解 \(\alpha_1,\cdots,\alpha_n\), 首先, 齐次化
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1
\end{bmatrix}
}
\cdot
{
\begin{bmatrix}
\alpha_{1}\\
\vdots\\
\alpha_{n}\\
-1
\end{bmatrix}
}=\boldsymbol{0}
\]
不过这样, 左边这个大矩阵就是 \(n\times (n+1)\) 的了, 直觉上, 我们自然希望补一行出来, 同时让矩阵仍保持原本的模样.
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}
\cdot
{
\begin{bmatrix}
\alpha_{1}\\
\vdots\\
\alpha_{n}\\
-1\\
\end{bmatrix}
}=
{
\begin{bmatrix}
0\\
\vdots\\
0\\
P_n
\end{bmatrix}
}
\]
其中, \(P_n=-R_0+\sum_{k=0}^{n-1} (\alpha_{k+1}R_{n-k}-R_0)\), 尽管它具体长啥样其实不太重要.
现在关键的问题是求 Toeplitz 矩阵的逆. 接下来, 采用递推的思维: 用低一阶的结果递推高一阶的结果. 这里, 低一阶 be like
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} \\
R_1 &\ddots &\ddots &\vdots \\
\vdots &\ddots &\ddots &R_1 \\
R_{n-1} &\cdots &R_1 &R_{0} \\
\end{bmatrix}
}
\cdot
{
\begin{bmatrix}
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1\\
\end{bmatrix}
}=
{
\begin{bmatrix}
0\\
\vdots\\
0\\
P_{n-1}
\end{bmatrix}
}
\]
递推的思维里, \(\alpha_1^{(n-1)},\cdots,\alpha_n^{(n-1)}\) 都是已知的. 现在需要用这些低阶的解, 构造出高阶的情况.策略是这样的:
\[
{
\begin{bmatrix}
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1\\
\end{bmatrix}
}\to
{
\begin{bmatrix}
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1\\
0\\
\end{bmatrix}
}\\
\Rightarrow
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}\cdot
{
\begin{bmatrix}
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1\\
0\\
\end{bmatrix}
}=
{
\begin{bmatrix}
0\\
\vdots\\
0\\
P_{n-1}\\
Q_n\\
\end{bmatrix}
}
\]
其中, 左上角的 \(n\times n\) 矩阵是低阶 Toeplitz.
接下来想办法用 Toeplitz 带来的性质: 假设有个 Toeplitz 矩阵 \(A\), \(A_{ij}\) 的值仅依赖于 \(i-j\). 于是:
\[
A
{
\begin{bmatrix}
x_{1}\\ \vdots\\ x_{n}
\end{bmatrix}
}=
{
\begin{bmatrix}
y_{1}\\ \vdots\\ y_{n}
\end{bmatrix}
}\Rightarrow
A
{
\begin{bmatrix}
x_{n}\\ \vdots\\ x_{1}
\end{bmatrix}
}=
{
\begin{bmatrix}
y_{n}\\ \vdots\\ y_{1}
\end{bmatrix}
}
\]
嗯 ... ... 看起来颠倒两个向量以后, 局势也没有变得更好. 不过对上述构造进行一个小小的更改便好. 我们不将左上角看作低阶 Toeplitz, 而是右下角! 对于 Toeplitz 矩阵, 这是完全可行的! 那么就有:
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}\cdot
{
\begin{bmatrix}
0\\
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1
\end{bmatrix}
}=
{
\begin{bmatrix}
Q_n\\
0\\
\vdots\\
0\\
P_{n-1}
\end{bmatrix}
}
\]
这看起来, 两个列向量就有点对称的意味了. 上下颠倒罢!
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}\cdot
{
\begin{bmatrix}
-1\\
\alpha_{n}^{(n-1)}\\
\vdots\\
\alpha_{1}^{(n-1)}\\
0
\end{bmatrix}
}=
{
\begin{bmatrix}
P_{n-1}\\
0\\
\vdots\\
0\\
Q_n
\end{bmatrix}
}
\]
有没有一种冲动, 想把它们加起来. 不过直接相加得不到太好的结果, 我们可以乘个待定系数 \(\lambda\) , 引入一个新的自由度:
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}\cdot
\left(
{
\begin{bmatrix}
0\\
\alpha_{1}^{(n-1)}\\
\vdots\\
\alpha_{n}^{(n-1)}\\
-1
\end{bmatrix}
}+\lambda
{
\begin{bmatrix}
-1\\
\alpha_{n}^{(n-1)}\\
\vdots\\
\alpha_{1}^{(n-1)}\\
0
\end{bmatrix}
}\right)=
\left(
{
\begin{bmatrix}
Q_n\\
0\\
\vdots\\
0\\
P_{n-1}
\end{bmatrix}
}+\lambda
{
\begin{bmatrix}
P_{n-1}\\
0\\
\vdots\\
0\\
Q_n
\end{bmatrix}
}
\right)
\]
我们可以让 \(Q_n+\lambda P_{n-1}=0\), 从而构造出高一阶中符合规范的形式. 而且, 关键的是, \(Q_n\) 和 \(P_{n-1}\) 都仅仅依赖于已知的\(\alpha^{(n-1)}\).
于是有:
\[
{
\begin{bmatrix}
R_{0} &R_1 &\cdots &R_{n-1} &R_n\\
R_1 &\ddots &\ddots &\vdots &R_{n-1} \\
\vdots &\ddots &\ddots &R_1 &\vdots \\
R_{n-1} &\cdots &R_1 &R_{0} &R_1\\
R_n &R_{n-1}&\cdots &R_1 &R_0
\end{bmatrix}
}\cdot
{
\begin{bmatrix}
-\lambda\\
\alpha_1^{(n-1)}+\lambda\alpha_{n}^{(n-1)}\\
\vdots\\
\alpha_n^{(n-1)}+\lambda\alpha_{1}^{(n-1)}\\
-1
\end{bmatrix}
}=
{
\begin{bmatrix}
0\\
0\\
\vdots\\
0\\
P_{n-1}+\lambda Q_n
\end{bmatrix}
}
\]
因此有:
\[
\begin{aligned}
&\lambda=-Q_n/P_{n-1}\\
\\
&\alpha^{(n)}_{k}=
\begin{cases}
-\lambda, &k=1\\
\alpha_{k-1}^{(n-1)}+\lambda\alpha_{n-k+2}^{(n-1)},&k=2,\cdots,n+1
\end{cases}\\
\\
&P_n=P_{n-1}+\lambda Q_n=(1-\lambda^2)P_{n-1}\\
\end{aligned}
\]
这显然就太好解了.
该过程称为 Levinson-Durbin 迭代. 在线性预测编码 (LPC) 中, 我们不需要传递原始信号, 而只需要传递这些简单易算的参数就好, 大大降低了带宽. 该过程在音频处理中是非常核心的.
从滤波器的角度看 LPC
考虑一个 \(X_n\overset{H}\longrightarrow Y_n\) 的系统, 如果希望 \(H\) 起到 LPC 的作用, 该如何办呢? 首先, 考虑一下 \(H\) 的基本架构, 有限冲激响应 FIR. 实际上, 这就是个很自然的递推过程. 做一个延时单元是一阶, 再做一个成两阶, 再做一个成三阶, 如此往复. 不过, 之前讨论中的 "翻转" 从何体现?
我们换个角度来看 "翻转": 从 Lattice 架构的角度.
tbc...