简介Lomb-Scargle periodogramLSP是处理非均匀时间序列周期检测的经典工具广泛应用于天文、地球科学及工程监测领域。这份Matlab实现资源包含lombscargle.m主函数与inputtolomb.m示例脚本共2个m文件压缩包仅7KB轻量易用适合需要快速上手LSP算法的科研人员、工程师及学生。函数通过数据预处理、线性谐波建模、最小二乘拟合与功率谱计算等步骤能够针对不等间隔数据准确估计周期信号示例文件则清晰演示了时间戳、观测值和频率范围的输入格式降低了算法应用门槛。已有495人学习下载资源虽小但结构完整可帮助读者从代码层面理解LSP核心原理并直接调用到自身项目中尤其适合用于非等间距采样数据的频率分析、多周期检测及显著性检验等场景是一份实用且高效的算法参考脚本。1. Lomb-Scargle在MATLAB里的落地不均匀采样序列的频谱分析不再玄学做时间序列分析的人八成都在均匀采样数据上用过FFT一旦数据中间缺了几段、观测时刻本身就不均匀FFT的结果就容易掺入大量由采样窗口造成的虚假周期怎么看都不对劲。Lomb-Scargle周期图恰恰是为这类“非均匀采样”场景设计的它直接在原始时间戳上做加权最小二乘拟合把频谱估计从均匀网格的约束里解放出来。这里要拆解的lombscargle.m就是一个MATLAB实现适合做天文光变曲线、振动监测、生理信号、环境传感器数据这类时间戳稀疏且不规整的序列能解决“数据本身不规矩但你还想找周期”的核心诉求。新手可以在半小时内跑通函数调用老手则能在这里找到频率网格参数、归一化方式和显著性检验的边界。2. 为什么FFT会翻车采样窗口、谱泄漏和Lomb-Scargle的数学替代方案2.1 非均匀采样对FFT的破坏力从栅栏效应到掩星伪峰先看一个最典型的场景某课题组拿到一段大气监测数据时间戳不规则平均每15分钟一条记录但偶尔有一个小时的空缺夜间数据稀疏白天加密。直接用FFT分析周期得到的频谱在0.1Hz附近会出现成片的伪峰这些峰的形状和真实物理周期无关而是采样窗口的傅里叶变换叠加上去的。均匀采样的FFT假定时间序列是等间隔的一旦时间戳偏离等间隔信号本身虽然没有变但等效于原信号乘上了一个不均匀的采样窗口函数在频域里就是原始频谱和窗口频谱的卷积。数据越稀疏窗口频谱的旁瓣越高卷积之后真实周期峰会被旁瓣污染严重时直接把主峰淹没。栅栏效应是另一个层面。FFT只能在离散频率网格上取值网格间隔由总时长决定。如果真实周期恰好落在两条谱线之间FFT谱峰幅度会被低估相邻网格的泄漏会让峰形变成“馒头状”。Lomb-Scargle不要求等间隔它的做法是在任意频率点上把数据投影到余弦基函数上通过解析公式算出该频率的谱功率因此频率网格可以由用户单独控制不必受时长倒数的限制。这也是为什么Lomb-Scargle对非均匀序列几乎是默认首选它把“采样不均匀”这件事从错误来源变成可建模的变量。2.2 Lomb-Scargle的数学核心加权最小二乘拟合和时延参数Lomb-Scargle周期图的统计基础是在单一频率下用正弦和余弦函数的线性组合去拟合数据然后看拟合相对残差平方和的改善程度。具体来说对每个测试频率f构造如下模型y(t) A·cos(2πf(t-τ)) B·sin(2πf(t-τ)) ε其中τ是时延参数用来让拟合对时间平移不敏感。当τ取特定值时cos项和sin项的交叉求和项归零使得谱估计在时域平移下保持不变这是Lomb-Scargle对经典周期图的一个重要修正。经典公式里谱功率P(f)与拟合优度的关系可以用归一化的形式写出来P(f) (1 / 2σ²) · [ (Σ(y_i - ȳ)cos2πf(t_i-τ))² / Σcos²2πf(t_i-τ) (Σ(y_i - ȳ)sin2πf(t_i-τ))² / Σsin²2πf(t_i-τ) ]这里σ²是数据的方差ȳ是均值。乘以1/2σ²是为了让谱值在纯噪声情况下服从指数分布这样就能给出明确的显著性阈值。实现lombscargle.m时这个归一化系数是第一个容易出错的地方有的版本除以N有的版本除以2σ²两个版本算出来的峰值尺度完全不同做显著性判断时结论也会跟着变。2.3 频率网格的参数逻辑ofac、hifac和奈奎斯特下限对均匀采样信号奈奎斯特频率是采样率的一半超过这个频率会混叠。非均匀采样没有严格意义上的最高频率因为最短的时间间隔决定了你可以探到多高的频率但稀疏区域的采样会在高频段注入噪声。因此实际实现中常用两个参数控制频率网格过度采样因子ofacoversampling factor和最高频率因子hifachigh-frequency factor。ofac控制频率网格相对于经典基频1/T的密度。T是序列总时长。ofac4表示在0到f_max之间每段基频间隔内再多插3个频点。hifac控制f_max相对平均奈奎斯特频率的位置。f_max一般取hifac·f_nyquistf_nyquist 1/(2·Δt_mean)Δt_mean是平均采样间隔。为什么ofac不能直接给到很高因为频点过多会让相邻频率的谱值高度相关看起来谱峰很光滑但实际分辨率并没有提升。频点多了计算量按线性增长而过采样超过8以后谱峰识别反而容易把同一个峰拆成相邻的两三条细线给后续的峰提取制造麻烦。常见做法是先设ofac4、hifac1做快速扫描找到候选峰之后再在峰附近用ofac16局部加密确认峰位置和宽度。lombscargle.m会把这两个参数当成可选输入并给出默认值这是工程上比较理想的设计。2.4 等价自由度与显著性不是所有峰都值得信拿到一个谱峰值以后下一步是想知道它到底是不是真实周期。Lomb-Scargle的谱值在零假设纯白噪声下服从指数分布所以给定一个谱值P它的误报概率False Alarm Probability可以写成FAP 1 - (1 - exp(-P))^MM是独立频率数近似等于频率网格点数但对过采样网格相邻频点相关实际独立数要按1.5~2倍的有效带宽估算。更稳妥的做法是用蒙特卡洛置换检验把时间戳固定把观测值随机打乱几千次重新计算每个频率的最大谱值看真实数据谱峰落在分布哪个分位点上。lombscargle.m如果只输出功率谱不提供显著性估计那实际使用时需要自己在外部补一个这样的置换模块。检验把戏看后续章节的具体实现。3. lombscargle.m的函数结构与调用从输入参数到谱峰提取的完整实现3.1 函数原型与参数说明在MATLAB里一个典型的lombscargle.m实现会设计成下面的函数签名function [freq, pwr, fap] lombscargle(t, y, ofac, hifac, weights) % LOMBSCARGLE - Lomb-Scargle periodogram for unevenly sampled data % % 输入: % t - 时间戳列向量 (秒/小时/天均可单位一致即可) % y - 观测值列向量与t等长 % ofac - 过采样因子推荐4~8 % hifac - 高频因子推荐1~2取1.0时最高频等于平均奈奎斯特频率 % weights- 权重向量可省略默认全为1 % 输出: % freq - 频率网格 (与y单位对应比如t单位是天则freq单位是1/天) % pwr - 归一化谱值纯噪声下服从指数分布 % fap - 各频点的误报概率估计 (基于指数分布近似)这段代码的逻辑分三层。第一层是参数检查t和y必须是等长的列向量ofac和hifac必须是正数如果缺了weights就用ones(N,1)补齐。第二层构建频率网格最低频从1/(T·ofac)开始最高频到hifac·f_nyquist频率步长是1/(T·ofac)。第三层对每个频点计算时延τ、频谱功率和误报概率。权重参数是多数简化版实现忽略的但对非平稳噪声数据很有用——比如观测误差已知时权重按误差平方倒数取值能显著压制低质量数据对谱估计的影响。3.2 核心计算循环时延、余弦求和和谱值更新下面是lombscargle.m里最核心的计算循环正常运行的核心逻辑是直白的N length(t); t t(:); y y(:); y y - mean(y); variance var(y); % 用分母N-1的方差归一化时需要 if nargin 5 || isempty(weights) weights ones(N, 1); end T max(t) - min(t); freq (1/(T*ofac) : 1/(T*ofac) : hifac*0.5/mean(diff(sort(t)))).; pwr zeros(size(freq)); for k 1:length(freq) f freq(k); tau atan2(sum(weights .* sin(4*pi*f*t)), ... sum(weights .* cos(4*pi*f*t))) / (4*pi*f); % 时延tau使余弦项和正弦项正交 cos_omega cos(2*pi*f*(t - tau)); sin_omega sin(2*pi*f*(t - tau)); cs sum(weights .* y .* cos_omega); ss sum(weights .* y .* sin_omega); cc sum(weights .* cos_omega .* cos_omega); sss sum(weights .* sin_omega .* sin_omega); pwr(k) (cs^2/cc ss^2/sss) / (2*variance); end fap 1 - (1 - exp(-pwr)) .^ length(freq);参数说明集中在时延τ的计算上。这里用atan2处理四倍角的正弦余弦项求和再把角度除以4πf得到的是经典Lomb-Scargle定义里的时延参数。计算cos_omega和sin_omega时带上权重这意味着每个观测点的贡献按质量缩放低权重点在拟合里的影响力被压低。分母里的cc和sss是余弦基函数和正弦基函数各自的加权模方单独平方和再除以模方再除以2倍方差才能保证谱值在零假设下具有指数分布的性质。fap这一行用的是独立频点数近似和精确置换检验比偏乐观但作为前期筛选已经够用。3.3 从谱峰到周期提取峰值检测和置信区间函数输出频率和谱值之后实际业务里要从谱里提取候选周期。常用做法是先定一个显著性门槛比如FAP小于0.01才把谱峰当作候选。MATLAB里的findpeaks配合阈值筛选可以完成这个任务代码示例如下[pks, locs] findpeaks(pwr, MinPeakHeight, 0.0, ... SortStr, descend); candidate_mask fap(locs) 0.01; peak_freqs freq(locs(candidate_mask)); peak_periods 1 ./ peak_freqs; peak_heights pks(candidate_mask);这里findpeaks是MATLAB信号处理工具箱的函数locs返回的是峰在网格里的索引。之所以先用descend排序再按FAP过滤是因为高谱峰和低FAP并不完全等价——在过采样网格里谱峰可能被相邻点抬高但FAP是按全局独立频点数计算峰值稍低但宽度窄的峰反而更可信。得到峰频率之后把频率倒数换成周期时要注意单位匹配t的单位是秒周期就是秒t的单位是天周期就是天不要做额外缩放。3.4 输入数据预处理去趋势、去均值、去离群点Lomb-Scargle模型本身只假设数据是平稳噪声加上若干正弦成分实际观测数据里往往带着线性漂移或趋势项。趋势项能量很大会漏进低频段形成斜坡伪峰。常见的做法是在调用lombscargle.m之前先做一次去趋势方式有很多我一般会用多项式拟合法因为实现最简单trend_degree 1; coef polyfit(t, y, trend_degree); y_detrended y - polyval(coef, t); % 再做一次均值归零 y_detrended y_detrended - mean(y_detrended); % 之后传入lombscargle.m离群点也要小心。如果某个单点观测值比周围数据高出三五个标准差又恰好落在某个周期位的相位上很容易制造出假峰。要么在预处理阶段做中值滤波剔离群要么在调用时给低权重。权重传参设计的价值在这里体现把离群点的权重设成0.1其余正常点设成1比直接把点删掉更安全因为删点会改变时间戳的分布进一步扰乱窗口函数。4. 实战对比Lomb-Scargle、FFT和Welch法的边界与选型4.1 同一组模拟数据下的三种频谱估计为了搞清楚什么时候必须用Lomb-Scargle构建一个模拟场景来说明问题。生成一段数据真实信号是两个正弦波叠加一个周期是12.3小时一个是33.7小时叠加少量噪声。把它按不均匀时间戳采样再分别用FFT、Welch法、Lomb-Scargle分析。模拟代码如下rng(2024); N 500; t sort(rand(N,1)) * 240; % 240小时内随机时间戳 y 1.5 * sin(2*pi*t/12.3) 0.8 * sin(2*pi*t/33.7) ... 0.4 * randn(N,1); figure; subplot(3,1,1); fft_out abs(fft(y(1:256))); f_fft (0:255)/256 * 0.5 / median(diff(t)); % 例如 plot(f_fft, fft_out); title(FFT with first 256 points); subplot(3,1,2); [pxx, f_w] pwelch(y, [], [], [], 1/median(diff(t))); plot(f_w, pxx); title(Welch estimate); subplot(3,1,3); [freq, pwr] lombscargle(t, y, 4, 1); plot(freq, pwr); title(Lomb-Scargle);这个对比里FFT用的是截取的前256个连续点但时间戳并不均匀所以频谱里主峰周围会出现额外的旁瓣结构。Welch法需要均匀重采样或者加窗分段这里用median(diff(t))当作等效采样率高频段会引入混叠误差。Lomb-Scargle直接利用原始时间戳两个正弦峰都能在正确位置出现。这个对比清楚地说明数据时间戳不均匀时FFT的谱是“信号谱与采样窗口谱的卷积”Welch法需要在重采样过程中做出取舍而Lomb-Scargle绕开了这两层问题。4.2 FFT能用的情况等间隔补齐和插值的成本如果数据量非常大比如上百万个时间戳Lomb-Scargle的逐频点循环计算会明显变慢。MATLAB里常见的加速办法是用矩阵运算把多个频率合并到一批计算但即便这样Lomb-Scargle的时间复杂度仍是O(N·M)N是数据点数M是频率网格点数。FFT是O(N log N)在同样数据量下快得多。因此如果数据缺失点很少、时间戳偏差很小可以先用插值补齐成等间隔序列再做FFT或Welch速度优势明显。但插值也有代价。线性插值本质上是对信号做了一个低通滤波高频成分被平滑短周期峰会被压低。样条插值保留更多高频细节但会在缺数据段引入振铃效应产生新的伪周期。工程上的折中做法是缺失比例低于10%时用三次样条插值到等间隔网格缺失比例高于20%直接放弃插值走Lomb-Scargle。缺失比例在这中间时两种方法都试试看结果是否稳定。4.3 窗函数选择对Lomb-Scargle的意义均匀采样做FFT时窗函数设计非常讲究汉宁窗、海明窗各有各的旁瓣特性。Lomb-Scargle不直接加窗但权重向量实际上起到了类似窗的作用。默认全1权重相当于矩形窗旁瓣较高如果知道观测质量随时间段变化比如某台设备夜间漂移严重可以通过权重把夜间数据削弱等效于做了一次时域加权窗。具体做法是把权重设成和观测噪声方差成反比再配合lombscargle.m里的权重参数。这里要注意权重是乘在平方和公式里的不是简单的幅值缩放。如果权重设置不当比如归一化后权重总和过大谱值会被整体抬高误报概率变小假峰更容易混进来。看到谱值整体异常偏大时先查权重是否有归一化。4.4 选型建议表场景条件推荐方法理由时间戳完全均匀无缺失FFT 窗函数速度快谱分辨率高缺失率10%时间戳近似均匀插值补全后FFT速度优势明显缺失率20%或时间戳严重不均Lomb-Scargle避开采样窗口卷积问题需要严格显著性检验Lomb-Scargle 置换检验FAP估计可用指数分布近似数据量大且需要实时处理分块处理 Lomb-Scargle计算量可控避免一次性全算噪声非平稳存在明显质量差异Lomb-Scargle 权重向量权重压低低质量时段影响这张表在选型时很有参考价值优先考虑的是数据本身的结构而不是算法库的便利性。很多初学者一看到非均匀采样就无脑上Lomb-Scargle实际上缺失率很低的情况下插值FFT更快且结果更稳。反过来缺失率很高时强行插值会把谱峰抹平这时候Lomb-Scargle才是正解。5. 避坑配置频率网格、归一化和显著性检验的五个常见问题5.1 ofac调太大真实峰被拆成梳状伪峰现象lombscargle.m输出谱里一个本该很锐利的峰周围出现连续三到四根接近的细峰形状像梳子峰值最大的那根不是真实周期。原因频率网格过度采样相邻频点间隔远小于谱峰固有宽度同一个物理峰被多个频点同时“命中”谱值形成局部多峰结构。这不是数据问题而是网格问题。解决全局扫描时ofac给到4不要超过8。找到候选峰后把频率范围缩窄到峰附近两三倍半宽范围内局部用ofac16或24加密重新计算谱值确认峰位和宽度再取局部峰值。没有局部加密前不要把候选峰直接用freq(locs)当作最终周期值。5.2 谱值大小不对比预期整体高出几个量级现象与文献里相同规模的Lomb-Scargle谱对比谱值整体偏大多数频点都超过10甚至100但峰位没变。原因归一化系数没有按公式用2倍方差处理而是用了总能量或没减均值。数据没去均值就直接算谱低频段会出现大而宽的伪峰且谱值的量纲全乱。解决检查y是否先减过均值。再用var(y, 1)还是var(y, 0)也需要注意。经典公式里的σ²用的是总体方差即分母N但MATLAB的var默认是N-1。如果N只有几十这两种取值差别不小频谱峰值也会差出几个百分点。建议在代码注释里写清楚用的是哪种方差并在文档里注明。5.3 时间单位没统一周期算错一个数量级现象t的单位是小时但代码里把freq直接乘了86400得到周期的单位莫名其妙变成天峰被压到极低频率。原因Lomb-Scargle里频率和周期互为倒数只要你保持t的单位一致频率的单位自动是t的倒数。多此一举的单位换算是最常见的翻车点。解决明确约定输入单位。例如t用小时表示周期结果就是小时t用天周期就是天。不要在函数内部做单位换算把单位换算放在调用前完成让lombscargle.m只处理单一单位体系。5.4 忽略采样窗口效应把宽谱峰当成趋势现象谱的低频端出现一片连续高功率平台峰形宽而平找不出明确的单峰。原因序列总时长内数据分布严重不均匀比如前三分之一密集、后三分之二稀疏等效窗口函数在低频段有较大旁瓣把信号能量摊成平台。常见于设备布设初期采样频繁、后期维护不足的情况。解决先画出时间戳的累计分布观察是否有明显分段。若有对密集段降采样或对稀疏段补权让时间戳在整个时长内更均匀。更稳妥的做法是把时间戳做抖动扰动验证谱峰是否稳定移动再来判断真实周期。5.5 显著性检验只算FAP没有做置换检验现象FAP给出0.001的峰替换掉数据顺序后仍然出现同样位置的高峰说明峰其实是时间戳结构的产物不是信号周期。原因FAP的指数分布近似假设频率网格独立且等权重对非均匀采样窗口不完全成立。时间戳的成簇会让某些频率天然拥有更高的拟合能力谱值偏大但和信号无关。解决不要只信指数近似至少做一轮蒙特卡洛置换检验。把y随机打乱N_perm次1000到2000次即可每次保留原时间戳重新计算全谱记录最大谱值用真实谱峰在置换分布里的百分位作为最终FAP。lombscargle.m只给近似结果时外部补这个检验代码。N_perm 1000; max_pwr_perm zeros(N_perm, 1); for i 1:N_perm y_perm y(randperm(N)); [~, pwr_perm] lombscargle(t, y_perm, ofac, hifac); max_pwr_perm(i) max(pwr_perm); end fap_exact sum(max_pwr_perm max(pwr_actual)) / N_perm;这里randperm是每次把观测值顺序随机打乱时间戳原地不动模拟的是“信号不存在但采样结构相同”的零假设。得到fap_exact以后和fap近似值比较如果差了十倍以上说明指数分布假设在你的时间戳结构下不可靠必须用置换值。6. 时段钟检查把周期扫描变成置信区间验证的三个进阶做法6.1 为什么不直接信任单次谱峰单次Lomb-Scargle得到的峰有时会骗人。数据非平稳、噪声有结构、观测段中有仪器切换都可能导致一次性的伪峰。有一回做某监测项目的数据首次扫描在28.5小时处得到非常强的峰FAP不到0.001看起来板上钉钉。后来把数据按时间切成四段分别做谱分析发现28.5小时峰只出现在第二段其他三段完全没有原因是那段时间传感器出了间歇性故障伪周期来自设备反复跳变。从那以后周期类分析里我拿到任何明显峰都会强制走一遍时段钟检查。6.2 三维时段钟实现频率随时间窗口的变化轨迹时段钟的做法是把整个时间序列划分成若干个滑动窗口每个窗口内单独做谱分析再把窗口位置、频率、谱值三个量放在一起观察。MATLAB实现代码并不复杂window_length 200; % 每个窗口内的数据点数 step_length 20; % 滑动步长窗口重叠率约90% window_starts 1:step_length:(N - window_length); freq_scan []; time_center []; pwr_scan []; for i 1:length(window_starts) idx window_starts(i):window_starts(i)window_length-1; [f_win, p_win] lombscargle(t(idx), y(idx), 6, 1); [pks, locs] findpeaks(p_win, MinPeakHeight, 2.0, ... SortStr, descend); if ~isempty(locs) freq_scan [freq_scan; f_win(locs(1))]; time_center [time_center; mean(t(idx))]; pwr_scan [pwr_scan; pks(1)]; end end % 画三维散点图: time_center × freq_scan × pwr_scan scatter3(time_center, freq_scan, pwr_scan, 36, pwr_scan, filled);窗口长度选200个点是为了保证每个窗口内低频段仍有足够周期数比如你要验证28.5小时周期总时长必须覆盖至少三到四个周期。min peak height设为2.0是因为Lomb-Scargle谱值在纯噪声下均值为1高于2说明谱值比噪声均值高一倍只有这样的峰才值得进入候选。90%的重叠率让窗口之间轨道连续真实周期会表现为一条在频率轴方向平稳的亮带伪峰则会随机散落、无连续性。设置完之后三维图的颜色轴用谱值高亮区一眼就能判断周期是否贯穿全程。6.3 相位折叠图验证拟合残差应该长什么样时段钟确认峰位稳定之后还有最后一个验证步骤把数据按候选周期折叠成相位看波形和正弦模型是否吻合。具体做法是把每个时间戳对主周期取模得到相位按相位排序画散点图叠加上正弦拟合曲线。best_period 28.5; % 来自时段钟确认的周期 phase mod(t / best_period, 1); phase_sorted sort(phase); y_sorted y(phase_sorted); fit_params [ones(N,1), cos(2*pi*phase), sin(2*pi*phase)] \ y; y_fit fit_params(1) fit_params(2)*cos(2*pi*phase) ... fit_params(3)*sin(2*pi*phase); figure; plot(phase, y, k., MarkerSize, 4); hold on; plot(sort(phase), y_fit(sort(phase)), r-, LineWidth, 1.5);如果周期是真的相位折叠图应该呈现一条清晰的正弦曲线散点在曲线附近随机散布且离散步幅与噪声水平一致。如果周期是假的折叠图散点会呈团状分布或出现严重的周期性偏差拟合后的残差还带有显著的高频结构。这个操作不需要额外工具箱基础的矩阵左除就能完成参数拟合工程成本很低。6.4 收在习惯上这三个方法配合起来的价值在于Lomb-Scargle给的是谱估计时段钟给的是稳定性证据相位折叠给的是波形一致性验证。单看任何一个都可能被非平稳性和窗口结构误导三个叠加之后候选周期才算真正落地。现在每次做非均匀采样序列分析我都会强制走一遍先看时间戳累计分布再调ofac和hifac做全谱扫描然后时段钟确认峰位的连续性最后相位折叠看一眼残差四步下来才敢把周期写进报告。这套流程也推荐给正在用lombscargle.m的同行希望帮到你。本文还有配套的精品资源点击获取