Fast Fundamental Frequency Estimation: Making a Statistically Efficient Estimator Computationally Efficient，但注重的并不是这篇论文的核心内容，而是它引述和整理此前工作的部分，此篇论文提出的新观点是对以前的工作进行优化，这点我们没有涉及。

A 信号模型与损失函数

A1 信号模型

• 对基频 $f_0$ 用数字角频率代表 $\omega_0 = \frac {2 \pi f_0}{f_s}$，并对由 $L$ 个整数倍谐波成分（频率分别为 $f_0, 2f_0, \cdots, L f_0$）叠加而成的波形进行这样的建模：

  $$ x[n] = \sum_{k=1}^L (a_k \cos (k \omega_0 n) - b_k \sin (k \omega_0 n)) + e[n], \\n \in [n_0, n_0 + N), \ \ \ \omega_0 \in [0, \pi) $$

  • 其中 $a_k, b_k$ 是第 $k$ 个谐波成分的系数（根据辅助角公式可知此谐波的幅度 $\sqrt{ a^2_k + b^2_k }$ ，正弦初相 $\phi_0$ 满足 $\tan \phi_0 = -\frac {a_k}{b_k}$）。
  • $e(n)$ 为误差，期望它是一个加性噪声，符合高斯分布。
  • $\omega_0 \in [0, \pi)$ 的考量自然是 $0 < f_0 < \frac 12 f_s$。

• 可以把这个信号模型写成矩阵表示：

  $$ \boldsymbol x = \boldsymbol Z(\omega_0) \boldsymbol \alpha + \boldsymbol e $$

  • 其中 $\boldsymbol Z$ 只由 $\omega$ 决定，而 $\boldsymbol \alpha$ 表示的是所有的系数：

    $$ \begin{align*}

    \boldsymbol Z(\omega) &=

    \begin{bmatrix}

    \cos \omega n_0 & \cos 2\omega n_0 & \cdots & \cos L\omega n_0 &

    \sin \omega n_0 & \sin 2\omega n_0 & \cdots & \sin L\omega n_0

    \\

    \cos \omega (n_0+1) & \cos 2\omega (n_0+1) & \cdots & \cos L\omega (n_0+1) &

    \sin \omega (n_0+1) & \sin 2\omega (n_0+1) & \cdots & \sin L\omega (n_0+1)

    \\

    \vdots & \vdots & & \vdots & \vdots & \vdots & & \vdots\\

    \cos \omega (n_0+N-1) & \cos 2\omega (n_0+N-1) & \cdots & \cos L\omega (n_0+N-1) &

    \sin \omega (n_0+N-1) & \sin 2\omega (n_0+N-1) & \cdots & \sin L\omega (n_0+N-1)

    \end{bmatrix} \\

    \boldsymbol \alpha &= \begin{bmatrix} a_1 \\ a_2 \\ \vdots \\ a_L \\ -b_1 \\ -b_2 \\ \vdots \\ -b_L \end{bmatrix}

    \end{align*} $$

A2 损失函数

• 自然引出损失函数 $L(\omega_0, \boldsymbol \alpha)$ 和优解 $(\omega_0^, \boldsymbol \alpha^)$：

  $$ L(\omega_0, \boldsymbol \alpha) = \|\boldsymbol e\|^2 = \|\boldsymbol x - \boldsymbol Z(\omega_0) \boldsymbol\alpha\|^2, \\ (\omega_0^, \boldsymbol \alpha^) = \argmin_{\omega_0\in [0, \pi), \ \boldsymbol \alpha \in \R^{2L}} L(\omega_0, \boldsymbol \alpha) $$

B 基础求解方式

B1 优化 $\boldsymbol \alpha$ 的部分

• 有全微分加链式法则求解（这里假定 $\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0)$ 满秩即 $\boldsymbol Z(\omega_0)$ 列满秩，显然至少要求 $N \ge 2L$）：

  

  REMINDER

  • 离谱的链式法则

    • 如果我有一个变量 $\omega$，通过一个标量变元矩阵函数 $\boldsymbol z = \boldsymbol Z(\omega)$，然后再通过一个矩阵变元标量函数 $f(\boldsymbol z)$ 得到 $v = f(\boldsymbol Z(\omega))$，$v, \boldsymbol z, \omega$ 三变量之间的链式法则应该描述如下：

      $$ \frac{dv}{d\omega} = {\rm tr} \left( \left( \frac{\partial f(\boldsymbol z)}{\partial \boldsymbol z} \right)^{\rm T} \frac{d\boldsymbol Z(\omega)}{d\omega} \right) $$

• 可以知道当给定 $\omega_0$ 时，最优的系数一定是 $\boldsymbol \alpha^*(\omega_0) = (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0) \boldsymbol x$。

• 剩下的问题就是求解 $\omega_0$，但它与损失函数之间是一个非线性关系，求导也没用，$\omega_0$ 并不会在某阶导消失。只好暂且把最优的 $\boldsymbol \alpha^*(\omega_0)$ 回代到 $L$ 式子中，得到 $L_1(\omega_0)$ 为新的最小化对象：

  $$ \begin{align*}

  L_1(\omega_0) &= L(\omega_0, \boldsymbol \alpha^*(\omega_0)) \\ &= \|\boldsymbol x - \boldsymbol Z(\omega_0) (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0) \boldsymbol x\|^2 \\

  \end{align*} $$

  折腾了一会。以为可化简成下面这样。

  $$ \| \boldsymbol x \|^2 \| \boldsymbol I - \boldsymbol Z(\omega_0) (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0) \| ^2 $$

  然后才想起，对于相容的矩阵范数和向量范数（这里是矩阵二范数和向量二范数），仅约束一个上界：$\| \boldsymbol A \boldsymbol x \| _2 \le \| \boldsymbol A \|_2 \|\boldsymbol x\|_2$，对约束最小化有一定帮助，但不大，尤其是如果这么做，$\| \boldsymbol I - \boldsymbol Z(\omega_0) (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0) \| ^2$ 的部分就和具体信号 $\boldsymbol x$ 没关系了。

  具体请参见：

  [第A4章] 向量范数、矩阵范数与摄动分析

B2 继续化简 omega_0 的损失函数

• 记 $\boldsymbol Y(\omega_0) = \boldsymbol Z(\omega_0) (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0)$，可知这个形式是一个投影矩阵，即满足对称性 $\boldsymbol Y^{\rm T} = \boldsymbol Y$ 和等幂性 $\boldsymbol Y^2 = \boldsymbol Y$，而 $(\boldsymbol I - \boldsymbol Y)$ 就是其补投影矩阵，也是投影矩阵、也具有这些性质。

  • 关于为什么这个形式是投影矩阵、投影矩阵及补投影矩阵等具体的讨论请见：

    [第E3章] 正交投影矩阵

• 因而可以通过下面方法展开 $L_1(\omega_0)$：

  $$ \begin{align*}

  L_1(\omega_0) &= \| (\boldsymbol I - \boldsymbol Y(\omega_0)) \boldsymbol x \|^2 \\

  &= \boldsymbol x^{\rm T} (\boldsymbol I - \boldsymbol Y(\omega_0))^{\rm T} (\boldsymbol I - \boldsymbol Y(\omega_0)) \boldsymbol x \\

  &= \boldsymbol x^{\rm T}(\boldsymbol I - \boldsymbol Y(\omega_0)) \boldsymbol x \\

  &= \boldsymbol x^{\rm T} \boldsymbol x - \boldsymbol x^{\rm T} \boldsymbol Y(\omega_0) \boldsymbol x

  \end{align*} $$

• 因而作 $L_2(\omega_0) = \boldsymbol x^{\rm T} \boldsymbol Y(\omega_0) \boldsymbol x$，当 $\boldsymbol x$ 给定的时候，优化问题发生如下转换：

  $$ \argmin_{\omega \in [0, \pi)} L_1(\omega_0) = \argmax_{\omega \in [0, \pi)} L_2(\omega_0) $$

  • 请注意，虽然 $L_2(\omega_0) = \boldsymbol x^{\rm T} \boldsymbol Y(\omega_0) \boldsymbol x$ 是一个夹心糖式，但这里的变数是 $\boldsymbol Y$，而不是常见的 $\boldsymbol x$。

• 与此同时，如果是先得到最优系数 $\boldsymbol \alpha^(\omega_0) = (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0) \boldsymbol x$，还能通过 $L_2(\omega_0) = \boldsymbol x^{\rm T} \boldsymbol Z(\omega_0) \boldsymbol \alpha^(\omega_0)$ 得到 $L_2(\omega_0)$。

B3 基础算法

算法

枚举 $l = 1, 2, \cdots, L$。

通过某种方法枚举 $\omega_0 \in (0, \frac \pi l)$ （见下面的说明）。

计算 $\boldsymbol Z (\omega_0)$。

计算 $\boldsymbol Y(\omega_0) = \boldsymbol Z(\omega_0) (\boldsymbol Z^{\rm T} (\omega_0) \boldsymbol Z(\omega_0) )^{-1} \boldsymbol Z^{\rm T} (\omega_0)$。

计算 $L_2(\omega_0) = \boldsymbol x^{\rm T} \boldsymbol Y(\omega_0) \boldsymbol x$。