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

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

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 ← θ_init

for k = 1 to K do

// 1. 计算当前残差和雅可比矩阵

for n = 0 to N-1 do

r_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|| < ε then

break

end

θ_prev ← θ_new

end

θ_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准则等）：折中方案，保证充分下降且计算量小。
   在线搜索策略上稍微展开一点（例如，“通常初始 μ=1，若目标函数未下降，则通过回溯法（如 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~10
    ub = [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('观测', '模型');
end

function 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

建议： 在正文中简要提及实际实现时通常会使用数值稳定的方法来解法方程，例如：

Cholesky 分解 (如果 H 正定)： 最常用。

QR 分解： 更稳定。

添加正则化 (Levenberg-Marquardt)： (H + λI) Δθ = -g，这在 GN 法不收敛或 H 病态时非常有效。考虑到你的问题有物理约束，LM 可能是更好的选择，或者至少作为备选方案提及。