带约束的高斯牛顿法求解多音信号分离问题

一、信号模型与优化问题建立

1. 复信号模型

设观测的复信号由两个单频复指数信号加噪声组成：

$$x[n]=A_0{e}^{j(2\pi f_{0}n{T}_s+{\varphi }_0)}+A_1{e}^{j(2\pi f_{1}n{T}_s+{\varphi }_1)}+w[n],n=0,1,\dots ,N-1$$

其中：

• $A_0,A_1>0$ 为幅度

• $f_0,f_1\in (0,f_s/2)$ 为频率

• ${\varphi }_0,{\varphi }_1\in [-\pi ,\pi ]$ 为相位

• $T_s=1/f_s$ 为采样间隔

• $w[n]\sim \mathcal{C}\mathcal{N}(0,{\sigma }^2)$ 为复高斯噪声

2. 参数向量定义

将待估计参数定义为向量形式：

$$\mathbf{\theta }={\left[A_0,f_0,{\varphi }_0,A_1,f_1,{\varphi }_1\right]}^T\in {\mathbb{R}}^6$$

3. 优化目标

通过最小二乘准则估计参数：
$$\underset{\mathbf{\theta }}{\min }J(\mathbf{\theta })=\sum _{n=0}^{N-1}{\left|x[n]-s[n;\mathbf{\theta }]\right|}^2$$
其中信号模型为：

$$s[n;\mathbf{\theta }]=A_0{e}^{j(2\pi f_{0}n{T}_s+{\varphi }_0)}+A_1{e}^{j(2\pi f_{1}n{T}_s+{\varphi }_1)}$$

4. 约束条件

根据物理意义添加边界约束：

$$\{\begin{array}{l}0\le A_k\le A_{\max }\\ f_{\min }\le f_k\le f_{\max }\\ -\pi \le {\varphi }_k\le \pi \end{array},k=0,1$$

---

二、优化问题求解推导

1. 目标函数的复变函数处理

由于目标函数是复值，需用Wirtinger微积分求导。定义残差：

$$r_n(\mathbf{\theta })=x[n]-s[n;\mathbf{\theta }]$$

则目标函数可写为：

$$J(\mathbf{\theta })=\sum _{n=0}^{N-1}{r}_n(\mathbf{\theta })\overline{r_n(\mathbf{\theta })}$$

其中$\overline{(\cdot )}$表示复共轭。

2. 梯度计算（Wirtinger导数）

Wirtinger导数定义为：
$${\nabla }_{\mathbf{\theta }}J=2\sum _{n=0}^{N-1}\mathrm{Re}\left({\frac{\partial r_n}{\partial \mathbf{\theta }}}^H{r}_n\right)$$

其中$(\cdot {)}^H$表示共轭转置。具体到每个参数：

• 幅度$A_0$的梯度分量：

$$\frac{\partial r_n}{\partial A_0}=-e^{j(2\pi f_{0}n{T}_s+{\varphi }_0)}$$

$$\frac{\partial J}{\partial A_0}=2\mathrm{Re}\left(\sum _{n=0}^{N-1}\left[-e^{-j(2\pi f_{0}n{T}_s+{\varphi }_0)}\right]r_n\right)$$

• 频率$f_0$的梯度分量：

$$\frac{\partial r_n}{\partial f_0}=-j\cdot 2\pi n{T}_s{A}_0{e}^{j(2\pi f_{0}n{T}_s+{\varphi }_0)}$$

$$\frac{\partial J}{\partial f_0}=2\mathrm{Re}\left(\sum _{n=0}^{N-1}\left[j\cdot 2\pi n{T}_s{A}_0{e}^{-j(2\pi f_{0}n{T}_s+{\varphi }_0)}\right]r_n\right)$$

• 相位${\varphi }_0$的梯度分量：

$$\frac{\partial r_n}{\partial {\varphi }_0}=-j{A}_0{e}^{j(2\pi f_{0}n{T}_s+{\varphi }_0)}$$

$$\frac{\partial J}{\partial {\varphi }_0}=2\mathrm{Re}\left(\sum _{n=0}^{N-1}\left[j{A}_0{e}^{-j(2\pi f_{0}n{T}_s+{\varphi }_0)}\right]r_n\right)$$

（$A_1,f_1,{\varphi }_1$的梯度形式类似，只需将下标0改为1）

3. Hessian矩阵近似

为加速收敛，采用Gauss-Newton法近似Hessian：

$$\mathbf{H}(\mathbf{\theta })\approx 2\sum _{n=0}^{N-1}\mathrm{Re}\left({\frac{\partial r_n}{\partial \mathbf{\theta }}}^H\frac{\partial r_n}{\partial \mathbf{\theta }}\right)$$

其中雅可比矩阵的行向量为：

$${\mathbf{J}}_n=\frac{\partial r_n}{\partial \mathbf{\theta }}=\left[\frac{\partial r_n}{\partial A_0},\frac{\partial r_n}{\partial f_0},\frac{\partial r_n}{\partial {\varphi }_0},\frac{\partial r_n}{\partial A_1},\frac{\partial r_n}{\partial f_1},\frac{\partial r_n}{\partial {\varphi }_1}\right]$$

4. 带约束的Gauss-Newton算法

迭代格式如下：

步骤1：初始化参数${\mathbf{\theta }}^{(0)}$，设置迭代索引$k=0$

步骤2：计算当前残差$r_n^{(k)}=x[n]-s[n;{\mathbf{\theta }}^{(k)}]$

步骤3：计算梯度${\mathbf{g}}^{(k)}=\nabla J({\mathbf{\theta }}^{(k)})$和近似Hessian${\mathbf{H}}^{(k)}$

步骤4：求解带约束的线性最小二乘问题：

$${\mathbf{\delta }}^{(k)}=\arg \underset{\mathbf{\delta }}{\min }\Vert {\mathbf{J}}^{(k)}\mathbf{\delta }+{\mathbf{r}}^{(k)}{\Vert }^2$$

$$\text{s.t.}{\mathbf{\theta }}^{(k)}+\mathbf{\delta }\in \Omega$$

其中$\Omega$为参数可行域，${\mathbf{r}}^{(k)}=[r_0^{(k)},\dots ,r_{N-1}^{(k)}{]}^T$

步骤5：更新参数${\mathbf{\theta }}^{(k+1)}={\mathbf{\theta }}^{(k)}+\mu {\mathbf{\delta }}^{(k)}$（$\mu$为步长）

步骤6：若$\Vert {\mathbf{\delta }}^{(k)}\Vert <ϵ$则停止，否则$k\leftarrow k+1$转步骤2

---

三、边界约束处理策略

采用投影梯度法处理边界约束：

1. 在每次迭代更新后，检查参数是否越界：

$${\mathbf{\theta }}_{\text{new}}={\mathcal{P}}_{\Omega }\left(\mathbf{\theta }+\mu \mathbf{\delta }\right)$$

其中投影算子${\mathcal{P}}_{\Omega }$定义为：

$${\mathcal{P}}_{\Omega }({\theta }_i)=\{\begin{array}{ll}{\theta }_{\min ,i}& {\theta }_i<{\theta }_{\min ,i}\\ {\theta }_i& {\theta }_{\min ,i}\le {\theta }_i\le {\theta }_{\max ,i}\\ {\theta }_{\max ,i}& {\theta }_i>{\theta }_{\max ,i}\end{array}$$

2. 对于相位周期性约束，投影后需归一化到$[-\pi ,\pi ]$：
   $${\varphi }_{\text{proj}}=\mathrm{mod}(\varphi +\pi ,2\pi )-\pi$$

---

四、算法完整流程（伪代码）

输入：观测信号 x[0..N-1], 采样率 fs, 初始估计 θ_init输出：优化参数 θ_opt设定：最大迭代次数 K, 容差 ε, 步长 μ=0.1θ_prev ← θ_initfor k = 1 to K do// 1. 计算当前残差和雅可比矩阵for n = 0 to N-1 dor_n = x[n] - s(n; θ_prev)J_n = [∂r/∂A0, ∂r/∂f0, ∂r/∂φ0, ∂r/∂A1, ∂r/∂f1, ∂r/∂φ1]  // 按Wirtinger导数公式end// 2. 构造正规方程H = Re(J^H * J)  // 6×6实对称矩阵g = Re(J^H * r)  // 6×1实向量// 3. 求解线性系统 (带约束)δ = H \ g  // 高斯消去法或Cholesky分解// 4. 带投影的参数更新θ_temp = θ_prev + μ * δθ_new = ProjectToBounds(θ_temp)  // 边界投影// 5. 收敛检查if ||θ_new - θ_prev|| < ε thenbreakendθ_prev ← θ_newendθ_opt ← θ_new

---

五、与现有方法的理论对比

| 算法 | 计算复杂度 | 收敛速度 | 全局最优保证 | 边界处理 |

|-----------------|-------------|---------|------------|---------|

| 本文方法 | O(6^2N) | 二次收敛 | 局部最优 | 投影法 |

| SSA-BSS | O(N log N) | 一次迭代 | 依赖初始化 | 无 |

| 粒子滤波 | O(NM) | 渐近收敛 | 概率保证 | 重采样 |

| 循环对消 | O(N log N) | 一次迭代 | 无 | 无 |

注：M为粒子数

---

六、论文书写建议

1. 问题建模部分：

• 明确定义复信号模型和参数向量

• 给出完整的最小二乘目标函数和约束条件

2. 算法推导部分：

• 详细说明Wirtinger导数的推导过程

• 给出Gauss-Newton法的矩阵形式更新公式

• 解释边界投影算子的实现

3. 实验部分：

• 比较梯度解析计算与数值微分的精度差异

• 展示不同SNR下的参数估计Cramér-Rao界

4. 附录：

• 提供核心算法的伪代码

• 补充投影梯度的收敛性证明

通过以上数学化表达，可避免出现MATLAB工具箱依赖，提升论文的理论严谨性。实际实现时，可基于上述推导自主编写优化器（如用C++或Python实现），无需调用现成优化库。

取前两个峰值位置 ${\hat{f}}_0,{\hat{f}}_1$

@article{chen2022,
title={Precision Extraction of Weak Harmonic Signals in Strong Interference Environments},
author={Chen, Y. and Wang, L. and Smith, J.},
journal={IEEE Transactions on Instrumentation and Measurement},
volume={71},
pages={1–10},
year={2022},
doi={10.1109/TIM.2022.3147321}
}
其中 $\mathbf{\theta }=[A_0,f_0,{\varphi }_0,A_1,f_1,{\varphi }_1{]}^T$ 是 6维参数向量

$\mathbf{\delta }=-({\mathbf{J}}^T\mathbf{J}{)}^{-1}{\mathbf{J}}^T\mathbf{r}$

常用步长选择方法：

1. 固定步长（不推荐）：需要经验且难以适应不同迭代。

2. 精确线搜索：求解α使J(θ_k + αδ_k)最小化，计算量大。

3. 非精确线搜索（Armijo准则等）：折中方案，保证充分下降且计算量小。