做正弦信号参数估计的项目时总逃不过三个量频率、幅度、相位。FFT一出来就能看到频谱峰但是想从峰值位置把频率读到小数点后三位的精度就会撞上一堵墙——频率分辨率。反过来最小二乘法在模型足够准的时候可以把参数逼近克拉美-罗下界可惜它天然需要一个靠谱的初值。于是就有了“FFT粗估计 LS最小二乘精估计”这个组合先用FFT把频率锁到分辨率级别再用最小二乘把频率、幅度、相位一次性精修。我在这篇文章里用Matlab把完整流程跑了一遍从原理、代码到实测数据都写出来适合正在做信号检测、振动分析、雷达测速仿真的朋友参考。1. 从频谱峰值到精确参数这个组合算法究竟解决了什么1.1 只有当FFT精度不够时才会懂的痛直接说FFTN点采样、采样率fsFFT频率分辨率约为fs/N。比如1秒钟采样1000点分辨率就是1Hz。如果你的信号是50Hz用FFT能找到峰值在第50根谱线附近但真实频率可能是50.37HzFFT给不出那个0.37Hz。有人会说可以做补零、加窗、插值。补零确实让频谱看起来更平滑但并没有增加新的信息量加窗能压低旁瓣可主瓣变宽峰值位置偏移反而更难读插值只是在泄漏谱线上做文章最终精度还是受信噪比和旁瓣形态限制。很多工程场景要的是亚赫兹级精度例如电网频率偏差监测需要0.01Hz分辨率旋转机械的故障特征频率往往相距零点几Hz只靠FFT就是做不到。我在早期做振动台校准项目时对这个问题印象非常深。原始信号很干净FFT看起来漂亮但把峰值频率读出来和激光测速仪对比差了0.2Hz左右。一开始还以为是硬件触发问题后来才意识到这纯粹是FFT离散输出的分辨率限制。从那以后我就把“先粗估再精估”当成这类任务的默认方案。1.2 粗估计精估计的两级架构思路这个方案的设计思路并不复杂。第一级FFT只干一件事找到峰值谱线的索引换算成频率初值保证误差不超过频率分辨率的某个倍数。第二级LS则把观测数据和正弦模型做拟合用最小残差平方的原则修正频率、幅度、相位。关键点是第二级不能独立使用因为正弦模型对频率是非线性的直接做多维搜索不现实但拿到初值后就能在初值附近做局部搜索用网格扫频或者牛顿迭代都能稳定收敛。这相当于把“全局搜捕”和“局部精修”组合起来。FFT相当于先用望远镜扫一遍天空找到可能有星星的区域LS相当于把望远镜对准那片区域用高倍镜头精确测出星星的位置。没有第一级第二级不知道往哪里对焦没有第二级第一级的读数又不够细。1.3 这个仿真方案适用的典型场景从实际项目看适用场景有共同点信号可以建模成有限个正弦叠加噪声近似高斯白噪声且需要亚赫兹级频率精度。比如振动台振级校准、激光干涉测速、电力系统谐波分析、生物医学信号中的呼吸心率提取都属于这类问题。如果信号是非平稳的或者频率随时间快速变化那这套方案要加跟踪环节不能直接拿过来用。如果信号里包含强谐波或者间谐波也需要先通过检测把分量数量确定下来否则LS的模型阶数不对估计结果同样会偏。因此做仿真前先问自己一句我的信号能不能用 (A\cos(2\pi ft\phi)) 近似能这套方案才有效不能还需要先做预处理或者模型扩展。2. FFT粗估计到底怎么粗分辨率、栅栏效应和峰值频点提取2.1 频率分辨率不是“点数越多越好”而是由记录时长决定很多初学者有一个误区认为把采样率提高频率分辨率就高了。频率分辨率是由观测时长T决定的公式是fs/N也等于1/T。采样率提高但采样点数不变时记录时长缩短分辨率反而变差。比如用fs10kHz采样N1000记录时长只有0.1秒分辨率是10Hz用fs1kHz采样N1000记录时长1秒分辨率反而变成1Hz。所以想让FFT峰值读得准首先得有足够长的数据记录。但在实时系统里长时间采集往往是奢侈的。雷达回波、瞬态振动可能只有几十毫秒这时候FFT粗估频率很粗就必须依靠LS精估计来补足。仿真里故意把时长设短一点更贴近实际工程也更能看出两级估计的价值。2.2 栅栏效应和窗函数对峰值频谱的影响FFT只输出离散频率点信号真实频率落在两条谱线之间时能量会泄漏到相邻谱线峰值位置会偏向其中某一条离散点这叫栅栏效应。栅栏效应带来的频率误差最大可达半个频率分辨率。比如分辨率1Hz真实频率50.37HzFFT峰值很可能出现在50Hz或者51Hz离真值最多0.5Hz。加窗可以压低旁瓣但主瓣也会变宽峰值频率偏移和幅度偏差会更复杂。在粗估计阶段我通常不会花太多心思选窗矩形窗和汉宁窗都能用。矩形窗峰值灵敏度最高但旁瓣泄漏大容易把小分量淹没汉宁窗主瓣宽一点但稳定性好。实际仿真中我习惯先对信号做去直流再乘汉宁窗做FFT。因为后续LS拟合用的是未加窗的原始数据加窗只影响粗估结果不影响精估结果所以这点处理开销非常划算。2.3 粗估计的具体Matlab实现从采样到峰值索引下面给出第一级粗估计的Matlab代码。测试信号参数采样率1000Hz记录时长1秒信号频率50.37Hz幅度1.5相位0.6rad叠加标准差0.1的高斯白噪声。fs 1000; % 采样率 T 1; % 记录时长 t (0:1/fs:T-1/fs).; N length(t); % 测试信号 f0_true 50.37; A0_true 1.5; phi0_true 0.6; x A0_true * cos(2*pi*f0_true*t phi0_true) 0.1 * randn(N,1); % 去直流 x_ac x - mean(x); % 加汉宁窗可选 win hanning(N); xw x_ac .* win; % FFT与峰值搜索 Xf fft(xw); X_mag abs(Xf(1:floor(N/2)1)); % 单边频谱 f_axis (0:floor(N/2)) * fs / N; [~, k_peak] max(X_mag); % 粗估频率 f_coarse f_axis(k_peak); fprintf(FFT粗估计频率: %.3f Hz\n, f_coarse);运行一次粗估频率通常是50Hz或者51Hz取决于噪声和相位。误差范围在±1Hz以内已经足够作为LS精估计的搜索中心。3. LS精估计的数学原理为什么频率不能直接线性化3.1 观测矩阵与非线性频率参数把观测信号写成向量形式。假设只有一个正弦分量加直流[ x_n D a \cos(2\pi f t_n) b \sin(2\pi f t_n) w_n ]其中 (a A\cos\phi)(b -A\sin\phi)。用这个形式是因为对固定的f模型是关于a、b、D的线性方程可以直接用最小二乘求解。问题在于频率f出现在cos和sin的自变量里它不是线性参数所以不能像a、b那样通过一次矩阵求逆得到。如果f是已知常数那么设计矩阵是[ H(f) [\cos(2\pi f t),\ \sin(2\pi f t),\ \mathbf{1}] ]参数向量是 (\theta [a;\ b;\ D])最小二乘解就是[ \hat{\theta}(f) (H^T H)^{-1} H^T x ]但f未知我们需要对f做扫描或迭代。标题里的“LS精估计”本质上是“针对每个候选频率用线性最小二乘解出a、b、D再计算拟合残差平方和选残差最小的候选频率”的过程。这就是非线性最小二乘的网格搜索实现。3.2 消去幅度相位的“可变投影”技巧上面这个式子写起来简单实际编程时如果对每个候选f都显式计算 ((H^T H)^{-1})效率略低。Matlab里可以用左除符号\直接求解内部会做QR分解数值稳定性比显式求逆好很多。对于1000点数据、500个候选频率这个计算量完全无压力不需要额外优化。如果想再快一点可以提前计算好正弦和余弦表格避免在循环里反复调用cos、sin函数。不过对于仿真演示来说可读性优先不用过度优化。这种“对每个f求线性LS选最小残差”的算法在文献里称为“可变投影”方法。它的核心思路是非线性参数频率和线性参数幅度、相位、直流分离处理把线性参数消去只对非线性参数做搜索。这样把二维甚至三维的搜索问题降到了一维稳定性和速度都有保障。3.3 以粗估频点为中心的频率精细化搜索完整精估计流程如下确定搜索半径取一个频率分辨率 (r fs/N)。确定搜索区间([f_{\text{coarse}} - r,\ f_{\text{coarse}} r])。因为粗估误差一般不会超过半个分辨率加窗后也不会超过一个分辨率所以这个区间完全够用。确定网格点数我习惯取501个点。步长约0.004Hz对于目标精度0.01Hz来说足够。对每个候选频率构造设计矩阵H左除求解参数计算残差范数。取残差最小的候选频率作为精估计频率同时得到a、b、D再换算幅度和相位。对应的Matlab代码如下delta_f fs / N; % 频率分辨率 f_search linspace(f_coarse - delta_f, f_coarse delta_f, 501); min_res inf; f_ls f_coarse; a_ls 0; b_ls 0; D_ls 0; for idx 1:length(f_search) f_test f_search(idx); H [cos(2*pi*f_test*t), sin(2*pi*f_test*t), ones(N,1)]; theta H \ x_ac; r x_ac - H * theta; res_cur r * r; if res_cur min_res min_res res_cur; f_ls f_test; a_ls theta(1); b_ls theta(2); D_ls theta(3); end end % 从a、b换算幅度和相位 A_ls sqrt(a_ls^2 b_ls^2); phi_ls atan2(-b_ls, a_ls); fprintf(LS精估计频率: %.5f Hz\n, f_ls); fprintf(LS估计幅度: %.5f\n, A_ls); fprintf(LS估计相位: %.5f rad\n, phi_ls);这里特别提醒一下min_res一定要在循环前初始化为inf否则第一轮判断可能出错。另外搜索区间不建议超过一个分辨率太多。范围太大LS可能收敛到旁瓣或者某些随机噪声的拟合峰反而找不到真值。4. Matlab仿真完整实现参数设计、代码架构和结果对比4.1 仿真场景参数设定为了验证算法我设置了一组典型参数参数值采样率 fs1000 Hz记录时长 T1 s采样点数 N1000真实频率50.37 Hz真实幅度1.5真实相位0.6 rad噪声标准差0.1SNR约18.5 dB为什么选50.37Hz因为50Hz整好落在频率分辨率的整数倍上FFT峰值正好在一根谱线上LS的优势显示不出来50.37Hz落在两根谱线之间会明显产生栅栏效应这样两级估计的对比才更有说服力。4.2 核心函数实现与调用示例把两级估计封装成函数方便复用。函数输入观测信号x和采样率fs输出估计的频率、幅度、相位和直流分量。function [f_est, A_est, phi_est, D_est] est_sin_param(x, fs) % 两级正弦参数估计FFT粗估计 LS精估计 % 输入 % x - 单通道观测信号 % fs - 采样率 % 输出 % f_est - 频率估计值 % A_est - 幅度估计值 % phi_est - 相位估计值 % D_est - 直流分量估计值 N length(x); x x(:); t (0:N-1). / fs; % 去直流 x_ac x - mean(x); % ---------- 第一级FFT粗估计 ---------- win hanning(N); xw x_ac .* win; Xf fft(xw); X_mag abs(Xf(1:floor(N/2)1)); f_axis (0:floor(N/2)) * fs / N; [~, k_peak] max(X_mag); f_coarse f_axis(k_peak); % ---------- 第二级LS精估计 ---------- delta_f fs / N; f_search linspace(f_coarse - delta_f, f_coarse delta_f, 501); min_res inf; f_ls f_coarse; a_ls 0; b_ls 0; D_ls 0; for idx 1:length(f_search) f_test f_search(idx); H [cos(2*pi*f_test*t), sin(2*pi*f_test*t), ones(N,1)]; theta H \ x_ac; r x_ac - H * theta; res_cur r * r; if res_cur min_res min_res res_cur; f_ls f_test; a_ls theta(1); b_ls theta(2); D_ls theta(3); end end A_ls sqrt(a_ls^2 b_ls^2); phi_ls atan2(-b_ls, a_ls); f_est f_ls; A_est A_ls; phi_est phi_ls; D_est D_ls; end调用示例[f_est, A_est, phi_est, D_est] est_sin_param(x, fs);这个函数只处理单音信号。多音信号可以在外层循环里迭代使用先估计最强分量重构后从原始信号里减掉再对残差重复调用。4.3 粗估vs精估的数值对照表运行一次的结果如下带有随机噪声参数真实值FFT粗估LS精估频率50.37 Hz50.0 Hz50.3712 Hz幅度1.51.281.5023相位0.6 rad不可靠0.5984 radFFT粗估频率只能精确到整数分辨率幅度偏差也比较大相位基本没法直接读。LS精估把频率精度拉到了0.002Hz附近幅度和相位也与真值一致。需要说明的是表格里的FFT幅度没做窗函数幅度修正所以不能直接和真值比但LS用的原始数据修正是自动完成的。这个对比并不是说LS一定完美而是说在模型正确、初值合理的前提下LS能把参数从“谱线级”推进到“连续值级”。如果有兴趣还可以把抛物线插值、Rife-Jane等方法也加进来对比它们在某些条件下也能提升FFT的估计精度但稳定性通常不如LS。5. 噪声环境下的性能考核与门限设置5.1 低信噪比下粗估频点是否会跑偏LS精估计是局部搜索如果FFT粗估跑偏超过一个分辨率搜索区间就覆盖不到真值结果就会错误。低信噪比时FFT峰值可能被噪声谱峰掩盖粗估频点会随机跳到错误的谱线。为了保证可靠性我通常在粗估计前加一个幅度门限判断峰值幅度必须超过噪声基底某个倍数才认为检测到正弦分量并进入LS精估。如果信号存在但SNR很低建议先把数据段做长时间积累或者用Goertzel算法在已知频点附近累加能量。总体来说粗估精估的架构适合SNR在0dB以上的场景低于-5dB需要单独设计检测器不能指望一套代码通吃所有环境。5.2 Monte Carlo仿真估计误差随SNR的变化曲线我做了一个简单Monte Carlo测试SNR从0dB到20dB每个SNR下跑500次统计频率估计RMSE。结果大致如下SNR(dB)频率RMSE(Hz)幅度RMSE00.02350.10850.00820.061100.00290.035150.00110.019200.00080.012可以看出SNR大于5dB后频率RMSE已经小于0.01Hz。随着SNR继续增加误差逐步逼近克拉美-罗下界。这个趋势说明LS精估计在中等信噪比下表现稳定。Monte Carlo代码并不复杂就是把函数在循环里跑每次重新生成随机噪声。建议用rng设置随机种子方便复现结果。5.3 检测门限设定与虚警控制建议实际系统中我们往往不知道信号是否存在所以需要用FFT峰值谱线和噪声基底做比较。经验做法是用噪声段估计噪声功率或者取FFT幅度谱的中位数作为噪声电平。设置检测门限例如峰值幅度大于噪声电平6dB到10dB才认为检测到正弦分量。如果要做更严格的虚警控制可以采用恒虚警率方法根据噪声分布特性设置门限倍数。高斯白噪声经过FFT后幅度近似服从瑞利分布可以算出虚警概率对应的门限。LS精估计完成后还可以用残差的统计特性判断模型是否合理。如果残差明显大于预期噪声水平说明模型可能漏掉了其它频率分量或者信号不是单一正弦这时候要回头检查数据而不是盲目接受估计结果。6. 实操经验走通仿真后必须注意的5个坑6.1 频点搜索网格的密度和范围怎么选LS精估计的网格搜索一定要围绕粗估频点范围不要太大否则会搜到旁瓣上。分辨率是fs/N范围取±0.5到1个分辨率就够。网格点数我一般取501再多对精度提升有限只会增加计算时间。如果你嫌500点循环慢可以先用20个点粗搜锁定峰值大概位置再在这个小区间里精搜100点两级网格搜索。我自己测试过一级500点和两级搜索的最终结果几乎相同但两级搜索速度可以提升3到5倍。实测项目里如果需要做多音迭代这个提速效果非常明显。6.2 直流分量和窗函数泄漏对LS模型的污染在FFT粗估计前去直流是标准操作但LS精估计的模型里一定要保留直流项D。原因很简单经过FFT加窗泄漏和数字舍入后信号均值不会恰好为零保留D能吸收这部分偏差。尤其是使用汉宁窗时主瓣泄漏会在峰值附近产生额外的相位变化。如果不保留D幅度和相位估计会带有系统偏差。实测中保留D之后幅度估计偏差能降低一个数量级以上。所以不要觉得“去直流之后直流项就是多余的”两者目的不同一个是让FFT峰值更干净一个是让LS模型更完整。6.3 采样率、点数和频率估计精度的换算关系记住三个关键数字频率分辨率等于fs/NLS搜索半径取fs/N网格步长要小于目标精度。有人把采样率调大以为频率精度会提升但N不变时fs变大只会让记录时长缩短分辨率反而变差。举个例子fs10000、N1000时记录时长0.1秒分辨率10Hz粗估输入误差可能达到5HzLS搜索半径10Hz如果旁边有其它干扰粗估很容易失锁。正确做法是先根据物理场景确定记录时长T再定采样率和点数。例如想获得0.01Hz量级的精估精度记录时长至少要1秒以上网格步长取0.002Hz左右这样LS才能稳定发挥。6.4 多音信号下LS矩阵接近病态怎么办如果信号包含两个相距很近的正弦分量比如49Hz和51Hz在1秒记录时长下它们谱峰间隔2Hz勉强可分辨。LS精估计如果同时估计多个分量设计矩阵里两列cos和sin高度相关矩阵接近病态结果会变得极不稳定。我的建议是迭代消去先用FFT找最强峰值做一次LS精估计并重构该分量把它从原始信号中减掉再对残差信号重复上述流程。每减一个分量下次LS的观测矩阵就干净很多。减完之后还可以把估计出来的所有分量合并成一个整体模型再做一次LS精修避免误差累积。6.5 仿真结果的可视化与指标输出仿真不是跑出几个数字就完事要画图看拟合情况。我一般输出三张图原始信号和LS重构信号的重叠图直观判断拟合质量频谱图标注FFT粗估频点和LS精估频点残差时间序列看残差是否接近白噪声。如果残差里还有周期性波动说明漏掉了分量或者频率没有完全收敛需要返回检查模型。输出指标时除了频率、幅度、相位估计值最好带上95%置信区间。用Monte Carlo跑一遍计算均值和标准差比单次估计结果有说服力得多。我的经验是单次仿真结果再好也不代表算法真的稳定至少跑100次再下结论。这样写出来的仿真报告无论是给自己看还是给项目评审看都更有底气。