滤波理论
滤波理论的核心问题：利用含噪观测信号估计期望信号，使估计误差尽可能小。
观测模型：
\[x(n)=s(n)+v(n)\]
其中：
- \(x(n)\)：观测信号
- \(s(n)\)：期望信号
- \(v(n)\)：噪声
滤波器输出：
\[y(n)=\hat{s}(n)\]
误差：
\[e(n)=s(n)-\hat{s}(n)\]
最优准则：
\[\min E[e^2(n)]\]
也就是使均方误差最小。
维纳滤波和卡尔曼滤波的对比
共同点：
- 都用于从含噪观测中估计信号。
- 都属于最佳线性估计方法。
- 都以均方误差最小为准则。
- 在平稳条件下，卡尔曼滤波的稳态结果与维纳滤波一致。
维纳滤波：
- 根据当前和过去的观测值估计当前信号。
- 解通常以滤波器系统函数 \(H(z)\) 或冲激响应 \(h(n)\) 给出。
- 本质是最佳线性滤波器。
- 主要适用于平稳随机过程。
- 需要已知信号与噪声的二阶统计量，如自相关函数、互相关函数、功率谱密度。
- 更适合标量滤波问题。
- 通常是非递推形式，需要解维纳-霍夫方程。
卡尔曼滤波：
- 用前一个状态估计值和当前观测值递推估计当前状态。
- 解通常以状态估计值形式给出。
- 本质是线性最优估计器。
- 适用于平稳和非平稳随机过程。
- 需要已知状态方程、量测方程以及噪声二阶统计量。
- 适合矢量状态估计问题。
- 是递推算法，适合实时处理。
维纳滤波典型形式：
\[\hat{s}(n)=\sum_{k}h(k)x(n-k)\]
卡尔曼滤波典型形式：
\[\hat{\mathbf x}_{n|n}=\hat{\mathbf x}_{n|n-1}+\mathbf K_n\left[\mathbf y_n-\mathbf H_n\hat{\mathbf x}_{n|n-1}\right]\]
总结：
- 维纳滤波：已知统计特性，直接设计一个最优线性系统。
- 卡尔曼滤波：已知状态空间模型，递推更新状态估计。
- 维纳滤波偏“滤波器设计”。
- 卡尔曼滤波偏“状态估计递推”。
离散时间维纳滤波
设输入观测信号为 \(x(n)\)，期望信号为 \(d(n)\)，滤波器冲激响应为 \(h(n)\)。
滤波器输出：
\[y(n)=\sum_{k=0}^{\infty}h(k)x(n-k)\]
估计误差：
\[e(n)=d(n)-y(n)\]
均方误差：
\[J=E[e^2(n)]\]
优化目标：
\[\min_h E[e^2(n)]\]
时域解（维纳霍夫方程）
因果维纳滤波器输出：
\[y(n)=\sum_{k=0}^{\infty}h(k)x(n-k)\]
误差：
\[e(n)=d(n)-\sum_{k=0}^{\infty}h(k)x(n-k)\]
最优条件：估计误差与所有用于估计的观测样本正交。
\[E[e(n)x(n-m)]=0,\qquad m=0,1,2,\cdots\]
代入误差表达式：
\[E\left[\left(d(n)-\sum_{k=0}^{\infty}h(k)x(n-k)\right)x(n-m)\right]=0\]
得到：
\[E[d(n)x(n-m)]=\sum_{k=0}^{\infty}h(k)E[x(n-k)x(n-m)]\]
定义互相关函数：
\[r_{dx}(m)=E[d(n)x(n-m)]\]
定义自相关函数：
\[r_{xx}(m-k)=E[x(n-k)x(n-m)]\]
维纳-霍夫方程：
\[r_{dx}(m)=\sum_{k=0}^{\infty}h(k)r_{xx}(m-k),\qquad m=0,1,2,\cdots\]
矩阵形式：
\[\mathbf R_{xx}\mathbf h=\mathbf r_{dx}\]
若只取有限阶 \(M\) 阶 FIR 滤波器：
\[y(n)=\sum_{k=0}^{M-1}h(k)x(n-k)\]
则：
\[\begin{bmatrix}r_{xx}(0) & r_{xx}(1) & \cdots & r_{xx}(M-1)\\r_{xx}(1) & r_{xx}(0) & \cdots & r_{xx}(M-2)\\\vdots & \vdots & \ddots & \vdots\\r_{xx}(M-1) & r_{xx}(M-2) & \cdots & r_{xx}(0)\end{bmatrix}\begin{bmatrix}h(0)\\h(1)\\\vdots\\h(M-1)\end{bmatrix}=\begin{bmatrix}r_{dx}(0)\\r_{dx}(1)\\\vdots\\r_{dx}(M-1)\end{bmatrix}\]
解为：
\[\mathbf h=\mathbf R_{xx}^{-1}\mathbf r_{dx}\]
最小均方误差：
\[J_{\min}=E[d^2(n)]-\sum_{k=0}^{\infty}h(k)r_{dx}(k)\]
有限维形式：
\[J_{\min}=r_{dd}(0)-\mathbf r_{dx}^T\mathbf R_{xx}^{-1}\mathbf r_{dx}\]
结论：
- 维纳-霍夫方程来自正交性原理。
- 求维纳滤波器就是求解一组线性方程。
- 时域解直观，但无限阶因果滤波器求解困难。
z域解
时域求解维纳-霍夫方程较困难，因此可转到 \(z\) 域求解。
观测信号功率谱分解：
\[S_{xx}(z)=\sigma_w^2B(z)B(z^{-1})\]
其中：
- \(B(z)\) 是因果、稳定、最小相位系统。
- \(w(n)\) 是白噪声。
- \(B(z)\) 可看作把白噪声 \(w(n)\) 变成观测信号 \(x(n)\) 的成形滤波器。
即：
\[x(n)=b(n)*w(n)\]
白化滤波器为：
\[\frac{1}{B(z)}\]
将观测信号白化：
\[w(n)=x(n)*b^{-1}(n)\]
设对白化信号的最优滤波器为 \(G(z)\)，则原系统滤波器为：
\[H(z)=\frac{G(z)}{B(z)}\]
因此问题转化为求 \(G(z)\)。
非因果维纳滤波器：
\[G_{\text{opt}}(z)=\frac{S_{ws}(z)}{\sigma_w^2}\]
由于：
\[S_{xs}(z)=B(z^{-1})S_{ws}(z)\]
所以：
\[S_{ws}(z)=\frac{S_{xs}(z)}{B(z^{-1})}\]
非因果维纳滤波器：
\[H_{\text{opt}}(z)=\frac{S_{xs}(z)}{S_{xx}(z)}\]
若 \(x(n)=s(n)+v(n)\)，且 \(s(n)\) 与 \(v(n)\) 不相关，则：
\[S_{xs}(z)=S_{ss}(z)\]
\[S_{xx}(z)=S_{ss}(z)+S_{vv}(z)\]
因此：
\[H_{\text{opt}}(z)=\frac{S_{ss}(z)}{S_{ss}(z)+S_{vv}(z)}\]
这是非因果维纳滤波器的频域形式。
因果维纳滤波器：
由于要求因果性：
\[h(n)=0,\qquad n<0\]
因此需要取因果部分。
因果解为：
\[G_{\text{opt}}(z)=\frac{1}{\sigma_w^2}\left[\frac{S_{xs}(z)}{B(z^{-1})}\right]_+\]
其中 \([\cdot]_+\) 表示取因果部分，即保留 \(z\) 反变换中 \(n\ge 0\) 的项。
原滤波器：
\[H_{\text{opt}}(z)=\frac{G_{\text{opt}}(z)}{B(z)}\]
所以：
\[H_{\text{opt}}(z)=\frac{1}{\sigma_w^2B(z)}\left[\frac{S_{xs}(z)}{B(z^{-1})}\right]_+\]
因果维纳滤波器设计步骤：
- 对观测信号功率谱做谱因式分解：
\[S_{xx}(z)=\sigma_w^2B(z)B(z^{-1})\]
- 构造：
\[\frac{S_{xs}(z)}{B(z^{-1})}\]
- 做 \(z\) 反变换并取因果部分：
\[\left[\frac{S_{xs}(z)}{B(z^{-1})}\right]_+\]
- 得到：
\[H_{\text{opt}}(z)=\frac{1}{\sigma_w^2B(z)}\left[\frac{S_{xs}(z)}{B(z^{-1})}\right]_+\]
结论：
\[H_{\text{opt}}(z)=\frac{S_{xs}(z)}{S_{xx}(z)}\]
- 因果维纳滤波器需要谱因式分解和取因果部分。
- \(z\) 域解的关键是白化观测信号。
维纳预测
预测问题：已知过去观测值，估计当前或未来信号值。
已知：
\[x(n-1),x(n-2),\cdots,x(n-P)\]
估计：
\[\hat{s}(n+N),\qquad N\ge 0\]
其中：
- \(N=0\)：一步当前预测。
- \(N>0\)：未来预测。
预测器形式：
\[\hat{s}(n+N)=\sum_{k=1}^{P}a_kx(n-k)\]
误差：
\[e(n)=s(n+N)-\hat{s}(n+N)\]
均方误差：
\[J=E[e^2(n)]\]
优化目标：
\[\min_{a_k}E[e^2(n)]\]
正交性原理：
\[E[e(n)x(n-m)]=0,\qquad m=1,2,\cdots,P\]
代入：
\[E\left[\left(s(n+N)-\sum_{k=1}^{P}a_kx(n-k)\right)x(n-m)\right]=0\]
得到预测的维纳-霍夫方程：
\[r_{sx}(N+m)=\sum_{k=1}^{P}a_kr_{xx}(m-k),\qquad m=1,2,\cdots,P\]
矩阵形式：
\[\mathbf R_{xx}\mathbf a=\mathbf r_{sx}\]
解为：
\[\mathbf a=\mathbf R_{xx}^{-1}\mathbf r_{sx}\]
最小预测误差：
\[J_{\min}=r_{ss}(0)-\mathbf r_{sx}^T\mathbf R_{xx}^{-1}\mathbf r_{sx}\]
预测能够实现的原因：
- 信号内部存在相关性。
- 平稳随机信号的自相关函数只与时间间隔有关。
- 数据间相关性越强，预测越准确。
- 白噪声前后样本不相关，因此无法预测。
- 预测越远，相关性越弱，误差通常越大。
结论：
- 维纳预测本质上仍是最小均方误差线性估计。
- 与维纳滤波的区别在于：滤波估计当前信号，预测估计当前或未来信号。
- 预测性能由信号的相关性决定。
连续信号维纳滤波
连续时间观测信号：
\[x(t)=s(t)+v(t)\]
线性滤波器输出：
\[y(t)=\int_{0}^{\infty}h(\alpha)x(t-\alpha)d\alpha\]
误差：
\[e(t)=s(t)-y(t)\]
优化目标：
\[\min_h E[e^2(t)]\]
由正交性原理：
\[E[e(t)x(t-\tau)]=0,\qquad 0\le \tau<\infty\]
代入：
\[E\left[\left(s(t)-\int_0^\infty h(\alpha)x(t-\alpha)d\alpha\right)x(t-\tau)\right]=0\]
得到连续维纳-霍夫方程：
\[R_{sx}(\tau)=\int_0^\infty h(\alpha)R_x(\tau-\alpha)d\alpha,\qquad 0\le \tau<\infty\]
其中：
\[R_x(\tau)=E[x(t)x(t-\tau)]\]
\[R_{sx}(\tau)=E[s(t)x(t-\tau)]\]
结论：
- 连续维纳滤波也是由正交性原理得到。
- 核心方程是积分形式的维纳-霍夫方程。
- 直接求解该积分方程较困难，通常使用频谱因式分解法或预白化方法。
预白化方法
预白化方法的思想：若输入信号是白色的，维纳滤波器求解会非常简单；因此先把输入信号白化，再进行滤波。
若输入信号 \(x(t)\) 是白色的：
\[\Phi_x(\omega)=1\]
其自相关函数为：
\[R_x(\tau)=\delta(\tau)\]
代入连续维纳-霍夫方程：
\[R_{sx}(\tau)=\int_0^\infty h(\alpha)\delta(\tau-\alpha)d\alpha\]
因此：
\[h(t)=\begin{cases}R_{sx}(t), & t\ge 0\\0, & t<0\end{cases}\]
结论：当输入为白色信号时，最佳滤波器冲激响应就是互相关函数的因果部分。
一般情况下，\(x(t)\) 不是白色信号。
设输入功率谱可因式分解为：
\[\Phi_x(\omega)=G_x^+(\omega)G_x^-(\omega)\]
其中：
- \(G_x^+(\omega)\) 对应因果稳定部分。
- \(G_x^-(\omega)\) 对应反因果部分。
白化滤波器：
\[W(j\omega)=\frac{1}{G_x^+(\omega)}\]
白化后信号：
\[y(t)=W(j\omega)x(t)\]
其功率谱为：
\[\Phi_y(\omega)=\Phi_x(\omega)|W(j\omega)|^2=1\]
于是 \(y(t)\) 为白色信号。
需要求 \(s(t)\) 与白化信号 \(y(t)\) 的互功率谱。
由：
\[y(t)=w(t)*x(t)\]
可得：
\[\Phi_{sy}(\omega)=W^*(j\omega)\Phi_{sx}(\omega)\]
又因为：
\[W(j\omega)=\frac{1}{G_x^+(\omega)}\]
所以：
\[\Phi_{sy}(\omega)=\frac{\Phi_{sx}(\omega)}{G_x^-(\omega)}\]
将其分解为因果部分与非因果部分：
\[\Phi_{sy}(\omega)=[\Phi_{sy}(\omega)]_++[\Phi_{sy}(\omega)]_-\]
对白化信号的最佳滤波器为：
\[G(j\omega)=[\Phi_{sy}(\omega)]_+\]
因此最终维纳滤波器为：
\[H(j\omega)=W(j\omega)G(j\omega)\]
即：
\[H(j\omega)=\frac{1}{G_x^+(\omega)}\left[\frac{\Phi_{sx}(\omega)}{G_x^-(\omega)}\right]_+\]
结论：
- 连续维纳滤波的预白化方法与离散 \(z\) 域方法思想一致。
- 先将输入 \(x(t)\) 白化。
- 对白化后的输入求简单维纳滤波器。
- 最终滤波器为白化滤波器与白化域最佳滤波器串联。
预白化法步骤：
- 对输入功率谱做因式分解：
\[\Phi_x(\omega)=G_x^+(\omega)G_x^-(\omega)\]
- 构造白化滤波器：
\[W(j\omega)=\frac{1}{G_x^+(\omega)}\]
- 计算白化后输入与期望信号的互功率谱：
\[\Phi_{sy}(\omega)=\frac{\Phi_{sx}(\omega)}{G_x^-(\omega)}\]
- 取因果部分：
\[[\Phi_{sy}(\omega)]_+\]
- 得到最终滤波器：
\[H(j\omega)=\frac{1}{G_x^+(\omega)}\left[\frac{\Phi_{sx}(\omega)}{G_x^-(\omega)}\right]_+\]
特点：
- 适用于输入功率谱为有理函数、可谱因式分解的情况。
- 把复杂的维纳-霍夫积分方程转化为谱分解和因果部分提取问题。
- 本质是“先白化，再滤波”。
卡尔曼滤波
卡尔曼滤波解决的问题：对随机动态系统的状态变量进行递推最优估计。
课件中的基本例子是一阶 AR 模型：
\[x(n)=ax(n-1)+w(n)\]
观测方程：
\[y(n)=x(n)+v(n)\]
其中：
- \(x(n)\) 是系统状态。
- \(y(n)\) 是观测数据。
- \(w(n)\) 是输入白噪声。
- \(v(n)\) 是测量引入的白噪声。
卡尔曼滤波的基本思想：用前一个状态的估计值和最近一个观测数据来估计状态变量的当前值。
状态方程和量测方程：
\[x_k=A_kx_{k-1}+w_{k-1}\]
\[y_k=C_kx_k+v_k\]
其中：
- \(x_k\)：第 \(k\) 步状态变量。
- \(y_k\)：第 \(k\) 步观测数据。
- \(A_k\)：状态变量之间的增益矩阵，可以随时间变化。
- \(C_k\)：状态变量与输出信号之间的增益矩阵，可以随时间变化。
- \(w_{k-1}\)：输入白噪声。
- \(v_k\)：观测白噪声。
噪声假设：
\[E[w_k]=0\]
\[E[v_k]=0\]
\[E[w_kw_j]=Q_k\delta_{kj}\]
\[E[v_kv_j]=R_k\delta_{kj}\]
其中：
\[\delta_{kj}=\begin{cases}1, & k=j\\0, & k\ne j\end{cases}\]
基本递推思想：
先不考虑输入噪声 \(w_k\) 和观测噪声 \(v_k\) 的影响，得到状态变量的预测值。
状态预测：
\[\hat{x}'_k=A_k\hat{x}_{k-1}\]
量测预测：
\[\hat{y}'_k=C_k\hat{x}'_k\]
代入状态预测：
\[\hat{y}'_k=C_kA_k\hat{x}_{k-1}\]
输出信号估计误差：
\[\tilde{y}_k=y_k-\hat{y}'_k\]
即：
\[\tilde{y}_k=y_k-C_kA_k\hat{x}_{k-1}\]
课件中的校正思想：用输出信号的估计误差 \(\tilde{y}_k\) 来校正状态变量的预测值。
校正公式：
\[\hat{x}_k=\hat{x}'_k+H_k(y_k-\hat{y}'_k)\]
即：
\[\hat{x}_k=A_k\hat{x}_{k-1}+H_k(y_k-C_kA_k\hat{x}_{k-1})\]
其中 \(H_k\) 是增益矩阵，本质上是一个加权矩阵。
卡尔曼滤波的关键：求 \(H_k\) 的最佳值，使状态估计误差的均方值最小。
状态估计误差：
\[\varepsilon_k=x_k-\hat{x}_k\]
误差均方值：
\[E[\varepsilon_k^2]\]
优化目标：
\[\min_{H_k}E[\varepsilon_k^2]\]
对于矢量状态，目标是使误差协方差矩阵尽可能小。
一阶标量情形：
状态方程：
\[x_k=ax_{k-1}+w_{k-1}\]
量测方程：
\[y_k=x_k+v_k\]
状态预测：
\[\hat{x}'_k=a\hat{x}_{k-1}\]
量测预测：
\[\hat{y}'_k=\hat{x}'_k\]
输出误差：
\[\tilde{y}_k=y_k-\hat{x}'_k\]
校正公式：
\[\hat{x}_k=\hat{x}'_k+H_k(y_k-\hat{x}'_k)\]
即：
\[\hat{x}_k=a\hat{x}_{k-1}+H_k(y_k-a\hat{x}_{k-1})\]
也可写成：
\[\hat{x}_k=(1-H_k)a\hat{x}_{k-1}+H_ky_k\]
含义：
- 若 \(H_k\) 大，则更相信当前观测 \(y_k\)。
- 若 \(H_k\) 小，则更相信模型预测 \(a\hat{x}_{k-1}\)。
误差方差递推：
设上一时刻估计误差方差为：
\[P_{k-1}=E[(x_{k-1}-\hat{x}_{k-1})^2]\]
预测误差方差：
\[P'_k=a^2P_{k-1}+Q\]
其中 \(Q\) 是输入噪声 \(w_k\) 的方差。
增益：
\[H_k=\frac{P'_k}{P'_k+R}\]
其中 \(R\) 是观测噪声 \(v_k\) 的方差。
更新后的误差方差：
\[P_k=(1-H_k)P'_k\]
因此完整递推为：
预测：
\[\hat{x}'_k=a\hat{x}_{k-1}\]
\[P'_k=a^2P_{k-1}+Q\]
增益：
\[H_k=\frac{P'_k}{P'_k+R}\]
校正：
\[\hat{x}_k=\hat{x}'_k+H_k(y_k-\hat{x}'_k)\]
误差方差更新：
\[P_k=(1-H_k)P'_k\]
矢量形式：
状态预测：
\[\hat{x}'_k=A_k\hat{x}_{k-1}\]
预测误差协方差：
\[P'_k=A_kP_{k-1}A_k^T+Q_k\]
量测预测：
\[\hat{y}'_k=C_k\hat{x}'_k\]
输出误差：
\[\tilde{y}_k=y_k-\hat{y}'_k\]
增益矩阵：
\[H_k=P'_kC_k^T(C_kP'_kC_k^T+R_k)^{-1}\]
状态校正：
\[\hat{x}_k=\hat{x}'_k+H_k(y_k-C_k\hat{x}'_k)\]
误差协方差更新：
\[P_k=(I-H_kC_k)P'_k\]
卡尔曼滤波递推流程：
第一步：状态预测
\[\hat{x}'_k=A_k\hat{x}_{k-1}\]
第二步：误差协方差预测
\[P'_k=A_kP_{k-1}A_k^T+Q_k\]
第三步：计算增益矩阵
\[H_k=P'_kC_k^T(C_kP'_kC_k^T+R_k)^{-1}\]
第四步：利用观测误差修正状态估计
\[\hat{x}_k=\hat{x}'_k+H_k(y_k-C_k\hat{x}'_k)\]
第五步：更新误差协方差
\[P_k=(I-H_kC_k)P'_k\]
卡尔曼滤波特点：
-
递推算法，在时域内设计滤波器。
-
不需要保存全部过去观测值，只需要前一时刻的估计值和当前观测值。
-
用状态方程描述状态变量的动态变化规律。
-
可以处理平稳随机过程，也可以处理非平稳随机过程。
-
适用于多维随机过程的估计。
-
误差准则仍是均方误差最小。
-
关键是计算最优增益矩阵 \(H_k\)。
与维纳滤波的关系：
卡尔曼滤波和维纳滤波都是线性最优估计器，准则都是均方误差最小。
卡尔曼滤波有一个过渡过程，其结果在开始阶段与维纳滤波不完全相同。
当卡尔曼滤波达到稳态后，其结果与维纳滤波相同。
结论：
\[\text{卡尔曼滤波稳态结果}=\text{维纳滤波结果}\]
原因：二者本质上都是以最小均方误差为准则的线性估计器。