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

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准则等）：折中方案，保证充分下降且计算量小。

对于大多数应用，内置的 lsqnonlin 配合 Levenberg-Marquardt 算法是最佳选择。只有在研究算法细节或需要特殊修改时，才需要自己实现高斯牛顿法。

代码：

function estimateSignalParameters()% 生成测试信号N = 128;n = (0:N-1)';w_true = [0.3*pi, 0.5*pi]; % 角频率在[0, pi]内c_true = [1.2-0.3i, 0.8+0.5i]; % 复振幅y = c_true(1)*exp(1i*w_true(1)*n) + c_true(2)*exp(1i*w_true(2)*n);y = y + 0.1*(randn(N,1) + 1i*randn(N,1)); % 添加噪声% FFT初始估计Y = fft(y);[~, idx] = sort(abs(Y(1:N/2)), 'descend');f_init = (idx(1:2)-1)/N; % 归一化频率w_init = 2*pi*f_init; % 角频率% 最小二乘估计初始振幅A = [exp(1i*w_init(1)*n), exp(1i*w_init(2)*n)];c_init = A \ y;% 初始参数向量 [ω1, ω2, Re(c1), Im(c1), Re(c2), Im(c2)]theta0 = [w_init(1), w_init(2), real(c_init(1)), imag(c_init(1)), real(c_init(2)), imag(c_init(2))];% 设置边界约束lb = [0, 0, -10, -10, -10, -10]; % 频率下限0，振幅范围-10~10ub = [pi, pi, 10, 10, 10, 10];   % 频率上限pi% 优化选项options = optimoptions('lsqnonlin', ...'Algorithm', 'trust-region-reflective', ...'Display', 'iter', ...'MaxFunctionEvaluations', 3000, ...'FunctionTolerance', 1e-8, ...'StepTolerance', 1e-10);% 运行约束优化[theta_opt, resnorm, residual, exitflag, output] = lsqnonlin(...@(theta) computeResiduals(theta, n, y), ...theta0, lb, ub, options);% 提取优化结果w_est = theta_opt(1:2);c_est = [theta_opt(3)+1i*theta_opt(4), theta_opt(5)+1i*theta_opt(6)];% 显示结果disp('真实参数:');disp(['频率: ', num2str(w_true)]);disp(['振幅: ', num2str(abs(c_true)), ' 相位: ', num2str(angle(c_true))]);disp('估计参数:');disp(['频率: ', num2str(w_est)]);disp(['振幅: ', num2str(abs(c_est)), ' 相位: ', num2str(angle(c_est))]);% 绘图验证s_opt = c_est(1)*exp(1i*w_est(1)*n) + c_est(2)*exp(1i*w_est(2)*n);figure;subplot(2,1,1);plot(n, real(y), 'b', n, real(s_opt), 'r--');title('实部对比'); legend('观测', '模型');subplot(2,1,2);plot(n, imag(y), 'b', n, imag(s_opt), 'r--');title('虚部对比'); legend('观测', '模型');
endfunction residuals = computeResiduals(theta, n, y)w1 = theta(1); w2 = theta(2);a1 = theta(3); b1 = theta(4);a2 = theta(5); b2 = theta(6);s = (a1 + 1i*b1)*exp(1i*w1*n) + (a2 + 1i*b2)*exp(1i*w2*n);r = y - s;residuals = [real(r); imag(r)]; % 实值残差向量
end