当前位置：首页 > news >正文

Fast Frequency Estimation Algorithm by Least Squares Phase Unwrapping

news 2025/7/25 7:20:10

I. 引言

单个含噪正弦信号的频率估计是一个研究已久的问题，并有多种应用[1]。在高斯白噪声假设下，最大似然(ML)频率估计器是Rife和Boorstyn [2]中提出的周期图估计器，其中傅里叶变换用于搜索周期图的最大值。周期图估计器被广泛认为是单频估计的最佳方法。在[2]中，粗略搜索基于快速傅里叶变换(FFT)，而精细搜索则基于非线性优化。为了简化非线性优化，[3]中研究了一种五峰值内插器，它可以找到峰值并产生接近ML估计器的估计。此外，信号子空间法(SA)也已用于估计单个含噪正弦信号的频率[4]，其中频率估计在小误差条件下近似无偏，其方差接近ML估计器的方差。

相位展开估计器是另一种单频估计方法，最早在[5]中提出。该方法展开复正弦信号的相位，然后通过线性回归（linear regression，LR）估计频率。然而，传统的相位展开算法在低信噪比(SNR)下会遇到模 $2π2\pi$ 模糊问题，以至于相位展开估计器的阈值处于相对较高的SNR [6]。

为了解决模 $2π2\pi$ 模糊问题，Clarkson将相位展开问题转化为寻找格中的最近点，并利用格点理论解决了该问题[7]。仿真结果表明，其估计性能接近周期图估计器，但作者没有提供理论证明。最近，McKilliam在最小二乘意义上发展了基于格的方法[8]，并证明了该估计器是强一致的，其估计方差在数值模拟中的高信噪比下接近克拉默-拉奥界(CRB) [2]。在[8]中，最小二乘相位展开估计器(LSPUE)显示出与经典周期图估计器相同的阈值。然而，基于格点理论的相位展开具有高达 $O(N^3\log N)$ 的计算复杂度，因此不能用于实际应用。最近，有人提出了一种基于相位滤波器的快速相位展开方法[9]，但是在频率估计问题中滤波器参数很难确定。

在这封信中，我们提出了一种快速频率估计算法，它结合了LSPUE的准确性和滤波方法的快速操作。所提出的方法在最小二乘意义上执行相同的相位展开，因此它具有与[8]中方法相同的性能，但计算复杂度低至 $O(N^2)$ 。

II. PROBLEM DESCRIPTION

假设观测到的包含 $N$ 个样本的序列具有以下形式：

$r_n = Ae^{j2\pi(fn+\theta)} + w_n, \quad n=0, 1, \dots, N-1 \quad \tag{1}$

其中 $f$ 是频率， $θ\theta$ 是初始相位， $A$ 是信号幅度。噪声样本 $w_n$ 是独立的、同分布的高斯复随机变量，均值为0，方差为 $2σ22\sigma^2$ 。在这封信中，我们试图估计频率 $f$ ，并假设 $f$ 在 $[- 1/2, 1/2)$ 范围内。 $r_n$ 的相位由下式给出：

$\angle r_n = 2\pi(fn + \theta + v_n) \quad (\text{mod } 2\pi) \quad \tag{2}$

其中 $v_n$ 是在 $[- 1/2, 1/2)$ 范围内定义的相位噪声。噪声项 $v_n$ 的分布如下[10]：

$f_v(x) = e^{-\gamma} + 2\sqrt{\pi\gamma}\cos{2\pi x}e^{-\gamma\sin^2{2\pi x}}Q(-\sqrt{2\gamma}\cos{2\pi x}) \quad \tag{3}$

其中 $γ=A22σ2\gamma = \frac{A^2}{2\sigma^2}$ 是信噪比(SNR)。Q函数定义为：

$\int_x^\infty \phi(t)dt \quad \tag{4}$

其中

$\phi(t) = \frac{1}{\sqrt{2\pi}}e^{-\frac{t^2}{2}} \quad \tag{5}$

是标准正态分布。对于 $x > 0$ ，我们有

$\int_x^\infty \phi(t)dt < \int_x^\infty \frac{t}{x}\phi(t)dt = \frac{\phi(x)}{x}$

$du=−t⋅dt)=1x2π∫−∞−x2/2eudu=1x2π[eu]−∞−x2/2=1x2π(e−x2/2−lim⁡a→−∞ea)=1x2π(e−x2/2−0)=1x(12πe−x2/2)=ϕ(x)x\begin{align*} \int_x^\infty \frac{t}{x}\phi(t)dt &= \frac{1}{x} \int_x^\infty t \cdot \phi(t) dt \\ &= \frac{1}{x} \int_x^\infty t \cdot \frac{1}{\sqrt{2\pi}}e^{-\frac{t^2}{2}} dt \\ &= \frac{1}{x \sqrt{2\pi}} \int_x^\infty t e^{-\frac{t^2}{2}} dt \\ &= \frac{1}{x \sqrt{2\pi}} \int_{-x^2/2}^{-\infty} e^u (-du) \quad (\text{令 } u = -t^2/2, \text{ 则 } du = -t \cdot dt) \\ &= \frac{1}{x \sqrt{2\pi}} \int_{-\infty}^{-x^2/2} e^u du \\ &= \frac{1}{x \sqrt{2\pi}} \left[ e^u \right]_{-\infty}^{-x^2/2} \\ &= \frac{1}{x \sqrt{2\pi}} \left( e^{-x^2/2} - \lim_{a \to -\infty} e^a \right) \\ &= \frac{1}{x \sqrt{2\pi}} \left( e^{-x^2/2} - 0 \right) \\ &= \frac{1}{x} \left( \frac{1}{\sqrt{2\pi}} e^{-x^2/2} \right) \\ &= \frac{\phi(x)}{x} \end{align*}$

并且

$\left(1+\frac{1}{x^2}\right)Q(x) > \int_x^\infty \left(1+\frac{1}{t^2}\right)\phi(t)dt = \frac{\phi(x)}{x}.$

$∫x∞(1+1t2)ϕ(t)dt=∫x∞(−ddt(ϕ(t)t))dt=−[ϕ(t)t]x∞=−(lim⁡t→∞ϕ(t)t−ϕ(x)x)=−(lim⁡t→∞e−t2/2t2π−ϕ(x)x)=−(0−ϕ(x)x)=ϕ(x)x\begin{align*} \int_x^\infty \left(1+\frac{1}{t^2}\right) \phi(t) dt &= \int_x^\infty \left(-\frac{d}{dt}\left(\frac{\phi(t)}{t}\right)\right) dt \\ &= -\left[ \frac{\phi(t)}{t} \right]_x^\infty \\ &= -\left( \lim_{t\to\infty} \frac{\phi(t)}{t} - \frac{\phi(x)}{x} \right) \\ &= -\left( \lim_{t\to\infty} \frac{e^{-t^2/2}}{t\sqrt{2\pi}} - \frac{\phi(x)}{x} \right) \\ &= -\left( 0 - \frac{\phi(x)}{x} \right) \\ &= \frac{\phi(x)}{x} \end{align*}$

因此，

$\frac{x}{1+x^2}\phi(x) < Q(x) < \frac{\phi(x)}{x}, \quad x>0. \quad \tag{6}$

考虑到 $Q (x) = 1 - Q (- x)$ ，我们有

$1+\frac{\phi(x)}{x} < Q(x) < 1+\frac{x}{1+x^2}\phi(x), \quad x<0. \quad \tag{7}$

将(6)和(7)代入(3)，我们发现对于较大的 $γ\gamma$ ， $f_v(x)$ 的下界和上界近似相等。在这种情况下，

$f_v(x) = \begin{cases} 2\sqrt{\pi\gamma}\cos{2\pi x}e^{-\gamma\sin^2{2\pi x}}, & |x| \le 1/4 \\ 0, & |x| > 1/4. \end{cases} \quad \tag{8}$

令 $xn=∠rn/(2π)x_n = \angle r_n / (2\pi)$ ，我们得到

$x_n = fn + \theta + v_n - u_n \quad \tag{9}$

其中 $un=[fn+θ+vn]u_n = [fn + \theta + v_n]$ ， $[x]$ 表示与 $x$ 最接近的整数。因此，我们有

$v_n = x_n + u_n - fn - \theta. \quad \tag{10}$

为了在最小二乘意义上使噪声项 $v_n$ 最小，我们定义目标函数为

$\theta, u_n) = \sum_{n=0}^{N-1}(x_n + u_n - fn - \theta)^2 \quad \tag{11}$

III. FAST LSPUE ALGORITHM

为了最小化 $\theta, u_n)$ ，困难在于 $u_n$ 是未知的。如果 $u_n$ 已知，我们可以很容易地用最小二乘法公式解决这个问题。在频率估计问题中，我们发现 $u_n$ 是 $n$ 的单调函数。具体来说，当 $f > 0$ 时， $u_n$ 是单调递增函数，当 $f < 0$ 时， $u_n$ 是单调递减函数。考虑到 $\in [-1/2, 1/2)$ ，我们有 $−N2≤f(N−1)+θ<N2-\frac{N}{2} \le f(N-1)+\theta < \frac{N}{2}$ 。因此，我们定义

$u_n^m = \left\lceil m\frac{n}{N} \right\rceil, \quad m = -\frac{N}{2}, -\left(\frac{N}{2}-1\right), \dots, \frac{N}{2}-1. \quad \tag{17}$

对于给定的 $m$ ， $um=[u0m,u1m,…,uN−1m]T\boldsymbol u^m = [u_0^m, u_1^m, \dots, u_{N-1}^m]^T$ 是一个在整数空间 $ZN\mathbb Z^N$ 中的向量。显然， $um\boldsymbol u^m$ 对应于频率为 $mN\frac{m}{N}$ 的信号。我们定义 $N$ 个向量，它们代表了 $ZN\mathbb Z^N$ 中最有可能的 $N$ 个解。我们首先检查这 $N$ 个解，然后在这些 $N$ 个解附近寻找更好的解。整个过程类似于[2]中描述的粗略搜索和精细搜索，但我们是在展开的相位中搜索最佳的 $u_n$ ，而[2]中的算法是在频域中搜索最大值。

现在，我们重新定义目标函数

$J^m(f, \theta) = \sum_{n=0}^{N-1}(y_n^m - fn - \theta)^2 \quad \tag{18}$

其中

$y_n^m = x_n + u_n^m \quad \tag{19}$

是 $x_n$ 的展开相位。由于 $u_n^m$ 是已知的，由[11]可得频率和相位

$\frac{12 \sum_{n=0}^{N-1}(n - (N-1)/2)y_n^m}{N(N-1)(N+1)} \quad \tag{20}$

$\theta(m) = \frac{2 \sum_{n=0}^{N-1}(2N-3n-1)y_n^m}{N(N+1)}.$

残差误差表示为 $Jm(f(m),θ(m))J^m(f(m), \theta(m))$ 。这是粗略搜索过程。显然， $f (m)$ 和 $θ(m)\theta(m)$ 是给定 $um\boldsymbol u^m$ 的最佳解，但它们不是全局最优解。现在，我们尝试在 $um\boldsymbol u^m$ 附近找到更好的解。这是精细搜索过程。首先，我们重构线性相位

$p_n = f(m)n + \theta(m). \quad \tag{21}$

然后，我们围绕线性相位 $p_n$ 展开相位 $x_n$ ，如下所示

$y^nm=pn+M[xn−pn](22) \hat{y}_n^m = p_n + M[x_n - p_n] \quad \tag{22}$

其中

$\begin{cases} M[X-1], & X > 1/2 \\ M[X+1], & X < -1/2 \\ X, & \text{otherwise}. \end{cases}$

在[9]中研究了(22)中的相位展开操作。这样， $y^nm\hat{y}_n^m$ 和 $x_n$ 满足 $y^nm=xn(mod1)\hat{y}_n^m = x_n \pmod 1$ 和 $∣y^nm−pn∣<1/2|\hat{y}_n^m - p_n| < 1/2$ ，这意味着 $y^nm\hat{y}_n^m$ 是 $x_n$ 围绕 $p_n$ 的相位展开。现在，我们重新定义目标函数

$J^m(f,θ)=∑n=0N−1(y^nm−fn−θ)2.(23) \hat{J}^m(f, \theta) = \sum_{n=0}^{N-1}(\hat{y}_n^m - fn - \theta)^2. \quad \tag{23}$

同样，我们通过最小二乘法公式求解频率和相位

$f^(m)=12∑n=0N−1(n−(N−1)/2)y^nmN(N−1)(N+1)(24) \hat{f}(m) = \frac{12 \sum_{n=0}^{N-1}(n - (N-1)/2)\hat{y}_n^m}{N(N-1)(N+1)} \quad \tag{24}$

$θ^(m)=2∑n=0N−1(2N−3n−1)y^nmN(N+1). \hat{\theta}(m) = \frac{2 \sum_{n=0}^{N-1}(2N-3n-1)\hat{y}_n^m}{N(N+1)}.$

残差误差表示为 $J^m(f^(m),θ^(m))\hat{J}^m(\hat{f}(m), \hat{\theta}(m))$ 。现在，我们比较残差误差 $Jm(f(m),θ(m))J^m(f(m), \theta(m))$ 和 $J^m(f^(m),θ^(m))\hat{J}^m(\hat{f}(m), \hat{\theta}(m))$ 。如果

$Jm(f(m),θ(m))≤J^m(f^(m),θ^(m)) J^m(f(m), \theta(m)) \le \hat{J}^m(\hat{f}(m), \hat{\theta}(m))$

则 $Jm(f(m),θ(m))J^m(f(m), \theta(m))$ 是给定 $m$ 的局部最优解，精细搜索结束。如果

$Jm(f(m),θ(m))>J^m(f^(m),θ^(m)) J^m(f(m), \theta(m)) > \hat{J}^m(\hat{f}(m), \hat{\theta}(m))$

我们按如下方式更新参数：

$f(m)=f^(m)θ(m)=θ^(m)Jm(f(m),θ(m))=J^m(f^(m),θ^(m)) \begin{aligned} f(m) &= \hat{f}(m) \\ \theta(m) &= \hat{\theta}(m) \\ J^m(f(m), \theta(m)) &= \hat{J}^m(\hat{f}(m), \hat{\theta}(m)) \end{aligned}$

然后重复从(21)到(24)的步骤。

同样地，我们计算所有 $m$ 的残差误差 $Jm(f(m),θ(m))J^m(f(m), \theta(m))$ ，并按如下方式找到最优的 $m$ ：

$m^* = \arg \min J^m(f(m), \theta(m)). \quad \tag{25}$

图1显示了整个算法的示意图。图1(a)显示了粗略搜索的流程，图1(b)显示了精细搜索的流程。显然，粗略搜索中有 $N$ 个循环，每个循环都包含一次用于寻找局部最优解的精细搜索。根据我们的数值模拟结果，每次精细搜索包含的迭代过程通常在五到八轮后收敛。在每次迭代中，相位展开和最小二乘计算需要 $O (N)$ 次算术运算。因此，整个FLSPUE算法的计算复杂度为 $O(N^2)$ 。

IV. SIMULATION RESULTS

在本节中，我们检验在第三节中开发的FLSPUE算法的性能。由于所提出的FLSPUE算法是LSPUE方案的一种快速实现，其性能应与[8]中的LSPUE相同。因此，我们仅比较了五种估计器的性能：ML估计器[2]、快速近似ML估计器（fast ML）[3]、SA[4]、LR[5]以及所提出的FLSPUE算法。理想的ML估计器是搜索周期图的最大值，但[2]中的非线性优化很复杂。在本信中，我们通过chirp z变换（CZT）算法[12]执行精细搜索，该算法计算由粗略搜索找到的峰值附近的离散傅里叶变换。然后，我们可以从CZT算法的计算结果中搜索最大值。CZT算法是基于FFT实现的，因此基于CZT的ML估计器的计算复杂度仍然是 $\log N$ 。

在这里插入图片描述

在我们的仿真中，单正弦信号的参数为 $f = 0.1$ 和 $θ=0\theta=0$ 。当 $N = 64$ 和 $N = 256$ 时的仿真结果显示在图2和图3中，每次仿真的信噪比(SNR)在–20和20 dB之间变化，每个SNR值都运行了10 000次试验。在图2和图3中，基于CZT的ML估计器被称为CZT ML。在SA算法中，观测信号表示为一个 $\times n$ 的数据矩阵，且SA算法在不同矩阵维度下的性能是不同的。在这封信中，我们为 $N = 64$ 选择 $\times 8$ 的矩阵，为 $N = 256$ 选择 $16 \times 16$ 的矩阵，因为使用这些参数可以获得最佳的阈值性能。