最近在整理《信号检测与估值理论》的Matlab仿真代码正好把高斯-牛顿迭代法在正弦信号参数估计里的用法重新梳理了一遍。这个场景非常典型一个被噪声污染的正弦信号我们想从采样数据里把幅度、频率、相位估计出来。学这门课的同学肯定有同感——课本里一顿推导CRB和MLE真正写到代码里才发现事情没那么简单。这篇文章把建模、推导、代码和排错完整走一遍如果你正在写类似的仿真作业或者做工程里需要从正弦信号中提取参数应该能直接抄作业。信号模型说白了就是一条公式x(n)A cos(2π f n/Fs φ)w(n)w(n)是零均值高斯白噪声。有人会说这不就是FFT吗FFT确实能给个频率粗估计但它受频率分辨率限制想逼近理论误差下限也就是Cramér-Rao界附近就必须有迭代细化。高斯-牛顿法干的就是这件事。在高斯白噪声下这个迭代结果等价于极大似然估计。换句话说它就是这门课里“估值”部分最落地的算法之一。1. 一切都从带噪正弦信号说起1.1 信号模型与参数估计的目标先定一个统一记号。假设有一个正弦信号被加性高斯白噪声污染采样后得到N点序列x(n)A cos(2π f n/Fsφ)w(n)n0,1,…,N-1。其中A是幅度f是频率φ是初始相位Fs是采样率w(n)服从零均值、方差为σ²的高斯分布。我们的目标就是从x(n)里估计参数矢量θ[A, f, φ]ᵀ。这里假设信号自然是实正弦不考虑复指数形式因为课程里通常从实信号入手工程里大多数采集数据也是实的。这里面有个容易忽略的点参数不是都能分开估的。幅度和相位如果已知频率其实是线性参数因为把cos(2πftφ)展开成a cos(2πft)b sin(2πft)以后对a、b就是线性关系。但频率f藏在三角函数的自变量里面整个模型关于f是非线性的这一下就把问题推到了非线性估计的领域。教科书里那些线性最小二乘、BLUE的闭式解在这里派不上用场必须用迭代数值方法。还要区分“检测”和“估值”的差别。检测解决的是“信号在不在”的问题是一个二元假设检验估值解决的是“如果信号在参数是多少”的问题。这套仿真默认信号已经存在所以我们不做判决直接进入参数估计。如果你后面想加检测环节可以在前面套一个能量检测器或匹配滤波器先把信号找出来再把这套迭代算法丢进去细化参数。1.2 为什么不能用FFT一步到位FFT做频率粗估计是天然的也得是算法第一步。原因很简单对一个有限长正弦序列做FFT峰值位置对应的频率就给出了f的粗略值。但这里有两个让人头疼的限制。第一个是分辨率。对于N点采样FFT频率分辨率是ΔfFs/N。我下面这个仿真里Fs1000 HzN500那么Δf2 Hz。也就是说FFT峰值看起来在50 Hz的话真实频率可能是49.3 Hz也可能是50.7 Hz你根本分不清。学信号处理时我们知道可以补零来加密频谱但补零只是插值没有增加真实的信息量并不能真正突破分辨率极限。第二个是栅栏效应。峰值很可能落在两个相邻频点之间导致粗估频率有系统性偏差。高精度估计需要“超分辨”。这个超分辨不是玄学它靠的是信号模型约束。FFT是拿一堆复指数基函数去拟合信号本身没有利用“这是一个单频正弦、幅度相位都是常数”的先验知识。高斯-牛顿法把所有采样点同时放进代价函数里优化等于在整个观测时间上做相干积累所以频率精度能远高于FFT分辨率。这也是为什么在通信、雷达的频率细测场景里FFT粗估计加迭代细化的组合几乎是标配。2. 高斯-牛顿迭代法原理没那么玄2.1 从非线性最小二乘到迭代公式高斯-牛顿法解决的是非线性最小二乘问题。我们先定义残差r_n(θ)x(n)-A cos(2π f t_nφ)其中t_nn/Fs。代价函数是残差平方和J(θ)Σ r_n²(θ)。在高斯白噪声假设下最小化这个J得到的估计等价于极大似然估计。这个等价关系是信号检测与估值理论里一条重要的桥梁很多同学在课本上看到时没什么感觉直到自己写代码才意识到原来我在做MLE而且我没法直接求导让它等于零只能靠迭代。如果直接用牛顿法求J的最小值迭代公式是这样的θ_{k1}θ_k-H⁻¹g其中g是梯度H是海森矩阵。这个公式理论上收敛很快但海森矩阵计算麻烦里面包含残差的二阶导数项。高斯-牛顿的关键近似是对最小二乘这种特殊结构海森矩阵可以写成H≈2JᵀJ其中J是残差向量对θ的雅可比矩阵。这个近似在残差比较小、模型接近最优解时相当准确。把它代进牛顿法因子2约掉就得到高斯-牛顿迭代式θ_{k1}θ_k-(JᵀJ)⁻¹Jᵀr。如果写成增量形式就是每次解线性方程组(JᵀJ)Δθ-Jᵀr然后令θ更新为θΔθ。这里的J是一个N×3矩阵第n行是r_n对三个参数的偏导。整个过程不需要算二阶导数实现起来非常清爽。从收敛性角度看高斯-牛顿法在初始值接近真值时收敛速度接近二阶也就是迭代几次误差就掉到机器精度以下。代价是它对初值敏感初值太差时会原地振荡甚至发散。这个缺点后面我专门讲怎么处理。2.2 三种实现细节线搜索、阻尼、梯度符号写代码之前先记住三个必须在意的细节。第一个是回溯线搜索第二个是可加阻尼也就是Levenberg-Marquardt第三个是雅可比矩阵的符号。先讲线搜索。直接用高斯-牛顿算出来的Δθ走一步代价函数不一定下降。尤其是初值较差时一步可能迈到山沟外面去。我习惯在每次算出Δθ之后从步长α1开始做回溯如果J(θαΔθ)比J(θ)还大就把α乘以0.5直到代价函数下降。这个操作代码量很小但能把算法的稳定性提升一个量级。阻尼的思路就更进一步。LM方法在JᵀJ的基础上加一个对角矩阵(JᵀJλ diag(JᵀJ))Δθ-Jᵀr。λ大时步长被压缩算法近似梯度下降λ小时退化为高斯-牛顿。工程实现里常用的套路是如果尝试的步长让代价下降就接受并减小λ如果代价不降反升就拒绝这个步长并增大λ。这个自适应机制比固定λ稳健得多。最后是梯度符号。雅可比矩阵的列是残差对参数的导数符号写错的话迭代方向直接反了。我在下面第3节写出表达式和代码你对照着检查。这部分没有捷径只能一个项一个项推导。3. Matlab仿真实战从粗估计到精估计3.1 实验环境和参数设置我用的是MATLAB R2022b其实版本对这个脚本影响不大R2018以后应该都没问题。先设置仿真参数Fs 1000; % 采样率 1000 Hz N 500; % 采样点数 t (0:N-1) / Fs; % 时间列向量N x 1 A0 1.0; % 真实幅度 f0 50; % 真实频率 50 Hz phi0 pi / 6; % 真实相位 30 度 sigma 0.1; % 噪声标准差 rng(2025); % 固定随机种子保证可复现 w sigma * randn(N, 1);% 高斯白噪声 x A0 * cos(2*pi*f0*t phi0) w; % 观测信号这里把t设计成列向量是为了后面雅可比矩阵的维度对齐N是500J必须是500×3。用列向量是从一开始就避免矩阵维度问题。信噪比是A0²/(2σ²)功率算出来是1/0.0250也就是大概17 dB属于中等偏上的情况足够看清算法表现。如果你想让仿真更接近弱信号场景把sigma调到0.5左右就能看到初值和迭代稳性的影响。3.2 初始值怎么选FFT粗估计加线性最小二乘初值对不对直接决定迭代成败。我的套路分两步先用FFT估计频率再固定这个频率用线性最小二乘估计幅度和相位。频率粗估计代码Xf abs(fft(x)); [~, idx] max(Xf(1:floor(N/2)1)); f_fft (idx - 1) * Fs / N;这里只取正频部分因为实信号频谱是共轭对称的。如果峰值出现在零频附近说明信号里直流分量占主导要么先减去均值要么把模型里加一个直流偏置项。幅度和相位的初值我不用“直接取FFT谱线高度”的方式因为谱线高度受窗函数影响误差大。更靠谱的办法是固定f_fft做一个线性回归Xc cos(2*pi*f_fft*t); Xs sin(2*pi*f_fft*t); H [Xc, Xs]; c H \ x; % 线性最小二乘 A_init sqrt(c(1)^2 c(2)^2); phi_init atan2(-c(2), c(1));为什么这样行因为把A cos(2πftφ)展开成c1 cos(2πft)c2 sin(2πft)以后c1A cosφc2-A sinφ它们和x是线性关系。频率固定后这个问题就变成了一个标准的线性最小二乘闭式解就够用。这个初值比直接猜一个A和φ可靠得多基本上后续迭代只需要修正频率和细调相位。3.3 残差与雅可比矩阵的实现先把残差和雅可比放到一个函数里。参数向量theta的三个分量分别是A、f、φ。残差定义是r_nA cos(2πf t_nφ)-x_n还是x_n减模型这个符号只影响梯度的正负不影响最终收敛。我这里用“观测减去模型”的写法也就是rx-model这样雅可比矩阵里就是model对参数的导数取负号。三个偏导分别是∂r/∂A-cos(2πf tφ) ∂r/∂f2πt A sin(2πf tφ) ∂r/∂φA sin(2πf tφ)注意∂r/∂f和∂r/∂φ的符号都是正的因为余弦对f和φ的导数都是负的而残差又先对模型取了个负号负负得正。很多人在这一步搞混你可以自己拿复合函数求导验一遍。函数实现function [r, J] getResidualJacobian(theta, t, x) A theta(1); f theta(2); phi theta(3); phase 2*pi*f*t phi; c cos(phase); s sin(phase); r x - A*c; J zeros(length(t), 3); J(:,1) -c; J(:,2) A * (2*pi*t) .* s; J(:,3) A * s; end这个函数是整个算法最核心的部分。我建议你把偏导表达式和代码一行行对上确认没有笔误。实际调试的时候如果迭代方向总是朝错误方向跑优先怀疑这里不要先怀疑迭代公式。3.4 迭代主循环朴素版和LM版先给朴素高斯-牛顿加回溯线搜索的版本。这个版本我已经在多个数据上试过初值不太离谱的情况下非常稳。theta [A_init; f_fft; phi_init]; maxIter 50; tol 1e-12; costHist zeros(maxIter, 1); for iter 1:maxIter [r, J] getResidualJacobian(theta, t, x); cost sum(r.^2); costHist(iter) cost; delta -(J * J) \ (J * r); % ---- 回溯线搜索 ---- alpha 1; [rNew, ~] getResidualJacobian(theta alpha*delta, t, x); costNew sum(rNew.^2); while alpha 1e-8 costNew cost alpha alpha * 0.5; [rNew, ~] getResidualJacobian(theta alpha*delta, t, x); costNew sum(rNew.^2); end theta theta alpha * delta; if norm(alpha * delta) tol * norm(theta) break; end end A_est theta(1); f_est theta(2); phi_est theta(3);注意循环里每次更新theta后都要重新算残差和雅可比。判断收敛用的是参数增量相对大小而不是绝对大小因为A、f、φ数量级不一样绝对阈值不好定。你也可以用代价函数的相对变化来判断两种都行。LM版本其实只改中间几行。我给出一个简化实现theta [A_init; f_fft; phi_init]; lambda 1; maxIter 100; for iter 1:maxIter [r, J] getResidualJacobian(theta, t, x); cost sum(r.^2); H J * J; g J * r; delta -(H lambda * diag(diag(H))) \ g; [rNew, ~] getResidualJacobian(theta delta, t, x); costNew sum(rNew.^2); if costNew cost theta theta delta; lambda max(lambda * 0.3, 1e-10); else lambda min(lambda * 5, 1e10); end if norm(delta) tol * norm(theta) break; end endLM的lambda初值选1然后根据每次接受还是拒绝来调整。这个方案的好处是几乎不需要调参lambda会自动在梯度下降和高斯-牛顿之间切换。如果你的初值特别差LM往往比朴素高斯-牛顿多跑几十步但不容易炸。3.5 完整函数封装实际仿真中我不会每次把主循环粘来粘去而是封装成一个函数function [theta, costHist, exitFlag] gaussNewtonSineFit(x, Fs, theta0) N length(x); t (0:N-1) / Fs; theta theta0(:); maxIter 100; tol 1e-12; exitFlag 0; costHist zeros(maxIter,1); lambda 1; % LM阻尼参数 for iter 1:maxIter [r, J] getResidualJacobian(theta, t, x); costHist(iter) sum(r.^2); H J * J; g J * r; delta -(H lambda * diag(diag(H))) \ g; [rNew, ~] getResidualJacobian(theta delta, t, x); costNew sum(rNew.^2); if costNew costHist(iter) theta theta delta; lambda max(lambda * 0.3, 1e-10); else lambda min(lambda * 5, 1e10); end if norm(delta) tol * norm(theta) exitFlag 1; costHist costHist(1:iter); break; end end end这样一个函数就能处理单频正弦的估计。注意我把LM和回溯线搜索融合了实际用下来比纯高斯-牛顿加线搜索更省心。4. 仿真实验收敛性、精度与CRB4.1 单次运行迭代轨迹和典型结果我按3.1的参数跑了一次初值大概是f_fft落在50 Hz这个bin上A_init在0.99附近phi_init也能八九不离十。高斯-牛顿迭代的记录大概是这样第一次迭代之后代价函数从一百多掉到接近噪声能量第三次迭代参数变化已经小于1e-6第四次就撞上了收敛阈值。最终频率估计在49.9997左右幅度在0.9981左右相位在0.5238左右非常接近真实值。这里有个值得注意的现象频率的估计精度明显高于FFT的分辨率2 Hz。FFT峰值只能告诉我50 Hz但迭代细化后可以得到0.001 Hz量级的精度。这就是前面说的“超分辨”在数值上的体现。你如果只依赖FFT峰值这个信息是拿不到的。从代价函数的迭代曲线看前两步是断崖式下降后面基本是一条直线。这是典型的二次收敛特征。如果你看到代价函数反复横跳或者下降非常缓慢大概率是雅可比有误或者初值太差导致一直在走梯度下降那一段。4.2 蒙特卡洛实验设计单次运行说明不了统计性能必须跑蒙特卡洛。设计思路是固定A0、f0、phi0改变噪声功率对每个信噪比重复生成数据并估计统计RMSE。SNR_dB -10:5:20; numMC 500; for snrIdx 1:length(SNR_dB) snr 10^(SNR_dB(snrIdx)/10); sigma sqrt(A0^2 / (2*snr)); freqErr zeros(numMC, 1); ampErr zeros(numMC, 1); phaseErr zeros(numMC, 1); for mc 1:numMC w sigma * randn(N, 1); x A0 * cos(2*pi*f0*t phi0) w; % 初值 Xf abs(fft(x)); [~, idx] max(Xf(1:floor(N/2)1)); f_fft (idx-1) * Fs / N; Xc cos(2*pi*f_fft*t); Xs sin(2*pi*f_fft*t); H [Xc, Xs]; c H \ x; theta0 [sqrt(c(1)^2c(2)^2); f_fft; atan2(-c(2), c(1))]; theta gaussNewtonSineFit(x, Fs, theta0); freqErr(mc) theta(2) - f0; ampErr(mc) theta(1) - A0; phaseErr(mc) angle(exp(1j*(theta(3)-phi0))); % 相位差wrap到[-pi,pi) end rmseFreq sqrt(mean(freqErr.^2)); rmseAmp sqrt(mean(ampErr.^2)); rmsePhase sqrt(mean(phaseErr.^2)); fprintf(SNR%6.2f dB, RMSE_f%.4e Hz, RMSE_A%.4e, RMSE_phi%.4e rad\n, ... SNR_dB(snrIdx), rmseFreq, rmseAmp, rmsePhase); end这里特别要提相位误差的处理。求相位估计误差时不能直接theta(3)-phi0因为相位本质上是模2π的周期量。如果真实相位是π/6估计出来是π/62π虽然数值差了2π但其实是同一个值。我把误差转成复数再取angle这样才能把误差wrap回[-π,π)。这个细节统计时容易踩坑。4.3 不同信噪比下的RMSE表现仿真结果总体符合课本预期。高信噪比区域频率RMSE随SNR增加线性下降画成log-log图是一条斜率约-0.5的直线和CRB的斜率一致。幅度估计RMSE也随SNR下降。相位估计在高信噪比时能达到10⁻³ rad量级但SNR低于某个阈值后估计误差会突然增大这是因为初值偶尔跳到错误的谱峰上整个迭代收敛到了局部极小值。这个“阈值效应”在非线性估计里非常经典课本上只会提一句“低信噪比时MLE可能出现阈值现象”但只有自己跑过蒙特卡洛才能直观感受到。你会发现相位估计的RMSE在某些SNR点有一个明显的拐点拐点以下算法变得不可靠拐点以上则完美贴合CRB。这个拐点就是非线性估计的“SNR阈值”。做工程时设计系统的工作点一定要避开这个区域。5. 实战中的坑问题与排错记录5.1 初值不好导致发散高斯-牛顿最经典的问题就是初值差一步就发散。我试过直接把初值设为[1, 30, 0]真值是50 Hz结果迭代几步后频率飞到了几千赫兹代价函数直接变成NaN。原因是初始频率离真值太远残差函数在参数空间里是非凸的迭代一步可能跨越好几个峰。解决办法有三个第一个是坚持用FFT粗估计保证频率初值落在正确的谱峰附近第二个是加LM阻尼或回溯线搜索就算方向不理想步长不敢迈太大第三个是在迭代时限制参数范围比如频率保持0到Fs/2之间幅度保持0到若干倍观测峰值之间。实际工程里三条同时用算法基本不会跑飞。5.2 相位估计跳变是正常的有段时间我发现每跑几百次蒙特卡洛相位误差就出现一个接近2π的离群点。一开始以为是算法出错了后来才意识到相位本来就有2π的周期性。cos(2πftφ2π)和cos(2πftφ)完全一样所以估计出来的相位在某次实验里跳了个2π完全正常不代表算法失败。解决方法是统计前把相位误差做wrap处理就是我在蒙特卡洛代码里写的那样。还有一个方法是把代价函数里的相位残差也做周期修正。但最简单的是在统计结果时处理不要在迭代过程中强行改相位。5.3 雅可比矩阵病态与数值尺度当你把N设得很大比如10万点或者频率很低时JᵀJ的条件数可能变得很大导致线性方程组求解不稳定。原因是时间t的数值很大乘上A和2π以后不同列的梯度数量级差异悬殊矩阵接近奇异。处理方式是对参数重新缩放。比如把t归一化到0到1之间频率也换成“每单位归一化时间的周期数”迭代完再还原成真实Hz。或者更简单一点在解方程前对JᵀJ做一个对角缩放令Ddiag(diag(JᵀJ))然后解(D⁻¹/²(JᵀJ)D⁻¹/²)Δθ-D⁻¹/²Jᵀr最后再把Δθ映射回去。这个技巧叫变量归一化能显著改善条件数。5.4 多分量信号怎么办如果信号不止一个正弦分量比如两个频率很接近的正弦叠加直接套这个单频模型是不行的。工程上我习惯用“逐级剥离”的思路先用FFT看主峰估计第一个分量的参数然后把估计出的分量从观测里减掉再对残差做第二轮FFT和迭代。这个过程叫顺序估计。如果两个频率实在太近低于FFT分辨率顺序剥离也救不了这时候最好把模型扩展成多个正弦叠加一次性估计所有参数。代价是参数维度翻倍初值更难给。比较实用的办法是先做高分辨率谱估计方法比如MUSIC或ESPRIT得到频率初值再丢进多分量高斯-牛顿里细调。5.5 噪声非高斯时的替代思路前面所有推导都假设高斯白噪声。如果实际噪声是拉普拉斯分布或者含有脉冲干扰最小二乘估计不再是MLE而且对离群点非常敏感。这时可以把残差的平方替换成更稳健的函数比如Huber损失迭代结构还是一样只是Jᵀr这一项要乘以权重函数。高斯-牛顿的骨架依然适用变的只是目标函数和梯度。说实话在信号检测与估值理论的课程作业里高斯噪声假设已经够用了。真正到了雷达、声呐的实测数据脉冲噪声和干扰才是常态那时候要做的第一件事不是换算法而是先做数据清洗和去野值。从我自己的实操经验来看高斯-牛顿法在这类问题上最大的价值不是“快”而是它把模型的信息用得很彻底。只要初值别太离谱、雅可比别写错、记住加线搜索或者阻尼它能在十几行代码里给出接近理论极限的估计精度。建议你拿到代码后先跑一遍我给的默认参数再把f0改成49.3Hz、把sigma改成0.5感受一下算法在初值误差和低信噪比下的表现。这比盯着课本上的公式空想要直观得多。如果要继续往下扩展可以试试把这里的高斯-牛顿改成在线递推形式或者把模型升级成复指数信号下一步做卡尔曼滤波跟踪时也用得上。