帶約束的高斯牛頓法求解多音信號分離問題

一、信號模型與優化問題建立

1. 復信號模型

設觀測的復信號由兩個單頻復指數信號加噪聲組成：

$A_0 e^{j(2\pi f_0 n T_s + \phi_0)} + A_1 e^{j(2\pi f_1 n T_s + \phi_1)} + w[n], \quad n=0,1,\dots,N-1$

其中：

$A_0, A_1 > 0$ 為幅度
$f_0, f_1 \in (0, f_s/2)$ 為頻率
$\phi_0, \phi_1 \in [-\pi, \pi]$ 為相位
$T_s = 1/f_s$ 為采樣間隔
$\sim \mathcal{CN}(0, \sigma^2)$ 為復高斯噪聲

2. 參數向量定義

將待估計參數定義為向量形式：

$\boldsymbol{\theta} = \left[ A_0, f_0, \phi_0, A_1, f_1, \phi_1 \right]^T \in \mathbb{R}^6$

3. 優化目標

通過最小二乘準則估計參數：
$\min_{\boldsymbol{\theta}} J(\boldsymbol{\theta}) = \sum_{n=0}^{N-1} \left| x[n] - s[n; \boldsymbol{\theta}] \right|^2$
其中信號模型為：

$\boldsymbol{\theta}] = A_0 e^{j(2\pi f_0 n T_s + \phi_0)} + A_1 e^{j(2\pi f_1 n T_s + \phi_1)}$

4. 約束條件

根據物理意義添加邊界約束：

$\begin{cases} 0 \leq A_k \leq A_{\max} \\ f_{\min} \leq f_k \leq f_{\max} \\ -\pi \leq \phi_k \leq \pi \end{cases}, \quad k=0,1$

二、優化問題求解推導

1. 目標函數的復變函數處理

由于目標函數是復值，需用Wirtinger微積分求導。定義殘差：

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

則目標函數可寫為：

$J(\boldsymbol{\theta}) = \sum_{n=0}^{N-1} r_n(\boldsymbol{\theta}) \overline{r_n(\boldsymbol{\theta})}$

其中 $\overline{(\cdot)}$ 表示復共軛。

2. 梯度計算（Wirtinger導數）

Wirtinger導數定義為：
$\nabla_{\boldsymbol{\theta}} J = 2 \sum_{n=0}^{N-1} \operatorname{Re} \left( \frac{\partial r_n}{\partial \boldsymbol{\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 + \phi_0)}$

$\frac{\partial J}{\partial A_0} = 2 \operatorname{Re} \left( \sum_{n=0}^{N-1} \left[ -e^{-j(2\pi f_0 n T_s + \phi_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 + \phi_0)}$

$\frac{\partial J}{\partial f_0} = 2 \operatorname{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 + \phi_0)} \right] r_n \right)$

相位 $\phi_0$ 的梯度分量：

$\frac{\partial r_n}{\partial \phi_0} = -j A_0 e^{j(2\pi f_0 n T_s + \phi_0)}$

$\frac{\partial J}{\partial \phi_0} = 2 \operatorname{Re} \left( \sum_{n=0}^{N-1} \left[ j A_0 e^{-j(2\pi f_0 n T_s + \phi_0)} \right] r_n \right)$

（ $A_1, f_1, \phi_1$ 的梯度形式類似，只需將下標0改為1）

3. Hessian矩陣近似

為加速收斂，采用Gauss-Newton法近似Hessian：

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

其中雅可比矩陣的行向量為：

$\mathbf{J}_n = \frac{\partial r_n}{\partial \boldsymbol{\theta}} = \left[ \frac{\partial r_n}{\partial A_0}, \frac{\partial r_n}{\partial f_0}, \frac{\partial r_n}{\partial \phi_0}, \frac{\partial r_n}{\partial A_1}, \frac{\partial r_n}{\partial f_1}, \frac{\partial r_n}{\partial \phi_1} \right]$

4. 帶約束的Gauss-Newton算法

迭代格式如下：

步驟1：初始化參數 $\boldsymbol{\theta}^{(0)}$ ，設置迭代索引 $k = 0$

步驟2：計算當前殘差 $r_n^{(k)} = x[n] - s[n; \boldsymbol{\theta}^{(k)}]$

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

步驟4：求解帶約束的線性最小二乘問題：

$\boldsymbol{\delta}^{(k)} = \arg \min_{\boldsymbol{\delta}} \| \mathbf{J}^{(k)} \boldsymbol{\delta} + \mathbf{r}^{(k)} \|^2$

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

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

步驟5：更新參數 $\boldsymbol{\theta}^{(k+1)} = \boldsymbol{\theta}^{(k)} + \mu \boldsymbol{\delta}^{(k)}$ （ $\mu$ 為步長）

步驟6：若 $\|\boldsymbol{\delta}^{(k)}\| < \epsilon$ 則停止，否則 $\leftarrow k+1$ 轉步驟2

三、邊界約束處理策略

采用投影梯度法處理邊界約束：

在每次迭代更新后，檢查參數是否越界：

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

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

$\mathcal{P}_{\Omega}(\theta_i) = \begin{cases} \theta_{\min,i} & \theta_i < \theta_{\min,i} \\ \theta_i & \theta_{\min,i} \leq \theta_i \leq \theta_{\max,i} \\ \theta_{\max,i} & \theta_i > \theta_{\max,i} \end{cases}$

對于相位周期性約束，投影后需歸一化到 $[-\pi, \pi]$ ：
$\phi_{\text{proj}} = \mod(\phi + \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

五、與現有方法的理論對比

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

注：M為粒子數

六、論文書寫建議

問題建模部分：

明確定義復信號模型和參數向量
給出完整的最小二乘目標函數和約束條件

算法推導部分：

詳細說明Wirtinger導數的推導過程
給出Gauss-Newton法的矩陣形式更新公式
解釋邊界投影算子的實現

實驗部分：

比較梯度解析計算與數值微分的精度差異
展示不同SNR下的參數估計Cramér-Rao界

附錄：

提供核心算法的偽代碼
補充投影梯度的收斂性證明

通過以上數學化表達，可避免出現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}
}
其中 $\boldsymbol{\theta} = [A_0, f_0, \phi_0, A_1, f_1, \phi_1]^T$ 是 6維參數向量

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

常用步長選擇方法：

固定步長（不推薦）：需要經驗且難以適應不同迭代。
精確線搜索：求解α使J(θ_k + αδ_k)最小化，計算量大。
非精確線搜索（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