信号和系统的特性和建模 这部分内容在 随机信号处理(理论)中有更基础的推导
信号和系统的特性描述 对于信号而言，其特性有 平稳性、周期性、功率谱密度等，它们用于建立信号的数学模型或描述其统计行为，并且我们通过自相关函数、傅里叶变换等方法提取信号的关键特性。对于系统建模而言，我们通过系统的时域或频域响应、系统的传递函数、极点与零点的分布等来提取系统等关键特性。
下面我们从估计离散信号自相关函数开始着手刻画信号的 统计特性 ，自相关函数反映了信号在不同时间延迟下的相似性，描述信号随时间变化的内在规律，是随机过程建模的基础
对于一个离散信号 x(n) ，其自相关函数定义为：
R_{xx}(\tau) = \mathbb{E}[x(n) x^*(n-\tau)]
这个公式非常好，很漂亮，但是也有一个实际问题，我们采样得到的 x(n) 是有限样本序列，我们不可能对全体信号进行期望运算，因此需要通过有限样本（长度为 N ）来近似估计真实的自相关函数。
期望操作 \mathbb{E}[\cdot] 的近似可以用 样本均值 替代。注意，用样本均值替代期望的合理性基于大数定律：在样本量 N \to \infty 时，样本均值会趋近于真实期望。因此我们最终给出:
离散信号自相关函数的估计公式
\hat{R}_{xx}(\tau) = \frac{1}{N - \tau} \sum_{n=\tau}^{N-1} x(n) \cdot x^*(n-\tau)
E\left[ \hat{R}_{xx}(\tau) \right] = \frac{1}{N - \tau} E\left[ \sum_{n=\tau}^{N-1} x(n) \cdot x^*(n-\tau) \right] = R_{xx}(\tau)
估计的自相关函数 \hat{R} _{xx}(\tau) 的期望值等于真实的自相关函数 R _{xx}(\tau) 。这说明，估计的平均值不会系统性地偏离真实值，也就是说明这个估计量是无偏的。再通俗来讲，\hat{R} _{xx}(\tau) 并不是每时每刻都等于真实的 R _{xx}(\tau) ，但是\hat{R} _{xx}(\tau) 所有时刻平均起来等于 R _{xx}(\tau) ，即这个估计方法是没有偏差的，可以用的。严谨一点叫做，取期望消除了样本随机性对结果的影响，体现了估计量的长期统计行为 。
证明:
E\left[ \hat{R}{xx}(\tau) \right] = E\left[ \frac{1}{N-\tau} \sum_{n=\tau}^{N-1} x(n) x^*(n-\tau) \right]
E\left[ \hat{R}{xx}(\tau) \right] = \frac{1}{N-\tau} \sum_{n=\tau}^{N-1} E\left[ x(n) x^*(n-\tau) \right]
在假设信号是平稳过程下，可以将真实自相关函数拿出来
E\left[ \hat{R}{xx}(\tau) \right] = \frac{1}{N-\tau} \sum_{n=\tau}^{N-1} R_{xx}(\tau)
自相关函数 R_{xx}(\tau) 与时间 n 无关
E\left[ \hat{R}{xx}(\tau) \right] = R{xx}(\tau) \cdot \frac{1}{N-\tau} \sum_{n=\tau}^{N-1} 1
E\left[ \hat{R}{xx}(\tau) \right] = R{xx}(\tau)
证明了 \hat{R}_{xx}(\tau) 是一个无偏估计量。
根据上面我们已经证明的离散信号自相关函数的估计公式，我们有以下结论:
假设 x(n) 是一个平稳过程，那么对于任何固定的 \tau ，这个估计值的期望值就是 R_{xx}(\tau) ，即：
E\left[ \hat{R}_{xx}(\tau) \right] = R_{xx}(\tau)
自相关函数的边界情况 我们有原式:
\hat{R}_{xx}(\tau) = \frac{1}{N - \tau} \sum_{n=\tau}^{N-1} x(n) \cdot x^*(n-\tau)
我们考虑自相关估计的一个特殊情况：当 \tau = N-1 时，只有一个样本 x(N-1) \cdot x^{*}(0) ，即
\hat{R}_{xx}(N-1) = x(N-1) \cdot x^*(0)
即，当 \tau = N-1 时，延迟值达到了信号长度的极限，仅剩下第一个样本 x(0) 和最后一个样本 x(N-1) ，自相关估计只能用这一个样本对来计算，没有太多统计意义。
因此我们可以说，当滞后值接近于信号长度的时候，其得到的自相关值是不可靠的，因为样本数量不足。
为此我们引入带偏估计
在带偏估计中，归一化系数是固定的 \frac{1}{N} 而不是延迟 \tau 调整的 \frac{1}{N - \tau} ，因此期望值会引入一个 \frac{N - \tau}{N} 的系数，使得估计量不再是无偏的。
带偏估计的自相关公式为：
\hat{R}_{xx}(\tau) = \frac{1}{N} \sum_{n=0}^{N-\tau-1} x(n) x(n-\tau)
带偏估计的期望值为：
假设 x(n) 是平稳过程，则：
E[\hat{R}_{xx}(\tau)] = \frac{N-\tau}{N} R_{xx}(\tau)
因此我们可见，其带偏估计的偏差来源就是这个因子 \frac{N-\tau}{N} ，这一因子随着 \tau 增大而减小，表明滞后值较大时，带偏估计会系统性地低估自相关值，适用于大样本估计。并且在边界这种极端情况下，也能工作。
Hankel 矩阵定义
Hankel 矩阵是一种特殊的矩阵，其特点是 沿反对角线上的元素相同 。在信号处理中，Hankel 矩阵是根据时间序列数据构造的，常用于表示信号的结构信息。
给定时间序列 x(0), x(1), x(2), \ldots, x(N+L-1) ，Hankel 矩阵的定义为：
H = \begin{bmatrix}
x(0) & x(1) & x(2) & \cdots & x(N-1) \\\\
x(1) & x(2) & x(3) & \cdots & x(N) \\\\
x(2) & x(3) & x(4) & \cdots & x(N+1) \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
x(L) & x(L+1) & x(L+2) & \cdots & x(N+L-1)\\
\end{bmatrix}
行的数量为 L+1 ，列的数量为 N 。
Hankel 矩阵通过序列 x(k) 的排列捕获信号的时序相关特性，尤其适合分析平稳信号的统计特性，在自回归模型和功率谱估计中用于估计信号的自相关矩阵
建模与傅里叶级数展开 离散时间信号 x(n) 傅里叶级数展开公式
x(n) = \sum_{k} A_k e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} n}
一个离散时间信号 x(n) 可以通过复指数形式的傅里叶级数展开成不同频率成分的叠加。其中， A_k 是第 k 个频率分量的复幅度，包含信号在对应离散频率 f_k 上的幅度和相位信息。
假设信号 x(n) 是一个平稳白噪声 ，即每个时刻的信号值彼此独立且服从相同的概率分布（通常为均值为零的高斯分布），即:
零均值
E[x(n)] = 0
自相关函数
R_{xx}(\tau) = \sigma^2 \delta(\tau)
表明信号的平均值为0，且信号在 \tau \neq 0 时无相关性，各时间点的值是独立的。
频域特性 -- 功率谱密度
白噪声的功率谱密度 S_{xx}(f) 是自相关函数 R_{xx}(\tau) 的傅里叶变换：
S_{xx}(f) = TF \left( R_{xx}(\tau) \right) = \sum_{\tau} R_{xx}(\tau) e^{-j 2 \pi f \tau}
代入 R_{xx}(\tau) = \sigma^2 \delta(\tau) ：
S_{xx}(f) = \sum_{\tau} \sigma^2 \delta(\tau) e^{-j 2 \pi f \tau}
因为 \delta(\tau) 仅在 \tau = 0 时非零：
S_{xx}(f) = \sigma^2
白噪声的功率谱密度 S_{xx}(f) 在所有频率 f 上都是均匀分布的常数 \sigma^2 ，即所有频率上的功率密度相同。换句话说，在白噪声的建模中，傅里叶级数揭示了白噪声信号在频域的均匀性，即没有特定的频率成分占主导。在时域上，自相关函数为 \sigma^2 \delta(\tau) ，说明信号无记忆、各时刻独立。在频域上，功率谱密度为常数 \sigma^2 ，说明信号在所有频率上具有相同的能量。
移动平均 (MA) 过程 移动平均（MA）过程是一种线性时间序列模型，其输出 x(n) 是当前和过去有限个输入信号 u(n) 的加权和：
x(n) = \sum_{j=0}^{q} b_j \cdot u(n - j)
• u(n) 是输入信号，通常假设为白噪声，零均值 E[u(n)] = 0 和有限方差 E[u^2(n)] = \sigma^2 。
• b_j 是 MA 过程的系数，表示每个输入信号对输出的权重。
• q 是 MA 模型的阶数（滤波器的记忆长度），输出信号 x(n) 依赖于 u(n) 的过去 q+1 个样本。
输出信号 x(n) 仅由有限长度的输入序列 u(n), u(n-1), \ldots, u(n-q) 决定，因此称为短期记忆过程 （processus à mémoire courte）
将输入信号 u(n) 和输出信号 x(n) 进行 z -变换，得：
X(z) = \sum_{j=0}^q b_j \cdot U(z) \cdot z^{-j}
因此得到了 MA 过程的传递函数：
H(z) = \frac{X(z)}{U(z)} = \sum_{j=0}^{q} b_j \cdot z^{-j}
H(z) = b_0 + b_1 z^{-1} + b_2 z^{-2} + \cdots + b_q z^{-q}
从传递函数中，我们可以看出，其 H(z) 仅包含有限项，滤波器只依赖于 q+1 个输入值，是一种 有限脉冲响应滤波器 (FIR Filter)。同时H(z) 是 z^{-1} 的多项式，假设具有 q 个零点 z_1, z_2, \ldots, z_q ，则H(z) 可以表示为：
H(z) = b_0 \prod_{j=1}^{q} (1 - z_j z^{-1})
其中，b_0 是增益常数，z_j 是系统的零点。
MA 模型的性质
**短期记忆性: ** MA 过程仅依赖于有限长度的输入信号，因此是 短期记忆过程
平稳性 : 输入信号 u(n) 是平稳白噪声，且 MA 过程只进行加权和操作，因此 MA 模型生成的输出信号 x(n) 也是平稳过程自相关函数: MA 过程的自相关函数 R_{xx}(\tau) 仅在 |\tau| \leq q 时非零(只关联到过去 q+1 个输入)功率谱密度: MA 过程的功率谱密度是有限脉冲响应滤波器作用于白噪声谱密度的结果，由滤波器传递函数的模平方决定
S_{xx}(f) = \sigma^2 \cdot |H(e^{j 2 \pi f})|^2
零点的复数表示 零点 z_j 是复平面上的点，通过复数的极坐标形式，可以将 z_j 表示为：
z_j = |z| \cdot e^{j \varphi_j}
相位角 \varphi_j 的计算 \varphi_j = \frac{2 \pi f_j}{f_{\text{ech}}}
因为物理频率 f_j 的范围是 [0, f_{\text{ech}}/2] （奈奎斯特采样定理的限制），归一化角频率 \varphi_j 的范围为 [0, \pi] （对应复平面上的上半圆）。
举个例子：
假设采样频率 f_{\text{ech}} = 1000 \, \text{Hz} ，零点对应频率 f_j = 250 \, \text{Hz} ，则：
\varphi_j = \frac{2 \pi f_j}{f_{\text{ech}}} = \frac{2 \pi \cdot 250}{1000} = \frac{\pi}{2}
零点的相位角 \varphi_j = \pi/2 ，表示零点位于单位圆上的正虚轴位置，滤波器在频率 f_j = 250 \, \text{Hz} 附近会显著衰减信号。若零点在远离单位圆的位置，则对频率响应的影响较小。
自回归 (AR) 过程 AR 模型定义 x(n) = -\sum_{i=1}^{p} a_i , x(n - i) + u(n)
• p 是 AR 模型的阶数，即使用的过去值的数量。
• a_i 是 AR 系数，代表过去信号对当前信号的影响。
• u(n) 是输入噪声项，通常假设为白噪声。
当前的输出 x(n) 是基于过去 p 个输出值的加权和，再加上一个随机噪声 u(n) 。AR 模型用于描述时间序列数据的内在依赖性，广泛应用于信号处理和时间序列分析。
传递函数 H(z) = \frac{1}{1 + \sum_{i=1}^{p} a_i \cdot z^{-i}}
进行因式分解
H(z) = \frac{1}{\prod_{i=1}^{p} (1 - p_i z^{-1})}
通过该传递函数，我们可以分析系统的频率响应和极点。
AR 模型的极点分布在复平面上会影响频率响应。如果极点接近单位圆，系统在相应频率会出现较强的响应；如果极点远离单位圆，系统的响应则较弱。
ARMA (Auto-Regressive Moving Average) 模型 ARMA 模型定义 x(n) = -\sum_{i=1}^{p} a_i , x(n - i) + \sum_{j=0}^{q} b_j , u(n - j)
在 ARMA 模型中：
• AR 部分 ：自回归部分是由前 p 个输出值 x(n - i) 的加权和构成，其系数为 a_i 。
• MA 部分 ：移动平均部分是由当前和前 q 个输入噪声 u(n - j) 的加权和构成，其系数为 b_j 。
• 该模型结合了自回归和移动平均的特性，使其能够更灵活地描述时间序列数据中的依赖关系。
ARMA 模型的传递函数可以通过 Z 变换表示为：
H(z) = \frac{\sum_{j=0}^{q} b_j \cdot z^{-j}}{1 + \sum_{i=1}^{p} a_i \cdot z^{-i}}
• 分子 ：分子部分是 MA 部分的 Z 变换形式，由 b_j 系数和 z^{-j} 项构成。
• 分母 ：分母部分是 AR 部分的 Z 变换形式，由 a_i 系数和 z^{-i} 项构成。
进一步将传递函数因式分解，得到：
H(z) = \frac{\prod_{j=0}^{q} (1 - z_j z^{-1})}{\prod_{i=1}^{p} (1 - p_i z^{-1})}
• z_j 表示 MA 部分的零点。
• p_i 表示 AR 部分的极点。
当零点接近单位圆时，滤波器会在相应频率范围内有显著衰减；当极点接近单位圆时，滤波器会在对应频率范围内有增强效果。
ARMA 滤波器的频率响应特性 即频谱图示为: 待补充
二阶移动平均 (MA) 模型的自相关函数计算 简单的二阶移动平均信号模型 x(n) = b_0 u(n) + b_1 u(n-1)
• u(n) 是输入信号，通常假设为白噪声。
• b_0 和 b_1 是滤波器的系数。
这个模型表示当前的输出 x(n) 是当前输入 u(n) 和前一个输入 u(n-1) 的加权和。
均值计算 E[x(n)] = b_0 E[u(n)] + b_1 E[u(n-1)] = 0
自相关函数：
R_{xx}(n_1, n_2) = E[x(n_1) \cdot x^*(n_2)]
\begin{align*}
R_{xx}(n_1, n_2) &= E\left[ (b_0 u(n_1) + b_1 u(n_1 - 1))(b_0 u^*(n_2) + b_1 u^*(n_2 - 1)) \right] \\
&= b_0^2 E[u(n_1) u^*(n_2)] + b_0 b_1 E[u(n_1) u^*(n_2 - 1)] + b_1 b_0 E[u(n_1 - 1) u^*(n_2)] + b_1^2 E[u(n_1 - 1) u^*(n_2 - 1)]
\end{align*}
为了简化表达式，引入滞后变量 \tau = n_1 - n_2
\begin{align*}
R_{xx}(\tau) &= b_0^2 R_{uu}(\tau) + b_0 b_1 R_{uu}(\tau + 1) + b_1 b_0 R_{uu}(\tau - 1) + b_1^2 R_{uu}(\tau) \\
&= (b_0^2 + b_1^2) R_{uu}(\tau) + b_0 b_1 \left( R_{uu}(\tau + 1) + R_{uu}(\tau - 1) \right)
\end{align*}
其中 R_{uu}(\tau) 来表示输入白噪声 u(n) 的自相关函数
• 第一部分 (b_0^2 + b_1^2) R_{uu}(\tau) ，表示信号 x(n) 的自相关主项。
• 第二部分 b_0 b_1 \left( R_{uu}(\tau + 1) + R_{uu}(\tau - 1) \right) ，表示由于滞后一单位时间所带来的影响。
假设 u(n) 是白噪声过程，其自相关函数 R_{uu}(\tau) = \sigma^2 \delta(\tau)
R_{xx}(\tau) = (b_0^2 + b_1^2) \sigma^2 \delta(\tau) + b_0 b_1 \sigma^2 \left( \delta(\tau + 1) + \delta(\tau - 1) \right)
• 当 \tau = 0 时，R_{xx}(0) = (b_0^2 + b_1^2) \sigma^2 。
• 当 \tau = \pm 1 时，R_{xx}(\pm 1) = b_0 b_1 \sigma^2 。
假设 x 是一个包含 N 个样本的随机向量：
x = \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{bmatrix}
每个元素 x_i 是信号在不同时间的值。
x 的自相关矩阵 E(x \cdot x^T) ，即：
E(x \cdot x^T) = E \left( \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{N-1} \end{bmatrix} \begin{bmatrix} x_0 & x_1 & \cdots & x_{N-1} \end{bmatrix} \right)
E(x \cdot x^T) = \sigma^2 \begin{bmatrix}
b_0^2 + b_1^2 & b_0 b_1 & 0 & \cdots & 0 \\\\
b_0 b_1 & b_0^2 + b_1^2 & b_0 b_1 & \cdots & 0 \\\\
0 & b_0 b_1 & b_0^2 + b_1^2 & \cdots & 0 \\\\
\vdots & \vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & 0 & \cdots & b_0^2 + b_1^2\\
\end{bmatrix}
• 主对角线上，每个元素都是 (b_0^2 + b_1^2) \sigma^2 ，对应于 R_{xx}(0) 。
• 在上下相邻的次对角线上，每个元素是 b_0 b_1 \sigma^2 ，对应于 R_{xx}(\pm 1) 。
二阶移动平均模型 (MA(2)) 的自相关特性：
• 主对角线 表示信号本身的方差。
• 次对角线 反映了相邻时刻之间的相关性，由 b_0 和 b_1 决定。
这种矩阵结构称为带状矩阵 ，仅在主对角线及其附近具有非零元素，适合表示具有有限“记忆”的信号模型。
TUMA（时变移动平均模型） 是移动平均模型的时变版本，其中滤波器系数随时间变化。模型定义如下：
x(n) = \sum_{j=0}^{p} b_j(n) \cdot u(n-j)
• b_j(n) 是时变的系数
• p 是模型的阶数。
这里的时变性体现在 b_j(n) 随时间 n 而变化，这使得 TUMA 模型能够捕捉信号的时变特性。TUMA 模型特别适用于描述具有非平稳性的信号，因为系数的时间依赖性使其能够更好地适应信号的变化。
另一种形式：
x(n) = \sum_{j=0}^{p} B_j(u(n)) \quad \text{这里不懂，可能抄错}
TVAR（时变自回归模型） （Time-Variant Auto-Regressive）是自回归模型的时变版本
x(n) = -\sum_{i=1}^{p} a_i(n) \cdot x(n-i) + u(n)
• a_i(n) 是时变的 AR 系数，它随时间 n 变化。
• p 是模型的阶数。
定义了系数 a_i(n) 的更新规则：
1. 第一种形式
a_i(n) = a_i(n-1) + N(n)
其中 N(n) 是一个BBGC 独立于 u(n)
2. 系数 a_i(n) 还可以写成如下形式：
a_i(n) = \sum_{j} \beta_{ij} b_j(n)
这意味着 a_i(n) 可以通过一组基函数 b_j(n) 的线性组合来表示，其中 \beta_{ij} 是权重。这种表示允许我们用基函数表示时变系数，从而降低模型的复杂度。
最终多维状态向量形式
\underline{x}(n) = -\sum_{i=1}^{p} A_i , x(n-i) + u(n)
• A_i 不再是单一标量，而是一个矩阵系数
长记忆过程（Processus à mémoire longue） 的模型使用分数阶积分来模拟长记忆的白噪声
滤波器的传递函数 H(z) 定义为：
H(z) = \frac{X(z)}{U(z)} = (1 - z^{-1})^{-d}
功率谱密度计算
\text{DSP}_{x(f)} = |H(z)|^2 \cdot \sigma_u^2
• \sigma_u^2 是输入白噪声 u(n) 的方差。
• 将 z 替换为 e^{j 2 \pi f / f_{\text{ech}}} ，得到频率域的响应 H(e^{j 2 \pi f / f_{\text{ech}}}) 。
代入 H(z) = (1 - e^{-j 2 \pi f / f_{\text{ech}}})^{-d} ：
\text{DSP}_{x(f)} = \left| \frac{1}{(1 - e^{-j 2 \pi f / f*{\text{ech}}})^d} \right|^2 \cdot \sigma_u^2
\text{DSP}_{x(f)} = \frac{\sigma_u^2}{\left| e^{-j \pi f / f*{\text{ech}}} \cdot 2j \sin\left( \frac{\pi f}{f_{\text{ech}}} \right) \right|^{2d}}
化简得到：
\text{DSP}_{x(f)} = \frac{\sigma_u^2}{2^{2d} \sin^{2d}\left( \frac{\pi f}{f*{\text{ech}}} \right)}
这个公式表明，输出信号 x(n) 的功率谱密度 \text{DSP}_{x(f)} 在频率 f \to 0 附近会趋于无穷大，这说明该系统具有长记忆性，低频成分较强。
ARFIMA (AutoRegressive Fractionally Integrated Moving Average) 模型 ARFIMA 模型的传递函数可以写为：
H(z) = \frac{\prod_{j=1}^{q} (1 - z_j z^{-1})}{(1 - z^{-1})^d \prod_{i=1}^{p} (1 - p_i z^{-1})}
• 分子部分 ：\prod_{j=1}^{q} (1 - z_j z^{-1}) 表示 MA（移动平均）部分，它包含了零点 z_j 。
• 分母部分 ：
• (1 - z^{-1})^d 表示分数阶积分项，通常 d 为分数，用于引入长记忆效应。
• \prod_{i=1}^{p} (1 - p_i z^{-1}) 表示 AR（自回归）部分，包含了极点 p_i 。
ARFIMA 模型是一种结合了自回归、分数阶积分和移动平均的模型，适用于描述具有长记忆性质的时间序列。通过分数阶积分项 (1 - z^{-1})^d ，模型能够捕捉时间序列的长记忆特性。
非线性模型（Non-linear Model）
一个与汽车模型相关的复杂动态过程
这个模型可以用以下表达式表示：
x(n) = \sum_{j=0}^{q} b_j u(n - j) + h(0, 0) u^2(n) + h(0, 1) u(n) u(n - 1) + h(0, 2) u(n) u(n - 2) + \dots
• 线性部分 ：\sum_{j=0}^{q} b_j u(n - j) ，这是传统的移动平均 (MA) 部分，表示系统对输入信号 u(n) 的线性响应。
• 非线性部分 ：包含 u(n) 的高阶项，例如 u^2(n) 、u(n) u(n - 1) 等。
• h(0, 0) u^2(n) 是输入 u(n) 的平方项，代表了二次非线性。
• h(0, 1) u(n) u(n - 1) 是输入和其滞后项的乘积，代表了跨时间的非线性耦合。
• 高阶项继续扩展为更高次方和不同滞后项的组合。
非线性项表示系统的响应不仅依赖于当前和过去的输入值，还依赖于这些输入值的非线性组合。
更高阶的非线性项例如：
h(1, 0, 0) u^3(n) + h(0, 0, 1) u^2(n) u(n - 1) + \dots
对于一个信号 x(n) = u^2(n)
x(n) = \cos^2\left(2 \pi n \frac{f_0}{f_{\text{ech}}}\right)
这个表达式是一个余弦平方信号，频率为 f_0 （归一化到采样频率 f_{\text{ech}} ）。我们可以利用三角恒等式（使用欧拉公式 ）将其展开。
已知 \cos^2(\theta) = \frac{1 + \cos(2\theta)}{2} ，同时将 \cos^2 转换为指数形式
可以将 x(n) 表达为：
x(n) = \frac{e^{j 4 \pi n \frac{f_0}{f_{\text{ech}}}} + 2 + e^{-j 4 \pi n \frac{f_0}{f_{\text{ech}}}}}{4}
矩阵形式的表示：
x(n) = \begin{bmatrix}
u(n) & u(n-1) & u(n-2) & \dots & u(n-q) & u^2(n) & u(n) u(n-1) & u(n) u(n-2) & \dots
\end{bmatrix}
\begin{bmatrix}
b_0 \\
b_1 \\
b_2 \\
\vdots \\
b_q \\
h(0,0) \\
h(0,1) \\
h(0,2) \\
\vdots
\end{bmatrix}
估计部分
两个不同类型：
1. Par bloc（分块估计） 将数据分成固定长度的块，然后在每个块上独立地计算所需的估计值
2. Échantillon par échantillon（逐样本估计） 对每个采样点独立地进行估计，不进行分块
估计模型 的输出信号 y(n) 的形式：
y(n) = \sum_{k} A_k e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} n} + b(n)
• \sum_{k} A_k e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} n} 表示一系列复指数信号的叠加，每个信号的频率为 f_k ，采样频率为 f_{\text{ech}} 。
• 这里 A_k 被视为一个随机变量 并且 E[A_k] = 0
• b(n) 是一个噪声项，假设为白噪声。
E[y(n)] = 0
自相关函数 R_{yy}(n_1, n_2) = E[y(n_1) y^*(n_2)]
\begin{align*}
R_{yy}(n_1, n_2) &= E \left[ \left( \sum_{k} A_k e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} n_1} + b(n_1) \right) \left( \sum_{l} A_l^* e^{-j 2 \pi \frac{f_l}{f_{\text{ech}}} n_2} + b^*(n_2) \right) \right] \\
&= E \left[ \sum_{k} \sum_{l} A_k A_l^* e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} n_1} e^{-j 2 \pi \frac{f_l}{f_{\text{ech}}} n_2} \right] + E[b(n_1) b^*(n_2)]
\end{align*}
噪声 b(n) 是白噪声，且 A_k 的均值为零并且彼此独立，这可以进一步简化为：
R_{yy}(n_1, n_2) = \sum_{k} E[|A_k|^2] e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} (n_1 - n_2)} + \sigma_b^2 \delta(n_1 - n_2)
• 第一项 \sum_{k} E[|A_k|^2] e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} \tau} 表示信号部分的自相关，包含各个频率分量的贡献，其中 \tau = n_1 - n_2 。
• 第二项 \sigma_b^2 \delta(\tau) 表示噪声部分的自相关，仅在 \tau = 0 时为非零。
因此，最终的自相关函数表示为：
R_{yy}(\tau) = \sum_{k} E[|A_k|^2] e^{j 2 \pi \frac{f_k}{f_{\text{ech}}} \tau} + \sigma_b^2 \delta(\tau)
定义向量 \underline{Y} 为：
\underline{Y} = \begin{bmatrix} y(n) \\\\ y(n+1) \\\\ \vdots \\\\ y(n+N-1)\\ \end{bmatrix}
自相关矩阵 R_y 定义为：
R_y = E[\underline{Y} \cdot \underline{Y}^H]
• \underline{Y}^H 是 \underline{Y} 的共轭转置（Hermitian 转置），在计算复数信号的自相关时特别重要。
• R_y 是一个 N \times N 的矩阵。
R_y = \begin{bmatrix}
r _{yy}(0) & r _{yy}(-1) & \dots & r _{yy}(-N+1) \\\\
r _{yy}(1) & r _{yy}(0) & \dots & r _{yy}(-N+2) \\\\
\vdots & \vdots & \ddots & \vdots \\\\
r _{yy}(N-1) & r _{yy}(N-2) & \dots & r _{yy}(0)\\
\end{bmatrix}
• 主对角线 上的元素为 r_{yy}(0) ，表示零滞后的自相关。
• 次对角线 上的元素为 r_{yy}(\pm 1) ，表示滞后为 \pm 1 的自相关。
• 该矩阵具有 Toeplitz 结构 ，即每一条对角线上的元素相同。
\underline{S} \underline{S}^H = \begin{bmatrix}
1 & e^{-j 2 \pi \frac{f_0}{f_{\text{ech}}}} & \cdots & e^{-j 2 \pi (N-1) \frac{f_0}{f_{\text{ech}}}} \\\\
e^{j 2 \pi \frac{f_0}{f_{\text{ech}}}} & 1 & \cdots & e^{-j 2 \pi (N-2) \frac{f_0}{f_{\text{ech}}}} \\\\
\vdots & \vdots & \ddots & \vdots \\\\
e^{j 2 \pi (N-1) \frac{f_0}{f_{\text{ech}}}} & e^{j 2 \pi (N-2) \frac{f_0}{f_{\text{ech}}}} & \cdots & 1\\
\end{bmatrix}
• 主对角线 ：所有元素都为 1，对应 e^{j 2 \pi \cdot 0} = 1 。
• 次对角线 ：依次为 e^{j 2 \pi \frac{f_0}{f_{\text{ech}}}} 、e^{j 2 \pi \cdot 2 \frac{f_0}{f_{\text{ech}}}} ，表示随着滞后增加的指数因子。
• 共轭对称性 ：由于 \underline{S} \underline{S}^H 是由 \underline{S} 和其共轭转置 \underline{S}^H 相乘得到，因此矩阵是共轭对称的，即 (\underline{S} \underline{S}^H) _{i,j} = \overline{(\underline{S} \underline{S}^H) _{j,i}} 。
加上噪声项后，自相关矩阵可以表示为：
R_y = E[|A|^2] \cdot \underline{S} \underline{S}^H + \sigma_b^2 I_N
• 信号部分 E[|A|^2] \cdot \underline{S} \underline{S}^H ：表示信号的自相关结构。
• 噪声部分 \sigma_b^2 I_N ：均匀分布在对角线上，反映了白噪声的自相关性。
矩阵 \underline{S} 表示为一个包含指数因子的向量：
\underline{S} = \begin{bmatrix}
1 \\\\
e^{j 2 \pi \frac{f_0}{f_{\text{ech}}}} \\\\
\vdots \\\\
e^{j 2 \pi (N-1) \frac{f_0}{f_{\text{ech}}}}\\
\end{bmatrix}
R_y = E[|A_1|^2] \cdot \underline{S} \underline{S}^H + \sigma_b^2 I_N
R_y \underline{S} = E[|A|^2] \cdot \underline{S} (\underline{S}^H \underline{S}) + \sigma_b^2 \underline{S}
设 \underline{S}^H \underline{S} 值为 N （是一个标量）（范数平方）
= (N E[|A|^2] + \sigma_b^2) \underline{S}
说明 R_y 的一个特征值为 N E[|A|^2] + \sigma_b^2
假设 \underline{\psi}_i \perp \underline{S}
也就是说，\underline{\psi}_i 与 \underline{S} 正交。
\underline{S}^T \underline{\psi}_i^* = 0 \quad \text{和} \quad \underline{S}^H \underline{\psi}_i = 0
\underline{\psi}_i 属于 R_y 的特征空间中，并且与 \underline{S} 的方向正交
原式
R_y \underline{\psi}_i = E[|A|^2] \cdot \underline{S} \underline{S}^H \underline{\psi}_i + \sigma_b^2 \underline{\psi}_i = \sigma_b^2 \underline{\psi}_i
自相关矩阵在 \underline{\psi}_i 方向上仅由噪声贡献，因为与信号成分 \underline{S} 正交，因此信号成分的部分会为零，留下的仅为噪声方差的贡献：
R_y \underline{\psi}_i = \sigma_b^2 \underline{\psi}_i
表明，可以通过正交特征向量来分离信号和噪声成分
3 \times 3 特征向量的线性组合：
\psi_3(0) \cdot 1 + \psi_3(1) \cdot e^{j \frac{2\pi}{f_{\text{ech}}}} + \psi_3(2) \cdot e^{-j \frac{2\pi}{f_{\text{ech}}}} = 0
• \psi_3(0) 、\psi_3(1) 、\psi_3(2) 是特征向量的各个分量。
• e^{j \frac{2\pi}{f_{\text{ech}}}} 和 e^{-j \frac{2\pi}{f_{\text{ech}}}} 表示复数指数因子，用于描述矩阵的频率特性。
z = e^{j \frac{2\pi}{f_{\text{ech}}}}
将 z 代入后，表达式变为：
\psi_3(0) + \psi_3(1) \cdot z^{-1} + \psi_3(2) \cdot z^{-2} = 0
维纳滤波器
接下来我们讨论维纳滤波器。维纳滤波器是设计最小均方意义上的滤波器的一个方法，目标是得到我们的期望响应。我们看到维纳滤波器实际上是一种分块的矩阵方法，这引出了相关矩阵的逆矩阵乘以输入和期望响应的互相关向量的计算。
然而，当矩阵规模较大时，直接求逆可能并不现实。因此，另一种方法是采用在线计算，而不是分块计算。我们将引入维纳滤波器的一个扩展，即一种自适应方法 ，它会随着时间的推移动态调整冲激响应，以考虑误差。随着时间推移，这种方法会逐渐收敛到维纳滤波器。换句话说，自适应滤波器，尤其是 LMS（最小均方）、NLMS（归一化最小均方）等。再往后，我们会介绍维纳滤波器的变种，即通道滤波方法。
状态向量定义为：
\underline{x}(t) =
\begin{bmatrix}
x(t) \\\\
\dot{x}(t) \\\\
\ddot{x}(t)\\
\end{bmatrix}
状态方程可以写成：
\dot{\underline{x}}(t) =
\begin{bmatrix}
\dot{x}(t) \\\\
\ddot{x}(t) \\\\
\overset{\cdots}{x}(t)\\
\end{bmatrix}
\dot{\underline{x}}(t) =
\begin{bmatrix}
0 & 1 & 0 \\\\
0 & 0 & 1 \\\\
0 & 0 & 0\\
\end{bmatrix}
\begin{bmatrix}
x(t) \\\\
\dot{x}(t) \\\\
\ddot{x}(t)\\
\end{bmatrix}
+
\begin{bmatrix}
0 \\\\
0 \\\\
1\\
\end{bmatrix}
w(t)
\dot{\underline{x}}(t) =
A \underline{x}(t) + B w(t)
状态向量 x(t) 的解为：
x(t) = \exp(A(t-t_0)) \cdot x(t_0) + \int_{t_0}^t \exp(A(t-\tau)) B w(\tau) d\tau
给出采样条件：
\begin{cases}
t = (k+1)T_{\text{ech}} \\\\
t_0 = kT_{\text{ech}} \\\\
t - t_0 = T_{\text{ech}}\\
\end{cases}
原式:
x(t) = \exp(A \cdot T_{\text{ech}}) \cdot x(kT_{\text{ech}}) + \int_{t_0}^t \exp(A(t-\tau)) B w(\tau) d\tau
矩阵指数 \exp(A \cdot T_{\text{ech}}) 的泰勒展开形式为：
\exp(A \cdot T_{\text{ech}}) = I + A \cdot T_{\text{ech}} + \frac{A^2 \cdot T_{\text{ech}}^2}{2!} + \cdots
矩阵 A 的方幂具有如下性质：
A =
\begin{bmatrix}
0 & 1 & 0 \\\\
0 & 0 & 1 \\\\
0 & 0 & 0\\
\end{bmatrix}, \quad
A^2 =
\begin{bmatrix}
0 & 0 & 1 \\\\
0 & 0 & 0 \\\\
0 & 0 & 0\\
\end{bmatrix}, \quad
A^3 = 0
将以上性质代入矩阵指数的展开：
\exp(A \cdot T_{\text{ech}}) =
\begin{bmatrix}
1 & 0 & 0 \\\\
0 & 1 & 0 \\\\
0 & 0 & 1\\
\end{bmatrix}
+
\begin{bmatrix}
0 & T_{\text{ech}} & 0 \\\\
0 & 0 & T_{\text{ech}} \\\\
0 & 0 & 0\\
\end{bmatrix}
+
\begin{bmatrix}
0 & 0 & \frac{T_{\text{ech}}^2}{2} \\\\
0 & 0 & 0 \\\\
0 & 0 & 0\\
\end{bmatrix}
\exp(A \cdot T_{\text{ech}}) =
\begin{bmatrix}
1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} \\\\
0 & 1 & T_{\text{ech}} \\\\
0 & 0 & 1\\
\end{bmatrix}
离散系统状态更新表达式为：
x(k+1) =
\begin{bmatrix}
1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} \\
0 & 1 & T_{\text{ech}} \\
0 & 0 & 1
\end{bmatrix}
x(k)
+
\int_{kT_{\text{ech}}}^{(k+1)T_{\text{ech}}}
\begin{bmatrix}
1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} \\
0 & 1 & T_{\text{ech}} \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
0 \\
0 \\
1
\end{bmatrix}
w(\tau) \, d\tau
将积分中的变量替换如下：
\tau_1 = \tau - kT_{\text{ech}}\\\\
\tau = \tau_1 + kT_{\text{ech}}\\\\
d\tau_1 = d\tau\\
x(k+1) = \Phi x(k) + \int_0^{T_{\text{ech}}}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau_1)^2}{2} \\
T_{\text{ech}} - \tau_1 \\
1
\end{bmatrix}
w(\tau_1 + k T_{\text{ech}}) d\tau_1
后半部分是 u(k)，即:
u(k) = \int_0^{T_{\text{ech}}}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau_1)^2}{2} \\\\
T_{\text{ech}} - \tau_1 \\\\
1\\
\end{bmatrix}
w(\tau_1 + kT_{\text{ech}}) d\tau_1
这部分代表了模型噪音，其均值为0
系统过程噪声协方差矩阵 Q
Q = \mathbb{E}[u(k) \cdot u(k)^T]
假设条件： w(t) 是不随时间相关的零均值白噪声。
u(k) 的计算基于积分结果，表达为：
u(k) = w(kT_{\text{ech}})
\begin{bmatrix}
\frac{T_{\text{ech}}^2}{6} \\\\
\frac{T_{\text{ech}}}{2} \\\\
T_{\text{ech}}\\
\end{bmatrix}
将 u(k) 带入协方差矩阵定义，得到：
Q =
\begin{bmatrix}
\frac{T_{\text{ech}}^6}{36} & \frac{T_{\text{ech}}^5}{12} & \frac{T_{\text{ech}}^4}{6} \\\\
\frac{T_{\text{ech}}^5}{12} & \frac{T_{\text{ech}}^4}{4} & \frac{T_{\text{ech}}^3}{3} \\\\
\frac{T_{\text{ech}}^4}{6} & \frac{T_{\text{ech}}^3}{3} & T_{\text{ech}}^2\\
\end{bmatrix}
另一种情况
Q = \mathbb{E}[u(k) \cdot u(k)^T]
\mathbb{E}[u(k) \cdot u(k)^T] = \mathbb{E} \left[ \int_0^{T_{\text{ech}}} \int_0^{T_{\text{ech}}}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau_1)^2}{2} \\\\
T_{\text{ech}} - \tau_1 \\\\
1\\
\end{bmatrix}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau_1)^2}{2} & T_{\text{ech}} - \tau_1 & 1
\end{bmatrix}
w(\tau_1 + kT_{\text{ech}}) w(\tau_2 + kT_{\text{ech}}) d\tau_1 d\tau_2 \right]
假设 w(t) 是零均值、方差为 \sigma_w^2 的白噪声，根据白噪声的正交性：
\mathbb{E}[w(\tau_1 + kT_{\text{ech}}) w(\tau_2 + kT_{\text{ech}})] = \sigma_w^2 \delta(\tau_1 - \tau_2)
利用这个假设，双重积分简化为单积分：
Q = \sigma_w^2 \int_0^{T_{\text{ech}}}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau)^2}{2} \\\\
T_{\text{ech}} - \tau \\\\
1\\
\end{bmatrix}
\begin{bmatrix}
\frac{(T_{\text{ech}} - \tau)^2}{2} & T_{\text{ech}} - \tau & 1
\end{bmatrix}
d\tau
Q = q \cdot
\begin{bmatrix}
\frac{T_{\text{ech}}^5}{20} & \frac{T_{\text{ech}}^4}{8} & \frac{T_{\text{ech}}^3}{6} \\\\
\frac{T_{\text{ech}}^4}{8} & \frac{T_{\text{ech}}^3}{3} & \frac{T_{\text{ech}}^2}{2} \\\\
\frac{T_{\text{ech}}^3}{6} & \frac{T_{\text{ech}}^2}{2} & T_{\text{ech}}\\
\end{bmatrix}
其中， q 表示噪声强度 \sigma_w^2 。
状态向量 x(k+1) 表达为：
x(k+1) =
\begin{bmatrix}
x(k+1) \\
\dot{x}(k+1) \\
\ddot{x}(k+1)
\end{bmatrix}
=
\begin{bmatrix}
1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} \\
0 & 1 & T_{\text{ech}} \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
x(k) \\
\dot{x}(k) \\
\ddot{x}(k)
\end{bmatrix}
+ u(k)
第三个状态 \ddot{x}(k+1) 的更新公式为：
\ddot{x}(k+1) = \ddot{x}(k) + \mu_3(k)
x(k+1) =
\begin{bmatrix}
x(k+1) \\\\
\dot{x}(k+1) \\\\
\ddot{x}(k+1) \\\\
y(k+1) \\\\
\dot{y}(k+1) \\\\
\ddot{y}(k+1)\\
\end{bmatrix}
\begin{bmatrix}
1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} & 0 & 0 & 0 \\\\
0 & 1 & T_{\text{ech}} & 0 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 1 & T_{\text{ech}} & \frac{T_{\text{ech}}^2}{2} \\\\
0 & 0 & 0 & 0 & 1 & T_{\text{ech}} \\\\
0 & 0 & 0 & 0 & 0 & 1\\
\end{bmatrix}
\begin{bmatrix}
x(k) \\\\
\dot{x}(k) \\\\
\ddot{x}(k) \\\\
y(k) \\\\
\dot{y}(k) \\\\
\ddot{y}(k)\\
\end{bmatrix}
+
\begin{bmatrix}
u_x(k) \\\\
\vdots \\\\
u_y(k) \\\\
\vdots \\
\end{bmatrix}
协方差矩阵 Q 被表示为分块矩阵形式：
Q =
\begin{bmatrix}
Q_x & 0 \\\\
0 & Q_y\\
\end{bmatrix}
其中:
Q_x 是关于 x 方向状态变量的协方差矩阵 ===> Q_x = \mathbb{E}[u_x(k) u_x(k)^T]
Q_y 是关于 y 方向状态变量的协方差矩阵 ===> Q_y = \mathbb{E}[u_y(k) u_y(k)^T]
矩阵的非对角线块为零，表示 x 和 y 两个方向的噪声是独立的。
测量向量 y(k) 的公式为：
y(k) =
\begin{bmatrix}
x_b(k) \\
y_b(k)
\end{bmatrix}
=
\begin{bmatrix}
1 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0
\end{bmatrix}
\begin{bmatrix}
x(k) \\
\dot{x}(k) \\
\ddot{x}(k) \\
y(k) \\
\dot{y}(k) \\
\ddot{y}(k)
\end{bmatrix}
+
\begin{bmatrix}
b_x(k) \\
b_y(k)
\end{bmatrix}
其中，b(k) 是测量误差项 ，其协方差矩阵 R 被定义为测量噪声的期望值：
R = \mathbb{E}[b(k) \cdot b(k)^T]=
\begin{bmatrix}
\sigma_{b_x}^2 & 0 \\\\
0 & \sigma_{b_y}^2\\
\end{bmatrix}
即协方差矩阵 R 是测量噪声 b(k) 的二阶统计特性，描述了测量噪声的分布和强度。
因此，我们得到综合线性状态方程表达式
x(k+1) = \Phi x(k) + u(k)
y(k) = H x(k) + b(k)
举一个非线性的例子，雷达定位:
目标点与估计位置之间的欧几里得距离 C\Delta T 被定义为：
C\Delta T = \sqrt{(x_A(k) - x(k))^2 + (y_A(k) - y(k))^2 + (z_A(k) - z(k))^2}+ \xi(k)
即:
\Delta = g(x(k)) + \xi(k)
这就是个非线性方程
输入信号模型：
s(k) = -\sum_{i=1}^{p} a_i s(k-i) + u(k)
测量信号 y(k) 表达式
y(k) = s(k) + b(k)
b(k) \sim \mathcal{N}(0, \sigma_b^2)
为了使用状态空间表达式，将信号 s(k) 表示为一个状态向量：
\underline{x}_s(k) =
\begin{bmatrix}
s(k) \\\\
s(k-1) \\\\
s(k-2) \\\\
\vdots \\\\
s(k-p+1)\\
\end{bmatrix}
\underline{x}_s(k+1) =
\begin{bmatrix}
-a_1 & -a_2 & \cdots & -a_p \\\\
1 & 0 & \cdots & 0 \\\\
0 & 1 & \cdots & 0 \\\\
\vdots & \vdots & \ddots & \vdots \\\\
0 & 0 & \cdots & 1\\
\end{bmatrix}
\begin{bmatrix}
s(k) \\\\
s(k-1) \\\\
s(k-2) \\\\
\vdots \\\\
s(k-p+1)\\
\end{bmatrix}
+
\begin{bmatrix}
1 \\\\
0 \\\\
0 \\\\
\vdots \\\\
0\\
\end{bmatrix}
u(k)
\underline{x}_s(k+1) = \Phi \underline{x}_s(k-1) + \underline{u}(k)
测量信号 y(k) 的表达式为：
y(k) =
\begin{bmatrix}
1 & 0 & \cdots & 0
\end{bmatrix}
\begin{bmatrix}
s(k) \\\\
s(k-1) \\\\
\vdots \\\\
s(k-p+1)\\
\end{bmatrix}
+ b(k)
\underline{x}(k) =
\begin{bmatrix}
a_1(k) \\
a_2(k) \\
\vdots \\
a_p(k)\
\end{bmatrix}
=
I_d \cdot
\begin{bmatrix}
a_1(k-1) \\
a_2(k-1) \\
\vdots \\
a_p(k-1)\
\end{bmatrix}
+
\begin{bmatrix}
w_1(k) \\
w_2(k) \\
\vdots \\
w_p(k)\
\end{bmatrix}
自回归信号的测量方程为：
y(k) = s(k) =
\begin{bmatrix}
-s(k-1) & -s(k-2) & \cdots & -s(k-p)
\end{bmatrix}
\cdot
\begin{bmatrix}
a_1(k) \\\\
a_2(k) \\\\
\vdots \\\\
a_p(k)\\
\end{bmatrix}
+ u(k)
维纳滤波器
误差信号：
e(k) = d(k) - \underline{H}_N^T \underline{X}_N(k)
滤波输出：
\hat{d}(k) = \underline{H}_N^T \underline{X}_N(k)
滤波器系数向量：
\underline{H}_N^T = [h(0), h(1), \ldots, h(N-1)]
\hat{d}(k) = h * x(k) = \sum_{n=0}^{N-1} h(n) \cdot x(k-n)
\hat{d}(k) =
\begin{bmatrix}
h(0) & h(1) & \ldots & h(N-1)
\end{bmatrix}
\begin{bmatrix}
x(k) \\
x(k-1) \\
\vdots \\
x(k-N+1)
\end{bmatrix}
\hat{d}(k) = \underline{H}_N^T \underline{X}_N(k)
\hat{d}(k) = \underline{X}_N^T(k) \underline{H}_N
目标函数：
J = \mathbb{E}[e^2(k)]
误差项展开：
J = \mathbb{E}[(d(k) - \underline{H}_N^T \underline{X}_N(k))^2]
J = \mathbb{E}[d^2(k) - 2d(k) \cdot \underline{H}_N^T \underline{X}_N(k) + \underline{H}_N^T \underline{X}_N(k) \cdot \underline{X}_N^T(k) \cdot \underline{H}_N]
J = \mathbb{E}[d^2(k)] - 2 \underline{H}_N^T \mathbb{E}[d(k) \cdot \underline{X}_N(k)] + \underline{H}_N^T \mathbb{E}[\underline{X}_N(k) \cdot \underline{X}_N^T(k)] \cdot \underline{H}_N
J = \mathbb{E}[d^2(k)] - 2 \underline{H}N^T R{dx} + \underline{H}_N^T R_x \cdot \underline{H}_N
关于 \underline{H}_N 的优化条件 ：
-2 \underline{R}_{dx} + 2 R_x \underline{H}_N = 0
解出 \underline{H}_N ：
\underline{H}_N = R_x^{-1} \underline{R}{dx}
这就是维纳滤波器基本形式
自回归模型 ：
x(k) = -\sum_{i=1}^{p} a_i x(k-i) + u(k)
x(k) = \hat{x}(k) + u(k)
系数向量 \underline{H}_N ：
\underline{H}_N =
\begin{bmatrix}
a_1 \\\\
a_2 \\\\
\vdots \\\\
a_p\\
\end{bmatrix}
状态向量 \underline{X}_N(k) ：
\underline{X}_N(k) =
\begin{bmatrix}
-x(k-1) \\
-x(k-2) \\
\vdots \\
-x(k-p)
\end{bmatrix}
自相关矩阵 R_x ：
R_x =
\begin{bmatrix}
R_{xx}(0) & R_{xx}(1) & \cdots & R_{xx}(p-1) \\
R_{xx}(1) & R_{xx}(0) & \cdots & R_{xx}(p-2) \\
\vdots & \vdots & \ddots & \vdots \\
R_{xx}(p-1) & R_{xx}(p-2) & \cdots & R_{xx}(0)
\end{bmatrix}
输入信号 ：
d(k) = x(k)
互相关向量 R_{dx} ：
\underline{R}_{dx} = -
\begin{bmatrix}
R_{xx}(1) \\
\vdots \\
R_{xx}(p)
\end{bmatrix}
系数向量 \underline{H}_N
\underline{H}_N = - R_x^{-1} R{dx}
\underline{H}_N = - \begin{bmatrix}
R_{xx}(0) & R_{xx}(1) & \cdots & R_{xx}(p-1) \\
R_{xx}(1) & R_{xx}(0) & \cdots & R_{xx}(p-2) \\
\vdots & \vdots & \ddots & \vdots \\
R_{xx}(p-1) & R_{xx}(p-2) & \cdots & R_{xx}(0)
\end{bmatrix}^{-1} \quad
\begin{bmatrix}
R_{xx}(1) \\
\vdots \\
R_{xx}(p)
\end{bmatrix}
这就是 Yule-Walker 方程
自相关函数 R_{xx}(\tau) 的定义 ：
对于 \tau > 0 ：
R_{xx}(\tau) = \mathbb{E}[x(k) x(k-\tau)]
R_{xx}(\tau) = \mathbb{E}\left[\left(-\sum_{i=1}^p a_i x(k-i) + u(k)\right) x(k-\tau)\right]
R_{xx}(\tau) = -\sum_{i=1}^p a_i \mathbb{E}[x(k-i) x(k-\tau)] + \mathbb{E}[u(k) x(k-\tau)]
忽略噪声项（独立性假设） ：
R_{xx}(\tau) = -\sum_{i=1}^p a_i R_{xx}(\tau-i)
对于 \tau = 0 ：
R_{xx}(0) = \sum_{i=1}^p a_i R_{xx}(-i) + \mathbb{E}[u(k) \cdot x(k)]
R_{xx}(0) = \sum_{i=1}^p a_i R_{xx}(-i) + \mathbb{E}[u(k) \cdot \left(-\sum_{i=1}^p a_i x(k-i) + u(k)\right)]
自相关函数的递归形式 ：
R_{xx}(\tau) = -\sum_{i=1}^p a_i R_{xx}(\tau - i) + \sigma_u^2 \delta(\tau)
R_{xx}(\tau) = -
\begin{bmatrix}
R_{xx}(0) & R_{xx}(1) & \cdots & R_{xx}(1-p) \\
R_{xx}(1) & R_{xx}(0) & \cdots & R_{xx}(2-p) \\
\vdots & \vdots & \ddots & \vdots \\
R_{xx}(p) & R_{xx}(p-2) & \cdots & R_{xx}(0)
\end{bmatrix}
\begin{bmatrix}
a_1 \\
a_2 \\
\vdots \\
a_p
\end{bmatrix}
噪声方差公式 ：
\sigma_u^2 = R_{xx}(0) + \sum_{i=1}^p a_i \cdot R_{xx}(-i)
\sigma_u^2 = \sum_{i=0}^p a_i \cdot R_{xx}(-i)
测量方程：
y(k) = x(k) + b(k)
自回归模型：
x(k) = -\sum_{k=1}^p a_k x(k-i) + u(k)
输出信号的自相关函数 R_{yy}(\tau) ：
R_{yy}(\tau) = \mathbb{E}[(x(k) + b(k))(x(k-\tau) + b(k-\tau))]
R_{yy}(\tau) = R_{xx}(\tau) + R_{bb}(\tau)
y 的自相关矩阵和x 的自相关矩阵的关系 ：
R_y = R_x + \sigma_b^2 I
R_y^{-1} R_y \neq R_x^{-1} r_x
(R_x - \sigma_b^2 I)^{-1} R_y \neq R_x^{-1} r_x
左右不相等，怎么办呢
一阶MA过程特性: 相关函数在滞后为 0 时，相关性不为零，在滞后为 ±1 时，相关性也不为零，也就是说只有主对角线和两边次对角线受到了影响，因此我们采用噪声补偿，搞定出现的 r_{xx}(0) ，怎么做呢？==> 序号值从 p+1 开始递增：
p+1, , p+2, , \ldots, 2p
自相关与互相关的关系：
\begin{bmatrix}
R_{xx}(p+1) \\
R_{xx}(2p)
\end{bmatrix}
=
\begin{bmatrix}
\hat{R}_{yy}(p+1) \\
\hat{R}_{yy}(2p)
\end{bmatrix}
=
\begin{bmatrix}
R_{xx}(p) & \cdots & R_{xx}(1) \\
R_{xx}(2p+1) & \cdots & R_{xx}(p)
\end{bmatrix}
\begin{bmatrix}
a_1 \\
\vdots \\
a_p
\end{bmatrix}
=
\begin{bmatrix}
\hat{R}_{yy}(p) & \cdots & \hat{R}_{yy}(1) \\
\hat{R}_{yy}(2p+1) & \cdots & \hat{R}_{yy}(p)
\end{bmatrix}
\begin{bmatrix}
a_1 \\
\vdots \\
a_p
\end{bmatrix}
我们称这些为修正后的 Yule-Walker 方程 。因为我们使用了不同的滞后值。因此，从理论上讲，我们可以用这些关系来求解参数 R 。在这种情况下，R_{xx}(0) 并没有参与了计算，所以理论上没有问题。
但在实践中我们无法得到精确的值。当滞后值增加时，样本平均的数量越少，方差会增大，可靠性就越差。
因此，从理论上看，这种方法非常完美，但在实践中，问题在于我们需要越来越大的滞后值。而滞后值越大，在估计自相关函数和自相关矩阵时，需要用更多的样本进行平均才能确保质量。
我们引入了维纳滤波器 ，并且将维纳滤波器应用于自回归过程（AR）
如果我想估计一个 MA（移动平均）过程的参数，首先要先看MA 过程和 AR 过程的不同之处，在频域中，AR 过程会出现共振峰 ，而 MA 过程会出现零点 ，表现为频率的衰减。所以我们可以取信号并绘制它的功率谱，并利用谱下限（lower gramian）来估计一个 MA 过程的参数。(对于一个 MA 过程，功率谱中应该会出现频率的衰减)。然后我们再取功率谱的倒数，也就是谱的逆。这样我们就可以得到共振峰。然后我们还需要与之关联的自相关函数。所以可以进行傅里叶逆变换，从而得到关联的自相关函数。最后重新应用核心计算，就可以得到需要的参数。
自适应滤波器
我们的核心目标是针对时间序列信号进行信号处理，最需要解决的是建模问题，尤其关注于系统的建模，也就是利用一系列方程或者状态空间表达式，来描述信号和系统。对于信号，我们主要关心的是其统计性质，广义平稳的信号（统计特性在时间内不随时间平移而变化）可以大大简化分析，同时需要关注信号的自相关函数(信号与自身在时间间隔上的依赖性)，如果自相关函数快速衰减，表明信号在时间上的相关性较短，易于建模，反之，长记忆过程的建模更加复杂。
我们在之前研究了
均值调整模型（Mean-Adjusted Models）
x(n) = \sum_{j=0}^{q} b_j \cdot u(n - j)
自回归过程（Autoregressive Process, AR）
x_n = a_1 x_{n-1} + a_2 x_{n-2} + \cdots + a_p x_{n-p} + \epsilon_n
自回归与均值调整结合 的过程
白噪声加长记忆特性 的过程
通过添加短相关性达到无限过程 的模型。
然后是是创建线性模型来分析时间序列信号，使用了最简单的模型之一，即 Volterra模型 ，以及最小二次自回归模型 等。
在完成建模部分后，问题就变成了如何处理这些方程，尤其是估计其中的参数。
我们可以将时间序列或信号结构化为矩阵形式，例如 Toeplitz矩阵 或者相关矩阵 （协方差矩阵）。这时，我们要用到一些在文献或学习中见到的线性代数工具，例如QR分解 、LDL分解 、正交投影算子 。还有一种 基于子空间的方法 ，通过分析相关矩阵的特征值 特征向量 和范数 来区分信号子空间 和噪声子空间，利用信号子空间和噪声子空间的正交性来分析信号。
接下来，我们探讨了线性滤波中的块处理问题，线性滤波的目标与传统的高通、低通或带通滤波不同(侧重于频率选择性)，线性滤波对输入信号进行处理，得到尽可能接近期望响应的输出，并最小化均方误差MSE 。我们使用有限冲激响应**（FIR）滤波器**，我们通过反转输入信号的相关矩阵，并乘以输入信号与期望响应的互相关向量来实现。
当信号中存在噪声时，噪声通常被建模为随机过程，其统计特性（如均值、方差、相关性）决定了其对信号和估计的影响，即噪声会显著影响参数估计，因此含噪信号的处理很重要。例如，当渐进模型中有噪声时，我们观察其影响并补偿。
下面我们实现用递归方式逼近解，也就是逐步解决估计问题，即自适应滤波。它是一种实时更新滤波器系数的技术，通过每次接收一个新样本，动态调整滤波器的参数，使得输出信号逐渐逼近期望响应，即随着迭代和样本数量的增加逐渐收敛到精确解。这种方法不需要计算输入信号相关矩阵的反转(块处理方法)，计算量较小，非常适合实时处理场景。
d_k : 期望响应（目标信号）y_k : 滤波器输出，是对 d_k 的估计值e_k = d_k - y_k : 误差，即期望响应和滤波器输出的差异\underline{H}_N(k - 1) 自适应滤波\underline{X}_N(k) : 输入向量，第 k 个样本的信号
我们将研究一个随时间变化的自适应滤波器，也就是说，随着时间的变化以及新样本的到来，这个滤波器将逐渐动态更新滤波器的系数，使其逐渐逼近最优解。
\underline{H}_N(k) = \begin{bmatrix} h_1(k) & \cdots & h_N(k) \end{bmatrix}^T
滤波器的冲激响应是有限的，可以看到这个冲激响应会随时间而变化
\underline{X}_N(i) = \begin{bmatrix} x(i) & \cdots & x(i-N+1) \end{bmatrix}^T
y(k) = \underline{H}_N^T(k-1) \cdot \underline{X}_N(k) = \underline{X}_N^T(k) \cdot \underline{H}_N(k-1)
e(k) = d(k) - \underline{H}_N^T(k-1) \cdot \underline{X}_N(k)
J(k) = \mathbb{E} \left[ e^2(k) \right]
J(k) =\mathbb{E} \left[ (d(k) - \underline{H}_N^T(k-1) \cdot \underline{X}_N(k))^2 \right]
J(k) = \mathbb{E} \left[ d^2(k) \right] - 2 d(k) \cdot \underline{H}_N^T(k-1) \underline{X}_N(k)
+ \underline{H}_N^T(k-1) \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_N(k-1)
现在要做的就是梯度下降法
\underline{H}_N(k) = \underline{H}_N(k-1) - \frac{\alpha}{2} \nabla J
从\underline{H}_N(k-1) 到\underline{H}_N(k) 需要加一个修正项，即步长因子乘负梯度，之所以叫做负梯度方向是因为，这里是减号 ，所以如果梯度是正(斜率为正)，那么还是减，如果梯度是负，说明跑过了，那么就得加一加，往右走走，不停在最小值附近摇摆，最后收敛。 \alpha/2 只是为了简化运算
所以现在要做的就是计算梯度\nabla J
\begin{align}
\nabla J &= \mathbb{E} \left[ -2 \underline{X}_N(k) d(k) + 2 \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_N(k-1) \right] \\\\
&= -2 \mathbb{E} \left[ \underline{X}_N(k) d(k) - \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_N(k-1) \right] \\\\
&= -2 \mathbb{E} \left[ \underline{X}_N(k) \left( d(k) - \underline{X}_N^T(k) \underline{H}_N(k-1) \right) \right] \\\\
&= -2 \mathbb{E} \left[ \underline{X}_N(k) e(k) \right]\\
\end{align}
经过计算后，我们得到了梯度，表达式，我们将其带入原梯度下降法公式中:
\begin{cases}
\underline{H}_N(k) = \underline{H}_N(k-1) + \alpha \mathbb{E} \left[ \underline{X}_N(k) e(k)\right] \\\\
e(k) = d(k) - \underline{X}_N^T(k) \underline{H}_N(k-1)\\
\end{cases}
我们发现\underline{H}_N(k) 有一个问题，实际中我们无法在每个时刻得到期望值，因此我们用瞬时值将其代替，原式可得:
\underline{H}_N(k) = \underline{H}_N(k-1) + \alpha \underline{X}_N(k) e(k)
因此我们获得了LMS算法的核心更新公式，其目的是收敛到这些参数的最小二乘意义下（最小均方算法LMS）的解，即滤波意义上的最优解，这个最优解在理论上 是Yule-Walker方程的解，求这个解不需要时间的演变，可一次直接获得。而自适应方法是利用每一个样本逐步更新参数估计，使得这些估计值最终收敛到基于误差意义上的参数估计，也就是Yule-Walker方程的解。
我们接下来会分析自适应滤波器与最优滤波器之间的偏差，这种偏差目前没有办法完全消除，也就是算法可能在最优解附近波动，但无法完全收敛到最优滤波器的参数。我们可以在某些特定时间点观察参数差异，如果相对较小，我们认可算法可能已经收敛，这个时候我们对后几次迭代参数估计只取均值或者中值(不要只取最后一次迭代的效果，会有瞬时波动的偏差影响)，这样就可以得到更稳健的滤波器系数。补充: 平均值可以平滑，中值可以对抗异常值。
对于其中的 α 学习率（步长因子），控制梯度下降的步幅大小。如果 α 的值很小，梯度下降的过程非常缓慢 ; 如果 α 的值很大，在最优解附近会强烈的振荡，甚至不稳定。
因此我们可以使用Yule-Walker 方程(Yule-Walker 方程在自回归 (AR) 模型中用于描述自相关函数与 AR 系数之间的关系) 验证参数 α 的合理性，其具体步骤如下
先假设一个 α ，然后根据理论AR模型来计算相关函数 r(k)
使用 Yule-Walker 方程，基于计算的 r(k) 和输入信号数据，反向估计 α
比较计算值和先前假设值
随后，我们可以利用 ALMS 滤波器（自适应最小均方滤波器）在不同条件下动态调整步长参数 α ，从而提高滤波器的性能。步长参数的动态调整使得 ALMS 更适合非平稳环境。这里关于ALMS不深入讨论。
我们现在观察自适应滤波器与最优滤波器之间的差异 \Delta \underline{H}_N(k)
\Delta \underline{H}_N(k) = \underline{H}_N(k) - \underline{H}_{\text{opt}}
\underline{H}_{\text{opt}} 是最优滤波器的系数，其表达式如下:
\underline{H}_{\text{opt}} = \mathbb{E} \left\{ \underline{X}_N(k) \underline{X}_N^T(k) \right\}^{-1}
\mathbb{E} \left\{ \underline{X}_N(k) d(k) \right\}
= \mathbf{R}_{xx}^{-1} \mathbf{R}_{xd}
通过分析 \Delta \underline{H}_N(k) ，我们可以评估自适应滤波器在每个时刻与最优滤波器之间的距离，从而查看收敛过程及波动特性，有助于优化算法参数 \alpha ，以实现更快的收敛和更稳定的滤波器性能。
回到前面，我们对\underline{H}_N(k) 有结论:
\underline{H}_N(k) = \underline{H}_N(k-1) + \alpha \underline{X}_N(k) e(k)
这只是第一种形式，反应的是和 e(k) 之间的关系等，现在我们回到它最原始的公式:
\underline{H}_N(k) = \underline{H}_N(k-1) + \alpha \underline{X}_N(k)
\left[ d(k) - \underline{H}_N^T(k-1) \underline{X}_N(k) \right]
我们已知:
\underline{H}_N^T(n-1) \underline{X}_N(n) = \underline{X}_N^T(n) \underline{H}_N(n-1)
因为它是个标量
因此展开:
\underline{H}_N(k) = \underline{H}_N(k-1)
- \alpha \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_N(k-1)
+ \alpha \underline{X}_N(k) d(k)
我们将\underline{H}_N(k) 这个结果带入到 \Delta \underline{H}_N(k) 中:
\begin{align*}
\Delta \underline{H}_N(k) &= \underline{H}_N(k) - \underline{H}_{\text{opt}} \\
&= \underline{H}_N(k-1) - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_N(k-1)
+ \alpha \underline{X}_N(k) d(k) - \underline{H}_{\text{opt}}
\end{align*}
\begin{align*}
\Delta \underline{H}_N(k) &= \left[ I - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \right] \underline{H}_N(k-1)
+ \alpha \underline{X}_N(k) d(k) - \underline{H}_{\text{opt}} \\
&= \left[ I - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \right] \Delta \underline{H}_N(k-1)
+ \alpha \underline{X}_N(k) \left[ d(k) - \underline{X}_N^T(k) \underline{H}_{\text{opt}} \right]
\end{align*}
从公式中可以看出，如果：当迭代在某一时间点 k-1 达到最优滤波器时，后续的迭代可不是最优滤波器，即:
\text{Si } \Delta \underline{H}_N(k-1) = 0 \implies \Delta \underline{H}_N(k) \neq 0
这种情况下我们认为LMS 并未完全收敛到最优解。
我们对 \Delta \underline{H}_N(k) 求期望:
\mathbb{E}[\Delta \underline{H}_N(k)] = \mathbb{E} \left[ \left( I - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \right) \Delta \underline{H}_N(k-1) \right]
+ \alpha \mathbb{E} \left[ \underline{X}_N(k) d(k) - \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_{\text{opt}} \right]
其中:
\mathbb{E} \left[ \underline{X}_N(k) d(k) - \underline{X}_N(k) \underline{X}_N^T(k) \underline{H}_{\text{opt}} \right]
= \mathbf{R}_{dx} - \mathbf{R}_{x} \cdot \mathbf{H}_{\text{opt}}
而我们知道最优滤波器公式为:
\mathbf{H}_{\text{opt}} = \mathbf{R}_{x}^{-1} \underline{\mathbf{R}}_{dx}
一代进去，后面等于0了，因此我们原式有
\mathbb{E}[\Delta \underline{H}_N(k)] = \mathbb{E} \left[ \left( I - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \right) \right]\mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
我们假设独立性，期望值的分解将变得简单:
\mathbb{E}[\Delta \underline{H}_N(k)] = \left( I - \alpha \underline{X}_N(k) \underline{X}_N^T(k) \right) \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
\mathbb{E}[\Delta \underline{H}_N(k)] = \left( I - \alpha \mathbf{R}_{x} \right) \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
我们利用线性代数中的性质
I = P P^{-1}
并且引入 R_x = PDP^\top 的特征分解形式，其中 P 是特征向量矩阵，D 是特征值对角矩阵，则：
\mathbb{E}[\Delta \underline{H}_N(k)] = \left( I - \alpha P D P^{-1} \right) \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
\mathbb{E}[\Delta \underline{H}_N(k)] = P \left( I - \alpha D \right) P^{-1} \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
通过左乘 P^{-1} 化简为：
P^{-1} \mathbb{E} \left[ \Delta \underline{H}_N(k) \right] = \left( I - \alpha D \right) P^{-1} \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right]
最终，这一公式描述了偏差的迭代收敛过程，其收敛速度由 (I - \alpha D) 的特征值及学习率 \alpha 决定。
我们假设:
P^{-1} \mathbb{E} \left[ \Delta \underline{H}_N(k) \right] = \underline{u}_N(k)
P^{-1} \mathbb{E} \left[ \Delta \underline{H}_N(k-1) \right] = \underline{u}_N(k-1)
原式变成:
\underline{U}_N(k)= \left( I - \alpha D \right) \underline{U}_N(k-1)
如果我们用 u_i(k) 表示向量 U_N(k) 的第 i 个分量，那么它满足以下方程：
u_i(k) = (1 - \alpha \lambda_i) u_i(k-1), \quad \forall i \in [1, N]
这是一个一阶差分的标量方程。通过逐步替换，可以从 u_i(0) 表示 u_i(k) ，如下：
u_i(k) = (1 - \alpha \lambda_i)^k u_i(0), \quad \forall i \in [1, N]
因此，证明 LMS 算法在均值上收敛到最优解 \mathbf{H} _{\text{opt}} 可以简化为证明误差 u _i(k) 收敛到零，\forall i \in [1, N] 。因此，\mathbf{R} _{xx} 的特征值 {\lambda_i} _{=1, \dots, N} 是实数且为正，因为 \mathbf{R} _{xx} 是一个对称的正定矩阵。这些特征值必须满足以下收敛条件：
|1 - \alpha \lambda_i| < 1, \quad \forall i \in [1, N]
所以学习率 \alpha 的范围为：
0 < \alpha < \frac{2}{\lambda_i}, \quad \forall i \in [1, N]
因此，LMS 收敛的必要且充分条件涉及输入信号自相关矩阵 \mathbf{R} _{xx} 的最大特征值：
0 < \alpha < \frac{2}{\lambda_{\max}}
只要 \alpha 在这个范围内即可保证了LMS算法在每个方向上的稳定性与最终的收敛性。理论看起来很美好，但是实际中我们必须要先计算相关矩阵，然后进行特征值分解，最后找到最大的特征值，如果找不到\lambda _{\max} 就选不了了，为了解决这个问题，我们提出 NLMS 模型
NLMS（归一化最小均方算法） 在LMS基础上对步长参数进行归一化，可以看作一个带约束的优化问题的解递归算法形式。
其滤波器更新公式为：
\underline{H}_N(k) = \underline{H}_N(k-1) + \frac{\alpha}{\beta + \|\underline{X}_N(k)\|^2} \underline{X}_N(k) e(k)
其中，0 < \alpha < 2 ，正则化标量 \beta 是为了防止向量 \underline{X}_N(k) 取较小值时导致分母为0。
e(k) = d(k) - \underline{X}_N^T(k) \underline{H}_N(k-1)
在实际使用NLMS算法时，需要多次运行并计算平均值 (消除瞬时波动的偏差影响)，即重复运行 NLMS 算法 N次，对所有运行的滤波器系数进行平均，如果只有一次运行，可以将过程分为多个部分，\underline{H}_N(k_1),\underline{H}_N(k_2)\quad ...\quad \underline{H}_N(k_M) ，每一部分的结果视为独立过程，重复步骤计算各段的平均值。
我们还可以对NLMS进一步扩展，变成仿射投影算法，考虑多目标信号和多观测信号的特性，与传统LMS不同，这里期望输出 d(k) 及输入\underline{X}_N (k) 等不再是标量，而是一个向量，其目标是最小化多目标误差的加权平方和 ，以此来调整滤波器，这种扩展使得算法能够处理多个输入信号，适合于高相关性信号或多目标估计场景
卡尔曼滤波器
卡尔曼滤波器发展起来后，因为其通用性和强大的功能，被广泛应用于不同的领域。不仅在信号处理领域可以看到它的应用，同时在自动控制领域也有广泛使用，因为卡尔曼滤波的方法基于系统的状态空间表示。
在我们通常的状态空间表达式建模过程中，测量噪声通常假设为白噪声，其协方差矩阵是已知的，但是实际情况下，噪声可能是彩色的，此时必须将其动态特性归纳到状态向量中进行建模
卡尔曼滤波在状态空间中由两组基本方程组成：状态方程和观测方程。
非线性情况：
x_k = f(x_{k-1}, u_k)
y_k = g(x_k, b_k)
线性情况：
在线性系统中，状态和观测方程具体化为：
x_k = \Phi_{k-1,k} x_{k-1} + G u_k + w_k
观测方程：
y_k = H_k x_k + v_k
其中：
\Phi_{k-1,k} 是状态转移矩阵，
G 是控制输入矩阵，
w_k \sim \mathcal{N}(0, Q) 是高斯分布的过程噪声，协方差为 Q 。
H_k 是观测矩阵，
v_k \sim \mathcal{N}(0, R) 是高斯分布的观测噪声，协方差为 R 。
卡尔曼滤波的核心目标是根据观测值 y_k 和模型动态，优化估计状态 x_k ，同时最小化估计误差的协方差。为了理论方便，我们以简单的线性情况为说明。
参数
先验
\hat{x}_{k|k-1} x_k 的先验估计，状态向量 X_k 在已知 k-1 时刻观测值的情况下的估计
\tilde{x}_{k|k-1}
\tilde{x}_{k|k-1} = x_k - \hat{x}_{k|k-1}
P_{k|k-1}
P_{k|k-1} = \mathbb{E}[\tilde{x}_{k|k-1} \tilde{x}_{k|k-1}^\top]
后验
\hat{x}_{k|k} X_k 在已知 k 时刻观测值的情况下的估计
\tilde{x}_{k|k}
\tilde{x}_{k|k} = x_k - \hat{x}_{k|k}
P_{k|k}
P_{k|k} = \mathbb{E}[\tilde{x}_{k|k} \tilde{x}_{k|k}^\top]
后验估计利用了当前观测值来更新状态，而先验估计则是基于系统模型对下一时刻状态的预测。误差协方差矩阵用于量化估计误差的分布范围。
卡尔曼滤波器的工作过程
卡尔曼滤波器的工作方式与LMS滤波类似，通过一个循环不断更新估计值。首先，在当前时刻 k-1 ，我们从 k-1 时刻的后验状态估计 \hat{x} _{k-1|k-1} 和协方差矩阵 P _{k-1|k-1} 开始，进入预测阶段 。
在没有新观测值 Y _k 的情况下，我们基于系统的动态方程对状态向量进行预测，计算 \hat{X} _{k|k-1} ，即状态向量在时刻 k 的先验估计，同时得到了先验误差协方差矩阵 P _{k|k-1} ，描述预测误差的大小。
更新阶段
当新观测值 Y_k 到达时，我们使用该新信息来修正之前的预测。此时，通过引入一个修正项，将预测值 \hat{X} _{k|k-1} 更新为后验估计 \hat{X} _{k|k} 。修正项的大小依赖于新观测值和预测值之间的差异，以及观测值的不确定性。对应的后验误差协方差矩阵 P _{k|k} 也会被更新，用以反映修正后估计的可靠性
两个步骤的交替运行，完成 k 时刻的预测和更新后，进入 k+1 时刻的预测。可以动态地估计状态向量，并逐步降低估计误差，实际上，预测的核心思想是，预测的误差由协方差矩阵来表征。当有一个新的观测值到来时，我们就可以利用这个新观测值，优化和改进预测的精度。因此，我们会修正预测值，从而得到一个新的状态估计值。这是一个阶段性的过程。
下面我们通过数学公式详细说明这两个过程
关于先验估计 \hat{x}(k|k-1) 和时刻 k-1 的后验估计 \hat{x}(k-1|k-1) 的关系，以及关于 P(k|k-1) 和 P(k-1|k-1) 的关系。
后验估计
\underline{\hat{x}}(k-1 | k-1) = \mathbb{E}[\underline{x}(k-1) | y(1), \dots, y(k-1)] = \int \underline{x} P_{x}(\underline{x} | y(1), \dots, y(k-1)) d\underline{x}
先验估计
\underline{\hat{x}}(k | k-1) = \mathbb{E}[\underline{x}(k) | y(1), \dots, y(k-1)]
\underline{\hat{x}}(k|k-1) = \Phi(k-1, k) \mathbb{E}[\underline{x}(k-1)|y(1), \dots, y(k-1)]
+ G(k) \mathbb{E}[u(k)|y(1), \dots, y(k-1)]
由于模型噪声是白噪声且均值为零，因此：
\mathbb{E}[\underline{u}(k)|y(1), \dots, y(k-1)] = 0
因此，状态向量在时刻 k 的预测可以通过状态向量在时刻 k-1 的估计推导如下：
\underline{\hat{x}}(k|k-1) = \Phi(k-1, k) \underline{\hat{x}}(k-1|k-1)
在这种情况下，仅需要转移矩阵\Phi(k-1, k)
我们定义 \underline{\tilde{x}}(k|k-1) 为时刻 k 的状态向量的先验估计误差 ：
\underline{\tilde{x}}(k|k-1) = \underline{x}(k) - \underline{\hat{x}}(k|k-1)
\underline{\tilde{x}}(k|k-1) = \Phi(k-1, k) \left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right] + G(k) \underline{u}(k)
先验误差的协方差矩阵定义为：
\begin{align*}
P(k|k-1) &= \mathbb{E} \left\{ \left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]
\left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]^T \right\} \\
&= \mathbb{E} \left[ \underline{\tilde{x}}(k|k-1) \underline{\tilde{x}}^T(k|k-1) \right]
\end{align*}
P(k|k-1) = \Phi(k-1, k) P(k-1|k-1) \Phi^T(k-1, k) + G(k) Q(k) G^T(k)
最终得到误差协方差矩阵的更新公式:
\begin{align*}
P(k|k-1) &= \mathbb{E} \left\{ \left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]
\left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]^T \right\} \\
&= \mathbb{E} \left[ \underline{\tilde{x}}(k|k-1) \underline{\tilde{x}}^T(k|k-1) \right]
\end{align*}
现在，我们进入修正阶段，对状态向量的估计进行更新，在这个阶段，我们从上一时刻的估计值出发，并添加一个修正项，它的特别之处在于，它利用了新观测值 \underline{y}(k) 的信息，状态向量线性估计公式如下：
\underline{\hat{x}}(k|k) = \underline{\hat{x}}(k|k-1) + K(k) \left[ \underline{y}(k) - H(k) \underline{\hat{x}}(k|k-1) \right]
K(k) 被称为滤波增益或卡尔曼增益\underline{y}(k) - H(k) \underline{\hat{x}}(k|k-1) 被称为 innovation ，实际观测值 \underline{y}(k) 与其预测值 \underline{\hat{y}}(k) 之间的差异
重新回顾估计误差 \tilde{x}(k|k) 的表达式
\underline{\tilde{x}}(k|k) = \underline{x}(k) - \underline{\hat{x}}(k|k)
并将估计值 \underline{\hat{x}}(k|k) 用上式替换：
\begin{align*}
P(k|k-1) &= \mathbb{E} \left\{ \left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]
\left[ \underline{x}(k) - \underline{\hat{x}}(k|k-1) \right]^T \right\} \\
&= \mathbb{E} \left[ \underline{\tilde{x}}(k|k-1) \underline{\tilde{x}}^T(k|k-1) \right]
\end{align*}
后验误差协方差矩阵:
P(k|k) = \left[ I - K(k) H(k) \right] P(k|k-1) \left[ I - H^T(k) K^T(k) \right] + K(k) R(k) K^T(k)
其中:
K(k) = P(k|k-1) H^T(k) \left[ H(k) P(k|k-1) H^T(k) + R(k) \right]^{-1}
